Numerical simulation of viscous accretion disks


Material Information

Numerical simulation of viscous accretion disks
Physical Description:
xi, 281 leaves : ; 29 cm.
Drimmel, Ronald Eugene, 1964-
Publication Date:


Subjects / Keywords:
Accretion (Astrophysics)   ( lcsh )
Disks (Astrophysics)   ( lcsh )
Astronomy thesis, Ph. D
Dissertations, Academic -- Astronomy -- UF
bibliography   ( marcgt )
non-fiction   ( marcgt )


Thesis (Ph. D.)--University of Florida, 1995.
Includes bibliographical references (leaves 277-280).
Statement of Responsibility:
by Ronald Eugene Drimmel.
General Note:
General Note:

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 002056074
oclc - 33662277
notis - AKP4081
System ID:

This item is only available as the following downloads:

Full Text







For my parents,
for the poor, and to
Sacred Truth and Beauty.


First, I must thank my advisor, Dr. James Hunter, Jr., who has always had high

expectations of me. He has not only given me guidance, but has taught me by example

many of the important qualities that a theorist should possess, including the qualities

of being careful, conservative, hard working, and honest. I thank him as well for his

conversation, and attempts to pass on to me some of his experience and physical intuition.

I also thank him as a teacher, especially for the lesson that one should take the time to

start from fundamental theory, and the importance of deriving even "basic" equations

before applying them to a problem. After my time in graduate school, I am certainly

the richer for our association. I must also thank the other members of my committee,

Drs. H. Kandrup, H. Smith, S. Dermott, and J. Klauder, for their helpful discussions,

questions, and recommendations.

I'd like to thank Chad Davies, officemate extraordinaire, who offered encouragement

both as a fellow scientist and friend. I must also thank Clayton Heller for his helpful

input and questions on numerical aspects of this work, as well as Nikos Hiotelis, who

made available to the Astronomy Department the original N-body code, TREECOD,

and assisted me in getting started. To the many other fellow astronomers and students,

whose conversation I have also benefited from, I give thanks. This work was also made

possible by a NASA GSRP Fellowship, which gave me generous support during most

of this research, which also was supported in part by the University of Florida and the

IBM Corp., through their Research Computing Initiative at the Northeast Regional Data

Center of Florida, without which the numerical work would not have been possible.

On a more personal and broader note, I would like to thank Jonathan Potter for his

friendship, and keeping me laughing when I needed it most, as well as all my other

friends, who have assisted me more than they will know. This includes Orlando Espin,

for teaching me that my love for astronomy and my love for the poor were not in

contradiction, but that I can serve all God's children with and through my work as an

astronomer. Also on a spiritual note, and another friend who needs mention, is Mary

Flanagan, who gently guided me on my way. I have been blessed as well by the entire

community of Saint Augustine Catholic Church, which has provided sustaining strength,

healing, and abundant joy. I give gracious thanks as well to Teilhard de Chardin, for

helping me see the big picture. To my parents, who have always believed in me, I am

most grateful. Lastly, I must give ultimate thanks to my God, Master of the Universe and

Lover of Souls, who created all in beauty, and who blessed me with the eyes to see it.



LIST OF TABLES ................

. . . iii

. vii

LIST OF FIGURES ..........................

ABSTRA CT .......................... ... .. ... ... .... xi


1. INTRODUCTION .................................... 1

2. METHODS ........ .... ................

Development of the Hydrodynamic Code ........
N-body Seed Code ...................
S P H . . . .
Variable Smoothing Length and Neighbor Searching
Hydrodynamic Equations ................
Integration Method ...................
Gravitational Softening .................
A ccretion .. ... ... ... .. .. ..
Summary of Code Parameters .............

Tests . . . .
M aclaurin Disk ................
Two-Dimensional Shocks ..........

Analysis of Disks ..................
Evaluation of Effective Shear Viscosity .
Evaluation of Modes and Their Frequencies

. ....... 6
. . 6
. .. . 6
. . 6
. . 10
. . .. 14
. . 16
. . 19
. . 20
. . 23
. . 25

. . . 26
. . . 26
. . . 27

. . . 30
. . . 30
. . . 33


Two-dimensional Disks ...............
M aclaurin Disk .................
Exponential Disks ...............
Inner Disks ....... .............


Param eters ......................

Evolution of Disks .................
General Features ................
Modal Evolution ................
Accretion ....................
Tests of Code Parameters .............

. . . 56

. . . 56

. . 60
. . 60
. . 62
. . 65
. . . 7 1

Orbital Evolution ..............
Formation Conditions ............
Encounter Model ..............

R results . . .
Code Development ..........
Modal Evolution ............
Accretion ....... .......
Formation of Satellites ........
Future W ork .................

A.T2DSPH ....................

B. Algorithm for T2DSPH ............

C.PROGRAM EXPD ..............

BIBLIOGRAPHY ................

. . 84
. . 85
. . 90
. . 95

. . 108
. . 108
. . 108
. . 111
. . 112
. . 113
. . 115

. . 119

. . 262

. . 263

. . 277


. 281


Code Parameters ...............

M odels . . .

Models at T=28 dynamical times ......

Accretion Rates (x 10-3) ..........

Normalized density amplitude growth rates

Satellites ........... ........

........and saturation levels...

and saturation levels.

. . .

Table 2-1:

Table 4-1:

Table 4-2:

Table 4-3:

Table 4-4:

Table 5-1:


Figure 2-1:

Figure 2-2:

Figure 2-3:

Figure 2-4:

Figure 2-5:

Figure 2-6:

Figure 2-7:

Illustration of a two-dimensional hierarchical tree cell structure with 17
particles. There are two levels of subcells below the root cell. 7

Surface density and velocity profile of the Maclaurin disk after 26
dynamical times.............. .......... ....... 27

Surface density of three two-dimensional accretion shocks. The bottom
accretion shock's surface density is offset from the top by -2, and the
middle by -1, from the top accretion shock. . ... 29

Normalized power spectrum of modes of model A2 (see Chapter 4) at
Time = 21 dynamical times. ....................... 35

Dynamic power spectrum of the normalized density for mode = 2 of
m odel A 2. . . . .. 36

Maximum of the power spectrum, in modes = 1 4, of the normalized
density in model A2. ........................... 37

The average Lombe periodogram of normalized density temporal
fluctuations along Cartesian axes for model A2 . 38

Figure 3-1: Initial distribution of particles for the Maclaurin Disk. ... 41

Figure 3-2:

The initial Toomre Q parameter (solid line) and tangential velocity
(dotted line) with respect to radius for a Mc/MD = 3 disk. .. 46

Figure 3-3: Initial distribution of particles in an exponential disk. . 48

Figure 3-4:

Figure 4-1:

Initial density profile of an exponential disk with MD = 0.25 and a scale
length of r, = 0.25. Each point corresponds to the estimated density at
each particle position, evaluated with the SPH method, while the solid
line is the intended density profile. . . .. 49

Particle distribution of Mc/MD = 3 models at Time = 10.0 dynamical
tim es. . . . . 73

Figure 4-2: Particle distribution of Mc/MD = 3 models at Time = 29.5 dynamical
tim es. . .. . . 74

Figure 4-3: Particle distribution of Mc/MD = 1 models at Time = 10.0 dynamical
times. ......... ..................... 75

Figure 4-4: Particle distribution of Mc/MD = 1 models at Time = 30.0 dynamical
tim es. . .. . . . 76

Figure 4-5: Maximum power in modes 1 through 4 for model B2. ... 77

Figure 4-6: Resonances in an accretion disk with Mc/MD = 3. Solid lines correspond
to the corotation resonance, the dotted lines to the Inner Lindblad
resonance, and the dashed lines to the Outer Lindblad resonance. 78

Figure 4-7: Mass accretion, as a fraction of initial disk mass, for models Al through
A 3. . . . . 79

Figure 4-8: Mass accretion, as a fraction of initial disk mass, for B models. .. 79

Figure 4-9: Mass accretion, as a fraction of initial disk mass, for C models. .. 80

Figure 4-10: Mass accretion, as a fraction of initial disk mass, for D models. 80

Figure 4-11: Mass accretion in models A2, A4, and A5, as a fraction of the initial
disk m ass. . . . . ... 81

Figure 4-12: Constant mass accretion rates for Mc/MD = 3 accretion disks 82

Figure 4-13: Constant mass accretion rates for Mc/MD = 1 accretion disks 82

Figure 4-14: Mass accretion of the encounter model (solid line) compared to the
mass accretion of model D2 (dotted line). . ... 83

Figure 4-15: Mass accretion, as a fraction of initial disk mass, for model D2 is
shown with three test models. ........................ 83

Figure 5-1: Evolutionary sequence of model D2 showing formation of satellite.. 96


Figure 5-2: Evolutionary sequence of model D2 continued. . ... 97

Figure 5-3: Final configuration of model D2, showing the position of the satellite at
earlier times. ............. .................. 98

Figure 5-4: Radial density profile and particle distribution of satellite D2-1. .99

Figure 5-5: Final configuration of model B2, showing the position of the satellites at
earlier tim es .. ... .. .. .. .. .. .. .. .. 100

Figure 5-6: Final configuration of model B3, showing the position of the satellites at
earlier tim es .. .. . . 101

Figure 5-7: Final configuration of model Bl, showing the position of the satellites at
earlier tim es .. .. .. .. .. .. .. ... .. ... .. 102

Figure 5-8: Final configuration of model DI, showing the positions of the satellites
at earlier times. Satellite 4, which has been reabsorbed by the disk, is
not show n. .. .. ... .. .. .. .. .. .. .. ... 103

Figure 5-9: Encounter model ................................ 104

Figure 5-10: Encounter model continued. . . ... 105

Figure 5-11: Encounter model continued. . . ... 106

Figure 5-12: Encounter model at Time = 16.0. Note the three satellites that have
formed in the tidal tail ............................ 107

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy




August, 1995

Chairman: James H. Hunter, Jr.
Major Department: Astronomy

A two- and three-dimensional numerical hydrodynamics FORTRAN code is devel-

oped to model viscous accretion disks, employing the method of smoothed particle hy-

drodynamics. The effective shear viscosity present in the code is evaluated. Using a

polytropic equation of state, models of self-gravitating accretion disks are evolved with

central mass to disk ratios of I and 3, with ratios of specific heats of 2 and 5/3, and

with an artificial viscosity parameter of 1, 0.5, 0.25. From these models it is found that

a characteristic mass accretion rate, constant with time, is maintained by the accretion

disks, and that mass accretion is inversely proportional to the strength of the local shear

viscous force. This is a consequence of the effect of local viscous forces upon the global

non-axisymmetric modes that primarily drive accretion in these disks. In addition, the

formation of satellites is observed in models which also develop a dominating m=l mode.

The relevance of these satellites to planet formation is speculated.


Accretion disks form an important class of astrophysical objects, occurring on many

scales in a variety of contexts. In quasars and AGNs they are believed to be the central

driving engine responsible for the extreme luminosities generated by these objects. In

a more common manifestation, they are believed to be a natural product of the star

formation process (Shu et al., 1987), and to be the progenitors of potential planetary

systems. As a structure that preprocesses infalling material before it accretes onto the.

central protostar, the accretion disk is referred to as a protostellar disk. Later, when mass

accretion has virtually ended and the disk mass is small compared to the young star at its

center, it is designated as a protoplanetary disk. This work, though it addresses accretion

disks in general, is focused on accretion disks with the general characteristics of massive

protostellar disks: the self-gravity of the disk will be important, and the central object,

representing the protostar, will be small compared with the spatial dimensions of the disk.

Though the existence of protostellar disks was first theoretically argued as a natural

consequence of the conservation of angular momentum, the observational evidence for

these objects has greatly improved. Initially the observational evidence was indirect.

Infrared (IR) excesses were observed around young stellar objects, particularly T Tauri

stars, though the star itself was unobscured by the material radiating at lower temperatures

(Rydgren & Zak, 1987). Ultraviolet (UV) emission lines are associated with these objects

as well. Other indirect evidence is the angular momentum regulation of T Tauri stars


with IR excesses (Edwards et al., 1993). The most plausible, and simple, of explanations

for these observations was that the material responsible for the excess IR emissions was

in a flat, disk-like structure, and that material accreting from the disk onto the star was

producing the UV emission. Many of these inferred disks have little mass, and are still

thought to be examples of disks at a later stage in evolution.

More massive disks were inferred from young stellar objects with much larger IR

excesses (Hillenbrand et al., 1992). The "flat-top" IR spectra of these objects were

initially explained by massive disks that were generating their own luminosity through

viscous processes (Adams et al., 1988). However, more recently the spectra of these

objects were shown as more probably resulting from a dusty, infalling envelope (Whitney

& Hartmann, 1993; Hartmann, 1993). Hence, to this point, the indirect evidence for

large gaseous disks around young stars has provided enticing, but not decisive, evidence.

However, with new and improved instrumentation, such as the JCMT-CSO submillimeter

interferometer and the serviced Hubble Space Telescope, direct evidence of gaseous disks

around young stellar objects has been observed. (Koerner et al., 1993; Lay et al., 1993)

Numerically many have shown, as per expectations (Lin & Pringle, 1990), that

massive disks do form from a collapsing cloud possessing angular momentum equivalent

to that observed in molecular cloud cores (Bodenheimer et al., 1990; York et al., 1993;

York et al., 1995). In such a self-gravitating disk, nonaxisymmetric modes play a crucial

role, as they effectively transport angular momentum, thereby driving the evolution of the

disk. In addition, ubiquitous shear viscous forces are an important transport mechanism in

a gaseous disk; nonaxisymmetric modes may dominate when present, but viscous forces


will always be present to some degree. In general the resulting angular momentum

transport will be outwards, causing the disk to become more extended, while the bulk

of the mass flows inward, resulting in accretion onto the central object (Lynden-Bell &

Pringle, 1974). The strength of the viscous force has been inferred from the lifetimes of

the disks themselves, and is often termed as anomalous viscosity, as the strength of the

viscous forces are much stronger than can be explained by the kinetic molecular viscosity

of the gas itself. The most common type of viscosity that has been suggested is turbulent

viscosity, though several mechanisms for inducing the turbulence have been suggested.

These mechanisms include convection induced by thermal gradients, which requires an

optically thick disk, and a magneto-rotational instability, which requires the presence of

a weak magnetic field (Balbus & Hawley, 1991; Hawley & Balbus, 1991). In this work

I do not attempt to shed light on the mechanism responsible for the turbulence, but will

nonetheless assume that turbulence is responsible for the shear viscosity in the accretion

disk. These two mechanisms of angular momentum transfer, the global nonaxisymmetric

modes and the local viscous processes, do not act independently, but affect one another.

The theory of accretion disks, motivated in large degree by their occurrence in some

binary systems, particularly cataclysmic variables, has been most fully developed for

the thin Keplerian disk (Pringle, 1981). In this regime the self-gravity of the disk is

neglected, and only pressure supports the vertical structure. In contrast, the theory of

self-gravitating accretion disks is less developed. This is in part due to the nonlinearity

of the problem, which dictates that numerical studies are necessary to understand these

disks. Early numerical work on self-gravitating accretion disks began with N-body


modeling (Cassen et al. 1981; Tomley et al., 1991). Papaloizou and Savonije (1991)

analytically investigate the instabilities that exist in two-dimensional self-gravitating

disks, and numerically follow the evolution of these unstable modes. Their numerical

method was to solve the hydrodynamic equations on a polar grid, which necessitated

the unphysical assumption of rigid inner and outer boundary conditions. They confirm

the importance of the nonaxisymmetric modes in redistributing the mass of the bounded

disk, though the subsequent evolution they observe characterizes real disks only in a broad

sense. In addition, their numerical code has an indeterminate amount of shear viscosity,

which also contributes to the angular momentum transport and mass redistribution that

they observe in their models.

Only recently have self-gravitating three-dimensional disks been investigated through

numerical simulation. The stability of three-dimensional tori has been investigated

both analytically and numerically (Papaloizou & Pringle, 1984; Papaloizou & Pringle,

1985; Zurek & Benz, 1986; Tohline & Woodward, 1992). Most recently Laughlin and

Bodenheimer (1994) have investigated the initial evolution of a disk that is formed from

a previous calculation of a collapsing gas cloud. In this later work ring modes that

had formed in a two dimensional collapse calculation were shown to be unstable in a

three-dimensional hydrodynamic simulation.

In this work I extend the numerical work of Papaloizou and Savonije (1991) by

providing a model that more closely resembles that of an astrophysical accretion disk,

but that is also restricted to two dimensions. In particular, the hydrodynamic evolution

of the disk is followed with a smoothed particle hydrodynamics (SPH) code, which


approximates a gas with a finite number of particles (Benz, 1990; Monaghan, 1992).

The evolution in phase space of these particles is determined by the Lagrangian form

of the hydrodynamic equations. As this method is not restricted to a fixed grid with

rigid boundaries, the disk is allowed to expand naturally in radius, and an inner boundary

condition allows accretion onto the central gravitating object. For the sake of continuity

with Papaloizou and Savonije (1991), the disks are initially given an exponential density

profile, though unlike them the disks are perturbed by introducing noise in the density

profile. In this way all unstable modes will be excited. The primary disk parameters

which are varied are the star to disk mass ratio, the ratio of specific heats, and the effective

shear viscosity of the disks. Like Papaloizou and Savonije (1991), a polytropic equation

of state is used to describe the gas. In particular, the effect of the shear viscosity on the

evolution of self-gravitating disks is investigated.

The following chapter describes the numerical method employed to model accretion

disks and its implementation in a FORTRAN code named T2DSPH. A listing of the

code itself is provided in Appendix 1. Tests of the code are also described, including

an evaluation of the effective shear viscosity that is present. Chapter 3 gives the details

of the initialization of the models that are evolved, an important and nontrivial aspect

of SPH. In Chapter 4 the initial parameters of the accretion disk models are first briefly

described, followed by a presentation of the results of the modeling. The evolution of the

nonaxisymmetric modes and the mass accretion observed in the models are discussed.

Chapter 5 addresses the formation of satellites that occur in some of the models, and

Chapter 6 summarizes the results and speculates on their significance.


Development of the Hydrodynamic Code

In this chapter the development of the numerical hydrodynamic code, named

T2DSPH, is described. A listing of the code is given in Appendix A, and an algo-

rithm of the code is given in Appendix B. Test simulations are also described, including

the evaluation of the effective shear viscosity present in the method. In addition, methods

for analysing the disks are discussed in the last section.

N-body Seed Code

The code developed to evolve accretion disks is itself evolved from a version

of Hernquist's N-body code called TREECOD (1987). This code uses a method of

calculating the gravitational forces on a system of particles, due to the collective influence

of the entire system of particles, using a method termed the hierarchical tree method. At

any particular time a system of N particles is spatially described, in Cartesian coordinates,

by a set of cells hierarchically organized. The root cell encompasses the entire system;

then succeeding levels of subcells are created. In two dimensions each cell will have four

subcells, while in three dimensions each cell has eight subcells. This is accomplished

by simply dividing a cell in half in each dimension (see Figure 2-1). The number of

particles in each cell is found, and successive levels of cells are created until subcells

have either one or no particles in it. Cells with a single particle can be considered the

leaves of the tree. The spatial aspect of the tree is expressed by allocating to memory the

position and dimension of each cell. Hence, the tree structure completely describes the

spatial structure of a system of particles. Further, by recording for each cell the pointer

to its parent cell, this tree may be ascended or descended.

Figure 2-1: Illustration of a two-dimensional hierarchical tree cell structure
with 17 particles. There are two levels of subcells below the root cell.

Moving from leaves to root, the mass and the center of mass of each cell is quickly

determined; in describing the spatial distribution of a system of N massive particles, the

mass distribution of the system is also described. This hierarchical tree structure can

now be utilized to calculate the gravitational forces on any single particle, and allows the

implementation of an advantageous approximation: the number of gravitation terms is

made significantly less than the N-1 terms a direct sum over all the other particles would

yield, by treating cells as gravitating particles. In other words, a single gravitational term

due to a given cell, its mass taken to be at the position of its center of mass, is used to
U U *

due to a given cell, its mass taken to be at the position of its center of mass, is used to

represent the collective gravitational influence of all the particles within the cell. Using

this approach, the gravitational force of a system of N particles on one of its members

can be approximated with a sum of Nt terms due to a set of Nt cells, where Nt << N.

The choice of cells used in the gravitational force calculation will in general be

different for each particle, and the choice of which cells to use will determine the accuracy

of the approximation. To determine the set of cells used for a given particle, the tree is

descended from root to leaves, with the descent continuing to the next level of subcells

until the cell subtends an angle, as "viewed" from at the particle, which is less than

a specified tolerance angle, 0. In effect this method will treat the influence of nearby

neighbors as a direct sum, but simplify the influence of particles at a distance by grouping

them into larger cells. This approximation is then based on the philosophy that the local

gravitational field is not sensitive to the detailed spatial distribution of particles at a

distance. The accuracy of the approximation is increased further by including higher

order terms in a multiple expansion of the mass distribution within the cells. Taking

the mass to be at the center of mass, as described above, we already have the monopole

expansion term. The next higher order correction term is the quadropole term; since the

center of mass is used as the expansion center the dipole term is zero.

Hernquist has found that for a system with 4096 particles that a tolerance angle of

0 < 0.8 radians gives a relative error in the ith component of the acceleration, A(6ai)/ai,

of less than 1%, where he defines the mean absolute deviation from the mean error,


and the absolute average acceleration,

_- arect (2.2)

He define the mean error as

6a= baij, (2.3)

where 6aij = atee adirect is the difference, in the ith component of the acceleration,

as calculated by the tree hierarchical and direct summation methods for particle j. This

was measured by Hernquist (1987) for a Plummer sphere, which has a density profile of

3M r2
p(r) = 0 (2.4)
4r (r2 + r5

where M is the mass of the system and ro is the scale length. For a typical two-

dimensional exponential disk, I measured an average relative error in the accelerations,

defined as
(ba oai)/arec (2.5)
For both an axisymmetric and a nonaxisymmetric disk, with a tolerance angle 0 = 0.8,

I found the average relative error to be of the order of 1%.

The computational time required to sample the gravitational force field of an N-body

system at N points is proportional to the number of gravitational terms to be calculated.

Therefore the computational time of the direct sum method is of the order of N2, whereas

Hernquist finds that the tree hierarchical system is of the order of Nlog(N). Models of

accretion disks presented in this work employ approximately 5000 to 10000 particles;


for these numbers of particles the increase in efficiency, as compared to using a direct

summation method, is over a thousand fold. Nevertheless, there are other methods

for estimating the gravitational field of an N-body system that are more computationally

efficient, such as grid methods that use fast Fourier techniques to evaluate the gravitational

field on a set of grid points. The gravitational force upon a single particle is then

interpolated from the nearest points. While such methods are faster, the tree hierarchical

method possesses several attractive advantages. Grid methods are limited in resolution

to the grid size, whereas the tree method's local resolution is limited by the local particle

distribution itself, or the gravitational softening parameter, since its approximation is only

applied to particles at a distance. In addition, tree methods are not limited to the size or

geometry of a grid, as the root cell is resized at every step and empty cells do not take up

memory. Hence tree hierarchical methods conserve mass exactly. These characteristics

make tree methods ideal for modeling violent encounters, as well as models that have

a broad range in particle densities.


Using this N-body code as a skeleton, a two-dimensional hydrodynamic code was

developed, called T2DSPH (Appendix A). The numerical hydrodynamic method chosen

to model accretion disks is known as smoothed particle hydrodynamics (SPH), and in

developing the code I followed much of the guidance given by Hernquist and Katz's

own combination of SPH and the tree hierarchical method (1989). A three-dimensional

version of this hydrodynamics code was then developed and dubbed T3DSPH.

In essence, "SPH is an interpolation method that allows any [continuous] function to

be expressed in terms of its value at a set of disordered points the particles," (Monaghan

1992). A finite number of particles is used to describe the density and bulk velocity field

of a fluid; hydrodynamical quantities needed to calculate the acceleration at the points

are estimated; and the accelerations are then applied to the particles using an appropriate

integrator. SPH, then, is a Lagrangian technique of solving the hydrodynamical equations,

employing a finite number of particles to approximate the fluid. A great advantage of such

a method is that it is quite general, being easily amendable to any choice of equation

of state. Also, a Lagrangian approach has an advantage over an Eulerian technique,

which must employ a grid, in that no symmetry is imposed or assumed, and further, the

simulation is not confined to a box of finite size. Reviews of this technique, which has

been popular in recent years to model a variety of astrophysical problems, are provided by

Monaghan (1992) and Benz (1990). The following discussion of the mathematical basis

of SPH follows the discussion given both in Benz (1990), and Hernquist and Katz (1989).

In SPH, quantities are estimated using the smoothed estimate,

(f(r)) = f(r') W(r r', h)dr', (2.6)

where W is a smoothing kernel and h a smoothing length. The kernel W is spherically

symmetric and normalized:

IW(r, h)dr = 1 (2.7)

If the kernel W is a function strongly peaked about r = 0 then the estimate (f) can be

written as a Taylor expansion:

(f(r)) = f(r) + c(V2 f)h2 + 0(h3) (2.8)

where the coefficient c = f u2h3W(u)du, (u = Ir r'|/h) is independent of h. The

term proportional to h has vanished because W is an even function in r. Hence it is said

that the smoothed estimate of the function f is second order accurate in h. If the quantity

f is known on a finite set of points described by a point distribution,
n(r) = 6(r- rj), (2.9)
then the integral expression for the estimate (f) (equation 2.3) can be written as a sum:
(f (ri)) = mjf(r) W,. (2.10)
j p(r,)

using the relation

(A(r) (A(r)) O(h2)
B(r) (B(r)) + (2.11)

I have introduced the abbreviation Wij = W(Iri rj h), and written (nj) in the form

p(rj)/mj, where mj is the particle mass and p(rj) is the density at the particle j. Using

equation (2.7), the density is estimated using
p(rj) = Zmj Wu. (2.12)

Because of the form of the estimate [equation (2.7)], it is most convenient to estimate

quantities with the general form pA:
(pA)i = WmAjWi. (2.13)

This convention is considered a "golden rule of SPH" implementation by Monaghan


Various smoothing kernels can be used; in this work the spherically symmetric kernel

r1 6u' + 6u3, if 0 < u < 1/2
W(r,h) = 2(1 u), if 1/2 < u < 1 (2.14)
S0, otherwise,
is used, where u = r/h. This kernel is a fourth order basis spline, normalized for three

dimensions. In two dimensions normalization yields the leading factor of 10/77rh2. I

should also note that the traditional way of writing the above kernel is with h --+ 2h,

so that W is nonzero for r < 2h rather than h, but the above form is more convenient

numerically. The kernel above is the one most widely used in the SPH community, and

has the advantage that it defines a local interpolation since, being nonzero only when

rij = |ri rjI < h, neighboring particles alone are used to make the estimate. Higher

order basis splines were tested but did not yield significantly better estimates of a given

function for equivalent smoothing lengths, and gave worse results when derivatives were


One of the powerful features of SPH is the ease with which the gradients of quantities

can be estimated:

(Vfi) = w(2.15)

(For the sake of brevity I will in general be writing fi for f(ri), and sums are to be

understood as being over neighboring particles only.) Note that the kernel is symmetric

(Wij = Wji), while its gradient is antisymmetric (Vr, Wij = -Vr, Wji). The divergence

of a vector quantity, such as velocity, can take a similar form. However, if the estimate

of a divergence is to be used in an equation of motion, it must be written in a form

that is antisymmetric in i and j in order to conserve angular momentum when particles

interact. Thus, for the quantity V v, the identity pV v = V pv v Vp is used to

derive the following estimate:

(V.v) = m,vij VWij, (2.16)

where vij = vi vj (see Monaghan 1992).

Variable Smoothing Length and Neighbor Searching

So far I have given the SPH formalism in the case where the smoothing length h

is a constant. In general, however, particles are given individual smoothing lengths hi,

which are adjusted, both in space and time, so that the SPH particles have approximately

the same number of neighbors. This refinement will insure that quantities throughout

a model are estimated with comparable accuracy, and enables SPH to adjust the local

resolution as the local number density changes during the course of a model's evolution;

denser regions will have finer resolution than diffuse regions. However, the form of the

kernel Wij, with its dependence on h, is now uncertain. Again, to insure that angular

momentum will be conserved, the kernel is rendered symmetric in i and j by defining

the smoothing length h = hij = (hi + hj)/2. Neighboring particles j will now be

used that satisfy the condition rij < hij. Another way to symmetrize the kernel is to

define Wij = (W(rij, hi) + W(rij,hj))/2, which is the form used by Hernquist and

Katz (1989).

Determining the set of neighboring particles for each particle (those satisfying the
condition rij < hij) is known as "neighbor searching", and can be one of the more time
intensive tasks in an SPH program. In the version of SPH developed for this work, the
neighbor lists are efficiently compiled using the already existing hierarchical tree data

structure used in the N-body portion of the code (see the beginning of this chapter).
First, the list of neighbors within hi is found by descending the tree, starting at the root,
and continuing to sublevels of those cells that intersect the neighborhood of particle i,
defined by hi. When the descent reaches a particle it is stored in the neighbor list of
particle i, provided it satisfies the condition rii < hij. At the same time, if ri > hj for
neighbor j particle i is added to the list of neighbors of particle j, since these particles

will not be found when the neighbor search is done for particle j.

The adjustment of the smoothing lengths at each time step is adapted from the
procedure outlined by Steinmetz and Muiller (1993), which uses a predicted smoothed
density, p!, to find a new smoothing length:
(77T 1/2
hi = A (2.17)

where ( is set to 1.3 and T is the average mass of the neighboring particles. The
predicted smoothed density is calculated by using

S(dp \smooth
Pi = pfmooth + At- (2.18)

( P) smooth P
-d-[ = P--( v. v), 2V" )i'" (2.19)


(p2V v) = mp(V v)Wi (2.20)

All smoothed quantities, indicated by superscript, are estimated using the smoothing

length h* = hi/0.8.

This method maintains the number of neighbors, Ns, between 16 and 24 on the

average, except for the outer (diffuse) portions of the disks, where Ns can drop below 10.

In order to give reliable estimates the interparticle distance between neighbors should be

less than 0.5hi for the kernel used [equation (2.10)]. This requires that Ns 13 in two

dimensions, and Ns > 34 in three dimensions. If there are large density contrasts in the

outer regions, the adjustment technique may also cause N5 for a particle to jump to some

large value. Therefore, as a practical matter, hi is adjusted if N5 is greater than 60 or less

than 10. As this adjustment takes place while the neighbor lists are being compiled, errors

in the neighbor lists of nearby particles will be introduced, which must be corrected.

Hydrodynamic Equations

To evolve the fluid, represented by a finite number of particles, the acceleration and

velocity of the fluid at each particle is applied to the particles themselves. The particles

can then be thought of as being parcels of fluid, or alternatively, as representative particles

moving with the fluid. In either case, such a system's dynamic evolution is described by

the Lagrangian form of the hydrodynamic equations:


dt p VP + asc VPi,
dt p,

where Pi is the pressure, a"isc is the artificial viscosity term, and V'Di is the acceleration

due to the self gravity of the system, evaluated using the hierarchical tree method. The

subscript i indicates that the above set of equations is applied to each particle. The

acceleration on each particle is evaluated by estimating the terms on the right hand side

of the equation of motion using the SPH formalism. Care must be taken that each term

is in a form which is antisymmetric in i and j, so that the mutual forces of particles

will be equal in magnitude, but opposite in direction. In conjunction with the forces

being central, the antisymmetry of the force terms will insure that angular momentum is

conserved. Following Hernquist and Katz (1989),

1 ( __
--VP, + ai = m(2 + I7ij VWij, (2.22)
Pi J PiPj

where the artificial viscous term, taken from Monaghan (1992), is

-Iij 1 p vij rij < 0
I =V (2.23)

hvij rij (2.24)

In the above equations rij=ri rj, and cij = (ci + cj)/2 is the average of the sound

speed at i and j. The term 7 2 = 0.01h2 is to prevent singularities. The viscosity

parameters a and 03 control the strength of the artificial viscous term. An artificial

viscous term is introduced into SPH to model the bulk viscosity necessary to reproduce

shocks. Without the dissipation that this term effects, there would be an "interpenetration"

of particles, the random kinetic energy would increase, and the particles would no longer

describe the bulk motion of a fluid. However, it also introduces an effective shear

viscosity. To reduce the effective shear viscosity Benz (1990) has introduced a switch in

the form of a multiplicative factor on #ij, (Iij p.ijfij). The factor fij is the average

of fi and fj, where

(V.v) V v)
fi = V)1i (2.25)
(V v)|( + I(V x v)|. + 0.001ci/hi'

and the curl of the velocity field is given by the estimate

(V x v) = I mnvij x VWij. (2.26)

This factor effectively reduces the shear viscosity, but at the expense of not being able

to vary the effective shear. In order to use the parameter a to control the effective shear

viscosity, I do not apply the factor fij to the a term in the numerator of IIbj. In the

following subsection I present models of shocks using this hybrid artificial viscous term.

The set of hydrodynamical equations is closed with the inclusion of an equation of

state. Most simply a polytropic form can be used: Pi = KEi7, E being the surface

density. The index 7 is equivalent to the ratio of specific heats, and the constant K is

the square of the isothermal (7 = 1) sound speed. More generally, the ideal gas law

can be used in the form

Pi = (7 1)piui, (2.27)

where ui is the specific internal energy of particle i. The equation that governs the

evolution of the specific internal energy is the first law of thermodynamics, du =

-Pdp/p2 + Tds, where all nonadiabatic effects are included in the change in specific

entropy, ds. The form of the thermal energy equation that is used in T2DSPH is the

same as in Hernquist and Katz (1989):

dt m- Im) viWij + n flijvijVIWi + H + C. (2.28)

The first term of this equation is the adiabatic term, the second term is the non-adiabatic

viscous heating, and the final terms represent other non-adiabatic heating, H, and cooling

processes, C, not associated with viscosity. If the polytropic equation of state is used, then

the above equation is not needed. However, in this case, this equation is still integrated

for each particle with only the adiabatic term included. This is done so that the amount

of energy lost, due to the inclusion of an artificial viscosity with a adiabatic equation of

state, can be measured by measuring the change in the total thermal energy. (See further

discussion of this problem in section 2 of this chapter under adiabatic shock tests.)

Integration Method

Each particle has an intrinsic time step 6ti determined by its velocity, acceleration

and a Courant-like condition:

6ti = Cmin (i, () i .vl+c,+l.2(c,+h maxj,)) (2.29)

where the parameter C is 0.3. The equations (2.17) are integrated using a time centered

leap frog integrator, in which the velocities and positions are alternately stepped forward

in time a half step out of synchronization. For an individual particle this can be written as

n+ 1/2 n +Vn At
rn+/2 = r + v 2

vn+l = Vn + anAtl (2.30)

rn+1 n+1/2 n+1 AO
r 1 r + v. 9
i i 2 '


the superscripts being the time step index. If the equations for all the particles are

integrated simultaneously then the particle with the smallest time step will determine the

time step for the entire system. In accretion disks with large central masses, the range

in time steps can be large. To increase efficiency, multiple time steps are used, and

particles are assigned an individual time step that is a power of two subdivision of the

largest time step of the system:

At sy Atsys = max(6ti). (2.31)

Each particle, then, is assigned a time bin ni so that Ati < 5bti. Particles are always

allowed to move to a larger time bin (smaller time step), but may only move to a smaller

time bin if that time bin is time synchronized with the particle's current time bin. Forces

are then calculated for each particle i only once per Ati in order to step their velocities,

while positions for all particles are stepped every Atpos = Atsys/2fnm'ax+1. Advancing

the positions at every half step is necessary as forces need to be calculated at every half

step. Care is also taken that an estimated velocity vi is determined for all particles every

Atpos based on vi and the most recent estimate of ai as the viscous forces depend

on the local velocity field. For further details of implementing multiple time steps, see

Herquist and Katz (1989).

Gravitational Softening

Gravitational softening is a common, and necessary, convention employed in N-body

codes. Real gravitating systems often possess orders of magnitude more gravitating

particles than can be modeled directly-the computational requirements are simply


too great. However, such a system can be approximated by using a smaller number

of particles which interact via a softened gravity, which damps the effect of nearby

encounters. In modeling stellar systems, such as galaxies, the ratio of the number of

stars in a system being modeled over the number of gravitating particles employed in an

N-body code, is often as many orders of magnitudes as the number of particles. When

using an N-body code to model a gaseous system this difference between model and

reality becomes essentially infinite, as a finite number of particles is now representing

a continuous fluid.

The conventional method of softening is to use a form for the gravitational force

between two particles which possesses a small constant term in the denominator:

F = (r+ 2)3/2 (22)

where G = 1 is assumed throughout this work, and e is the softening parameter. However,

an alternative form of gravitational softening can be derived within the SPH formalism.

The gravitational potential of a continuous mass distribution can be written as

[ p(r')dr'
^ J r -r- r'|

Using the SPH estimate of the density [equation (2.9)] we have

\- fW(r'-rj, e)dr'
T(r) = -- m v -r = mjI,. (2.34)

The integral Ij can be evaluated by noting that

V2, = -47rW(r- r.),


from Poisson's equation. Since the operator (V-) = Irr2 in spherical coordinates we


r2VIj = 4r W(r rj, e)r dr, (2.36)

or using V., instead of Vr, and Vr = VuVruj, we have
2 u"
V2I = 4r3W(u)u2du, (2.37)
where V represents Vr. From equation (2.29) we can write

V = mjVI, (2.38)

and using equation (2.31) we can write

V{ = m- 3 W(u)u2du Vuj, (2.39)

where Vuj = uj/uj. Using the kernel
3/2(2/3 u2 + u3/2), if 0 < u < 1
(u) = 1/4(2 u)3, if 1 < u < 2 (2.40)
W0, otherwise
(equivalent to the kernel [equation (2.10)] with h -- 2e), and evaluating the integral

in the previous equation, a polynomial expression for the gravitational acceleration on

particle i due to particle j can be written as aij = -mjrijg(rij), where

S [ +1 u] 0 < u<1
g(r) + 3 + 5 ] 1 < u < 2 (2.41)
1/r3 u > 2.
For the sake of completeness the contribution to the potential of a particle pair can be

written p = -mjf(rij), where
2 3 2 4 151 5 O -3U-- UU+ 0 U
f(r) [u2- u3 + u4 u5] + e 1< u < 2 (2.42)
1/r u > 2.

The above expressions are given by Herquist and Katz (1989) without derivation,

citing Gingold and Monaghan (1978). This form of the softened gravity has the attractive

feature that it is equivalent to Newtonian gravity when a particle pair's separation is

greater than 2e.


During the hydrodynamic simulation an inner region, within a radius Ra of the central

particle, is treated semi-analytically in order to circumvent the modeling of the central

object (protostar), to avoid small time steps, and to define a radius at which particles

are accreted onto the central object. The radius Ra ~ 0.05RD, and particles initially

within this radius are kept and redistributed within Ra at each time step, so as to insure a

continuous boundary at Ra, and thereby provide pressure support for the gas outside Ra.

If, instead, no density profile within Ra is prescribed, but particles are simply removed

from the simulation upon entering this inner region, the rate of accretion will be dependent

on the size of Ra. In real disks, on the other hand, the accretion rate is determined by

the combined action of any global non-axisymmetric modes present, and by viscous and

magnetic forces acting locally throughout the disk. Since the physical condition of the

disk determines the accretion rate, I wish the same to be true in the simulations.

While the outer particles in the disk (r > Ra) are evolved with SPH, inner particles

are given an axisymmetric distribution so as to insure that the density and velocity are

continuous across the boundary at Ra. This is done by first extrapolating from the outer

disk the density E and the radial velocity Vr at Ra, which is accomplished by doing

a linear regression on ln(Ei), and Vri of particles in a small region around the inner

boundary. An estimate of V at Ra is made by assuming that the potential at Ra is due

to the central mass and an exponential disk, with a scale length 0.25Ra, in hydrostatic

equilibrium. Now it is left to specify E(r), V(r), and Vr(r) for the rest of the inner

region. One of two different density distributions were used: a massless polytropic disk

in a Keplerian potential, and an exponential disk. The initialization of these disks is

described in Chapter 3.

The particles in the inner disk are used as neighbors in the SPH estimations for the

particles located outside Ra to provide pressure support. To calculate the gravitational

forces, however, these particles are ignored; the inner disk and the central object are

treated as a single particle, with mass Me = Mc(initial)+(accreted mass). The central

region remains centered on the central particle which moves only under the gravitational

influence of the outer disk, which is treated as a nongaseous particle during the simulation.

Particles outside Ra are allowed to accrete into this region, at which point they are

removed from the simulation and their masses are added onto Me.

Treating the inner region in this semianalytical way introduces an inconsistency, in

that the gaseous mass within Ra (Mi), determined from the density profile of equation

(2.39) or (3-40), is not the same as the mass initially within Ra. Not only that, but Mi

will change as the boundary values of E and V change. However, the only purpose of

this inner disk is to provide a continuous inner boundary for the outer disk, rather than

being an attempt to physically model the disk within Ra.


Summary of Code Parameters

A list of the parameters which control the performance of T2DSPH is listed for

convenience in Table 2-1. The first parameter, N, determines how well a continuous

Table 2-1: Code Parameters

Parameter Description Value
N Number of particles. 5131
e Gravitational smoothing parameter 0.0272RD
0 Tolerance angle 0.7 radians
C Courant number 0.3
Ra Radius of inner region 0.0544RD
Smoothing length adjustment parameter 1.3
a, 3 Artificial viscosity parameters -, 1.5

fluid is approximated by the numerical method. While a large number of particles results

in a better approximation, the cost is greater computational time per time step. The next

two parameters control aspects of the gravitational force calculations, while the Courant

number determines the size of the time integration steps. The fifth parameter defines the

radius about the central massive particle at which SPH particles are removed from the

simulation, and their masses added to the central particle. The parameter C determines the

average number of neighbors, which should be greater than 13 for two dimensions, and

greater than 35 for three dimensions. The last two parameters are the artificial viscocity

parameters; 3 is not varied, while a is given one of three values to introduce varying

degrees of effective shear viscocity. In this sense a is treated as a physical parameter of

the gas. The nominal values used for the simulations presented in Chapter 4 and 5 are also


given. The units of length are expressed as fractions of the initial disk radius, RD. Tests

verifying that the chosen values result in accurate simulations are presented in Chapter 4.


To validate the accuracy of the hydrodynamics code, or at least to increase confidence

in its results, a variety of test models were run which can be compared to analytical

expectations. The first model presented is that of a stable, rotating Maclaurin disk. A

series of two-dimensional adiabatic shock fronts were also modeled to specifically test

the hybrid viscosity employed in T2DSPH (see previous section of this chapter).

Maclaurin Disk

A test to the code, using both SPH and self-gravity, is a simulation of a Maclaurin

disk, the two-dimensional counterpart to a Maclaurin spheroid, which represents a

polytropic stable solution to Poisson's equation and hydrostatic equilibrium when 7y=3.

The surface density profile is E(r) = o/ 1 (r2/R2,), and is initialized by radially

stretching a regular grid of particles (see Chapter 3).

Figure 2-2 shows the surface density and velocity profile after 26 dynamical times

(Tdynamic = RD/V(RD)), the dotted line representing the analytical solution. The

particle velocities lie slightly below the analytical solution because the code employs

softened gravity. Because the disk is in solid body rotation the effective shear viscosity

of the code has no effect on the disk during the simulation, that is, angular momentum

transfer is not induced.

0.6 -

0.4 -

0.2 -


0 0.2 0.4 0.6 0.8 1 1.2

Figure 2-2: Surface density and velocity profile of the Maclaurin disk after 26 dynamical times.

Two-Dimensional Shocks

Two-dimensional adiabatic shocks, without self-gravity, were modeled as a means of

specifically testing the implementation of SPH in the program T2DSPH, and as a test of

the hybrid viscosity that is employed (see previous section). A two-dimensional sheet,

with an initial density equal to one, was initialized by placing particles on a regular

grid. Particles with a positive x coordinate are given an initial negative x component

velocity of Mach 1, and those with a negative x coordinate are given a Mach 1 velocity

in the opposite direction. This results in the formation of an accretion shock, having

two boundaries parallel to the y axis traveling away from each other. The sheet has a

finite extent, resulting in rarefaction waves traveling inward from the edges of the sheet.

However, for a time, the central region of the sheet will be uncontaminated by these

waves, allowing the simulation there to be compared with an analytical solution.

For all the tests presented here the polytropic equation of state was used, with 7 = 2.

From the jump conditions of an adiabatic shock, and using the polytropic equation of

state, a solution can be found for the shock velocity and post-shock density. However,

this will not be the desired solution. By employing an artificial viscosity, at the same

time that we use a polytropic equation of state, we have introduced an effective cooling,

because the kinetic energy dissipated by viscosity is not added to the system as heat.

This thermal energy would be added to the material if the specific thermal energies of

the particles, ui, were evolved:

dt 2 SmliJiJvi (2.43)
S/ non-adiabatic

In other words, in spite of our polytropic equation of state, we do not have a purely

adiabatic shock. This can be taken into account in our analytical solution by the addition

of a post-shock cooling term, Q (heat lost per unit mass), in the energy jump condition.

Therefore, using the polytropic equation of state, with K=1, the jump conditions are:

p1V1 = P2V2,

piv1 + p = p2v2 + p, (2.44)

1(7 -1 2 (7 1 2

which results in the final solution:

-b + Vb + 4 (7 -1)(v/2 Q) + pp 1

v2) =

P2 Vo v2
P2 Vo + V2


b = 2-vo + Q.
2 L?'o

The solution above is given in the post-shock reference frame, with vo and v2 representing

the pre-shock velocity and the shock front

velocity respectively.

10 0- .. .- .. Ix I

0 M
a N
0 0 0 -

00o ^ g "1n 00"- *'yrii11

~ *00~~00

a 0


Figure 2-3: Surface density of three two-dimensional accretion
shocks. The bottom accretion shock's surface density is offset from
the top by -2, and the middle by -1, from the top accretion shock.






In Figure 2-3 three simulations are shown. The bottom two, with a = 0.25, 1, in

ascending order, and /3 = 0, do not use the viscosity switch. In the third simulation,

shown at the top of the figure, the hybrid viscosity is employed, with a = 0.25, '3 = 1.5.

Although this solution is not as good as the first, it has much less effective shear viscosity.

Also, notice that there is no interpenetration of particles, as there is in the first case.

Analysis of Disks

Evaluation of Effective Shear Viscosity

In the SPH formalism shear viscous forces are not introduced directly, but instead

are a result of the artificial viscous term, av"'i, being evaluated over a finite region.

Hernquist's form of the artificial viscous term illustrates this more explicitly; being

dependant on (V v) it formally introduces no shear force but is purely a bulk viscous

force. However, what is used numerically is an estimate of (V v), which is evaluated

by looking at particles within a finite region (neighborhood) around the point of interest.

In addition, the total force on a particle is a sum of symmetric forces between itself and

neighboring particles at other radii, and the tangential components of these forces transfer

angular momentum. These errors in the estimation formally vanish as the smoothing

length approaches zero.

This implicit introduction of the shear viscous forces in SPH causes an unfortunate

situation for those attempting to use this method to model astrophysical disks, where

large velocity gradients can exist and shear viscous forces play an important role. In

SPH both bulk and shear viscous forces are effectively coupled together, rather than


being independently parameterized. In an attempt to decouple the bulk viscosity needed

to describe shocks and the effective shear introduced by the artificial viscous term, an

alternate form for the artificial viscous term was introduced in section 2. Here I wish

to address how the effective shear viscosity scales with the viscous parameter a or any

other quantities. The functional form of the artificial viscous term suggests that, like

the av viscosity of accretion disk theory (Shakura & Sunyaev 1973), the effective shear

coefficient scales as a,ch, c being the sound speed and h representing a characteristic

length. Artymowitz (1993) has suggested that for SPH av, 2 O.1a, with the characteristic

length being considered as the smoothing length. While the similar functional form of

the artificial viscous term and the viscosity of accretion disk theory is suggestive, it is no

guarantee that the implicit shear will behave similarly, being a by-product of estimation.

It is therefore desirable to determine not only the relative magnitudes of a and av, but

also whether the effective shear scales in the way that the artificial term suggests.

To investigate the effective shear associated with the artificial viscosity three models

were run with a = 1, 0.5, 0.25. An exponential disk was initialized, with the velocity field

now consistent with assuming that the disk is in a Keplerian potential with no pressure

support. These disks were evolved with SPH consistent with these assumptions, and they

remain axisymmetric in the absence of self-gravity. Since no pressure support is required

it was found to be advantageous to exclude the inner disk, and to simply remove particles

from the simulation when they approached within Ra of the gravitating particle located at

the origin. After a few dynamical times a radial velocity flow is established throughout

most of the disk as a result of the effective shear viscosity.

From standard accretion disk theory the radial velocity profile of an axisymmetric

thin disk, with no radial pressure support, is given by

vr = [r(r fr2) (Evr3'), (2.47)

where f is the angular velocity, and the primes denoting radial derivatives (Pringle 1981).

If the shear coefficient v is not a constant, then
v = (Er )' IVrEr(r 2)'dr. (2.48)
For a Keplerian potential f22 = GM/r3, which upon substitution results in
v = ( Er 1/2) I / vrdr. (2.49)
If in SPH the shear coefficient, v, scales as a,,ch, then a,, = v/ch will be constant

over the disk, and can be estimated by evaluating the above expression for v. The SPH

particles of the evolved disks are used to find the surface density and radial velocity

of annuli centered on the origin, and the integral in equation (2.48) is estimated by the

sum over the annuli

v.ri Ei r1 Ar, (2.50)

where Ar is the width of the annulli and is equal to .05RD. Now at each radius ri the

shear viscous parameter av can be estimated by using

a,(ri) = 2 (2.51)
1.3 (mp
remembering that the sound speed c = v/2/
h = 1.3 -mp/E, where mp is the mass of each SPH particle and assuming that the

smoothed density E* = E. Using this method on the three models described above

yielded a, approximately constant with radius, with values of 0.1, 0.05,0 .02 .01

for each model respectively. Much of the uncertainty arises from a time variability that

is probably due to radial oscillations. The fact that av is found to be approximately

constant over the disk confirms that the effective shear viscosity scales the same way as

the artificial (bulk) viscosity.

Evaluation of Modes and Their Frequencies

To follow the evolution of the nonaxisymmetric modes a method similar to that of

Papaloizou and Savonije (1991, hereafter PS) is adapted; a polar grid, concentric with

the center of mass of the system, is imposed on the disk with a radial extent from 0.1

to 1. The density is evaluated at 64 points, equally spaced in azimuth, for each of 25

equally spaced radii, and then normalized by dividing the density at each point by the

average density at each radius. The modes are identified by doing a Fourier analysis of

the normalized density in azimuth at each radius, and plotting the power spectrum as a

function of mode number and radius (Figure 2-4). Unlike PS, the density is not taken to

be the mass/area of individual cells defined by the grid, but instead the SPH formalism

is used to estimate the density at each grid point:

k = rmiki, (2.52)

where the sum is over all particles having the grid point k within their smoothing length

hi. This analysis can be done at each time that the position and smoothing lengths of the

particles are output by the SPH code, which is at every 0.5 dynamical times (see Chapter


3 on choice of units). Another useful way of presenting the results of this analysis is with

a contour plot of the dynamic spectrum: the power in a particular mode as a function

of radius and time (Figure 2-5). The evolution of a particular disk can be more simply

characterized by considering the maximum power, with respect to radius, of each mode

as a function of time (Figure 2-6).

To find the frequencies that are present, the SPH code outputs the estimated density

at 100 points along both the x and y axes of a Cartesian grid, with its origin at the center

of mass of the system, at every third system time step. After the simulation a Fourier

analysis is done in time, at each radius, for a specified time interval. Contour plots of

the power, as a function of radius and frequency, can be displayed separately for each

axis, or averaged together (Figure 2-7). Provided that there are a limited number of

modes, the contour plots of the modes and of the frequencies can be used together to

match frequencies with particular modes.





0 2 4 6 8 10
INT=0.05, BOT=0.01 TIME=21.0055

Figure 2-4: Normalized power spectrum of modes of
model A2 (see Chapter 4) at Time = 21 dynamical times.

0.8 V-*V

0.6 -2i

0.4 -


10 15 20 25 30
INT=0.05, BOT=0.05

Figure 2-5: Dynamic power spectrum of the normalized density for mode = 2 of model A2.

0.8 -- '-
._.. m=4

S0.6 I

\ IF
' ".. ,,

S' 1 \ '

/ 11 / /A
o ,I; / ,I I ,

0 10 20 30 40

Figure 2-6: Maximum of the power spectrum, in
modes = 1 4, of the normalized density in model A2.

15 I I I


5- ^ ^

0 0.2 0.4 0.6 0.8
BOT=13.2012=99.9% cont. Time = 20.0 to 30.0

Figure 2-7: The average Lombe periodogram of normalized
density temporal fluctuations along Cartesian axes for model A2.


In this chapter the initialization of the disks numerically evolved with SPH is detailed.

Initializing a disk for an SPH program involves finding a distribution of points in phase

space which suitably describes the particular density and velocity field of the disk. Except

for perturbations, all initial disks used in this work are axisymmetric and in radial

hydrostatic equilibrium. Global parameters describing a disk are the disk mass MD

and the radius of the disk RD. With the exception of the Maclaurin disk, and in the spirit

of modeling protostellar accretion disks, there will also be a central particle with mass

Me, giving a total mass for the protostar/disk system of MT = Mc + MD. Treating the

central object as a point mass is equivalent to assuming it is several magnitudes smaller

in size than the disk radius. This assumption is valid for protostellar disks that have radii

approximately 10 1000AU. Parameters of the gas itself are the ratio of specific heats,

7, and the isothermal sound speed K. Unless stated otherwise, dimensionless units of

mass, length, and time will be used for which MT = RD = G = 1. In this case the unit

of time then becomes the dynamical time: TD = RD/V(RD) = (R'/GMT)1/2. The

models may then be scaled to the desired dimensions. For instance, if MT = 1M and

RD = 100AU then the dynamical time is 159 years.

Two-dimensional Disks

Maclaurin Disk

A two-dimensional disk has already been presented in Chapter 2 as one of the

tests to the SPH code: the Maclaurin disk. This is a two-dimensional version of the

Maclaurin spheroid, which represents a polytropic stable solution to Poisson's equation

and hydrostatic equilibrium when -y=3. The surface density profile of the Maclaurin

disk is

E(r) = EV/1 (r2/R2), (3.1)

and is initialized by radially stretching a regular grid of particles. The central density

Eo is given by 3MD/2rRj.

This method of constructing an axisymmetric disk I call the stretched grid method,

and is equivalent to transforming the radial coordinate, R, of a particle in a uniform disk

to the new radial coordinate r via a coordinate transform function: r = A(R)R. To apply

this method the transform function, or stretching factor, A(R), must be found which will

transform a uniform disk into a Maclaurin disk. The desired function can be found by

considering an invariant of the transform, the mass within a given radius. That is, the

mass within a given radius of the uniform disk, MDR2/R2D, is equivalent to the mass

within r = AR of the Maclaurin disk, which is

MD 1 (1 R)3/2. (3.2)
\ "D/


Setting the above expression equal to MDR2/R2 replacing r with AR, and solving for

A gives the desired transform function:

A(R) = R 2/3 1/2(3.3)

Giving each particle in the uniform disk with radial coordinate Ri < RD a new radial

coordinate ri = A(Ri)Ri, and discarding any remaining particles, transforms a uniform

grid into a Maclaurin disk (see Fig 3-1).


0.5 h

>- 01-

-1 -0.5 0 0.5 1
Figure 3-1: Initial distribution of particles for the Maclaurin Disk.

This method, though quite general, was found to be of limited usefulness when it

was used to construct disks that possess stronger central mass concentrations, such as the

polytropic Keplerian and exponential disks described below. The undesirable feature of

the disks built with this method is the presence of artifacts of the grid from which the disk

was stretched: the direction to the nearest neighbor is not isotropic on the average, the

particles instead appearing to be "planted" in rows, and the edge of the disk is corrugated,

or scalloped. These artifacts are not severe in the Maclaurin disk shown in Fig 3-1, but

they become much more pronounced in centrally condensed disks.

Once the particles have been spatially distributed their velocities are determined.

The initial tangential velocity of each particle Vi in the Maclaurin disk is kVo(ri),

Vo(ri) = Cri/RD being the initial tangential velocity required if the disk were supported

purely by rotation, with C2 = 3rGMD/4RD. The rotation parameter k = V/Vo = 0.4 is

assumed to be a constant with respect to radius, and gauges the importance of rotational

support. Maclaurin disks with k < 0.5 will be stable against secular instabilities (Binney

and Tremaine, 1987). The tangential velocity Vo was found from the gravitational force

on each particle, evaluated using a direct sum over all other particles, rather than using

the analytical formula for Vo, because of the softened form of gravity employed. From

the requirement of radial hydrostatic equilibrium,

1 9P V2 V2
= (3.4)
E or r r

which follows from the radial component of the equation of motion in cylindrical

coordinates, and the assumptions of axisymmetry and no radial motion, the constant K in

the equation of state is found. It is determined to be a function of the disk parameters:

Gr R3
K= (1-k2) RD (3.5)

The above expression is found from assuming unsoftened Newtonian gravity, whereas

the code that will evolve the initial conditions uses softened gravity. Because the value

of the softening parameter is small, the error introduced by this inconsistency is also

small. Using k = 0.4 and MD = RD = 1 gives us K=2.8939.

Exponential Disks

Exponential disks were chosen to simulate accreting protostellar disks in large part

due to the previous theoretical and numerical work that has been done with such disks

(Papaloizou and Savoneji, 1991). This is the practical reason. That such a choice is

reasonable physically is somewhat born out by the numerical result that the exponential

profile is preserved over most of the disk for the majority of the duration of the numerical

evolution, at least in the case where the central star to disk mass ratio was equal to three.

The models with equivalent star and disk masses differ significantly in that a m= 1 mode

completely dominates the mass distribution of the disk and erases all axisymmetry. In any

event, angular momentum transport and mass accretion reinforce the centrally condensed

mass distribution.

To represent the mass distribution of these disks, particles representing the gas, all

of equal mass, are placed in concentric rings around the central particle at the origin.

The surface density profile,

E(r) = Eoexp(-r/r,), (3.6)

is achieved by appropriately spacing concentric rings, from the inside out, so that the

correct interparticle separation, (m/E(r))1/2, is nearly imposed. The scale length r, is


set to 0.25 in all models, while the ratio of the central to disk mass Mc/MD is either 3

or 1. The central density is then given by the expression

o = 1 -e-RDtl+ ~. (3.7)
27rr; \ r, /

Before the disk is built the number of total desired particles, Ntot, is specified. The

gaseous particles are then assigned a particle mass of mp = MD/Ntot. The first ring

of gaseous particles consists of six particles placed at equal angular intervals around

the central particle at a distance 0.5 x r(6mp), one half the radius within which a mass

of 6mp resides in. Here, and in what follows, the radius r(mr) is found by using the

Newton-Raphson method, with the expression for the mass that lies interior to a given

radius for the exponential disk:

m(r) = Eo27rr [1 e-r/r'(r/r, + 1)]. (3.8)

For subsequent rings, the ring radius rn is determined by finding the number of particles

in the ring, Nn, that minimizes the quantity

(rn rn-1)-2 (3.9)

where r2 = (r2 + r_1 ) /2, with each primed radii being the outer radius of the annulus

about the nth ring. This is an attempt to equate the distance between ring n and n + 1

with the distance between particles in ring n. To calculate the quantity above, the radius

of an nth ring with Nn particles is found from r_1, and the radius rn within which
there are Nr = E Nj particles. This iterative procedure is continued until Nr > Ntot. If
Nr > Ntot then Ntot is set equal to Nr, the particle mass mrn is reinitialized, and the disk


rebuilt. For an exponential disk this procedure will converge on a value of Ntot which

will result in a disk being built with Nr = Ntot at the outer most ring. In the case where

5000 particles are initially requested, the above algorithm results in an exponential disk

with 5130 particles. Though the above procedure for building an axisymmetric disk with

particles on concentric rings is in principle quite general, requiring only that the desired
function m(r) = 27 f E(u)udu be specified, it is not known whether the iteration will
converge for other density profiles.

An initial circular velocity V is assigned for each ring of particles according to the

equation of radial hydrostatic equilibrium, which leads to

V2 = Kr,-2d + V = -K r '-- + V,, (3.10)
dr r,

where V12/r for each ring is calculated by finding the radial gravitational acceleration for

each particle, through a direct summation, and averaging the radial components of the

acceleration of the particles in each ring. In doing this calculation softened gravity is used.

The value of the constant K, the square of the isothermal sound speed, is determined by

assuming radial hydrostatic equilibrium throughout the disk and specifying a minimum

value for Toomre's stability parameter Q = Kc/TrGE, where K = (- + 2 )1/2 i

the epicyclic frequency. Toomre's parameter is a local stability parameter, derived from

the linearized hydrodynamic equations of a local region of a rotating gas sheet, which

describes the stability of a two-dimensional disk (Toomre 1964). When Q > 1 the disk is

stable against axisymmetric perturbations, but the disk generally remains unstable to non-

axisymmetric modes for values of 1 < Q < 3. The epicyclic frequency, K, is numerically


evaluated for each ring n with the expression

1 (V2(rn+) V'2(r-)) V (r)
n= 2 (3.11)
rk ((n+1 rn-1) r.
Using the expression c2 = 7E'~1, Toomre's stability parameter can be written as

= r (3.12)

An expression for the constant K can be found from specifying Qmin using the previous

K- (Qmin7rG)2
-- (3.13)
/ r3-7y
7mink Kkk
For all the two-dimensional exponential disks, Qmin is set to 1.15. In Figure 3-2 is shown

the initial velocity profile and the value of the Toomre Q parameter with respect to radius.

0 0.2 0.4 0.6 0.8 1

Figure 3-2: The initial Toomre Q parameter (solid line) and tangential
velocity (dotted line) with respect to radius for a Mc/MD = 3 disk.

A density perturbation is then seeded into the disk by adding independent Gaussian

noise to the x and y coordinates of each particle, the displacement having a standard

deviation of 0.02 of the interparticle distance. This perturbation is added after the

tangential velocity V is determined for all of the rings. The subroutine which generates

the Gaussian noise is taken from Numerical Recipes (Press et al., 1992). This introduces

a density fluctuation with an rms amplitude 0.03 times the average density. The purpose

of introducing the noise is twofold. The first is to excite all possible modes, so as to

observe which modes grow fastest and those which tend to dominate the disk. This was

the primary reason for the introduction of noise. The second, more serendipitous, is to

avoid a flaw in the initializing method. If the SPH particles were distributed on concentric

rings, without noise, then the regularity of their distribution will suppress some possible

responses while enhancing others. For instance, it was found from evolving a disk with no

noise, but with an m=2 perturbation (see below), that the initial growing mode developing

early in the disk's evolution, while particles still lay on concentric rings, led to the rings

being distorted, and themselves becoming structures that excited a response. The folded

and compressed portions of the rings formed shock "wakes" to the arms.

The FORTRAN code which generates the initial positions and velocities of the

particles that describe the exponential disk is listed in Appendix C. The initial distribution

of points for the exponential disk is shown in Figure 3-3, while the initial radial profile

of the density is shown in Figure 3-4.

0.5 -

> 0-



-1 -0.5 0 0.5 1

Figure 3-3: Initial distribution of particles in an exponential disk.

In one model an m=2 density perturbation is introduced of the form E'(r, p) =

Eo(r)p(r)cos (my), where p(r) = Ap[7r(r rs)/(RD r,)], with Ap being 0.01, and Eo

representing the unperturbed density given by equation (3.5). This perturbation is induced

by displacing the particles along the circumference of their respective ring. To find the

displacements which will effect the desired density perturbation in a particular ring, I

consider the perturbed density in the form E'(k, x) = Acos (kx), where A = Eo(r)p(r),

k = m/r, and x = pr, an arclength. This form of the perturbed density can be considered

-1 "


0 0.2 0.4 0.6 0.8 1 1.2

Figure 3-4: Initial density profile of an exponential disk with MD = 0.25 and a scale length of
r, = 0.25. Each point corresponds to the estimated density at each particle position,
evaluated with the SPH method, while the solid line is the intended density profile.

as being the perturbed density of a one-dimensional sound wave at time equal zero,

which is more generally written as E'(x, t) = Acos(kx wt). Linearization of the one-

dimensional equation of continuity and motion gives the wave equation

= O' (3.14)

where Co is the sound speed equal to w/k. Using this equation and the given form of

the perturbed density, I find the perturbed velocity v' = cop(r)cos(kx wt). From this

perturbed velocity a displacement can be found:

Ax = v'dt = P(sin(kx wt). (3.15)

The desired displacements are those at t = 0, at which time an equation for the perturbed

particle position x is formed:

x = o -p sin(kx), (3.16)

Xo being the unperturbed position, and Ax = x Xo. This equation is solved for each

particle using the Newton-Raphson method.

To find the associated perturbed velocities, the velocities of the one-dimensional

sound wave employed above are not used. While the above derivation is sufficient for

finding the displacements, it is an oversimplification with regards to the gas dynamics

of the disk. Instead, the equation of radial hydrostatic equilibrium is used to find

the tangential velocity V(r, ;) by inserting the perturbed density. This results in the


V2() = V02 -7 -K S 1 p(r)cos(my)(1 + ( )cot[ r(r-r }.


Another model that was evolved includes an encountering particle, with a mass of 0.5,

that approaches the central particle with a minimum distance of roughly 1.0. Otherwise

the disk and central particle are initially identical to model D2. The initial separation

between the encountering particle and central particle is about 8.772. The energy of the

encountering particle is specified to be

1 GMme
E (3.18)
2 a

where M is the mass of the central object and disk (1.0), me is the mass of the

encountering particle (0.5), and a is the separation. With this energy the relative velocity


of the particle is v = V/3GM/a. A common parameter to describe an encounter is the

impact parameter b. Instead of specifying this parameter I have specified the minimum

separation amin = 1.0. The impact parameter can then be determined as

b= ,maamin (3.19)

7( GMme, 2
Vmax = 2E + (3.20)
amin me
The energy of the encountering particle can also be parameterized in terms of the velocity

of the particle if removed to infinity, Vo = W-FE. The specific angular momentum of
mespecific angularmomentumof
the encountering particle is then bVo. Using this, the components of the velocity can then

be determined: VT = bVo/a, and v2 = v2 v2, being the velocity tangential and parallel

to the separation vector respectively.

Inner Disks

As described in Chapter 2, a small inner region about the central particle is described

analytically, rather than being evolved with SPH. Two different inner disks are used.

If the disk mass within Ra is much less than Mc, it may be considered as a massless

Keplerian disk. Assuming axisymmetry in two dimensions, radial hydrostatic equilibrium,

and using a polytropic equation of state with 7 = 2, the density distribution within Ra

may be written

E(r) GM(1 k 2)) + E(Ra). (3.21)

It has also been assumed that k = V/Vo = V(Ra)/Vo(Ra) is a constant throughout the

inner disk, and Vo is the circular rotation without pressure support. The inner particles are


repositioned on concentric rings to describe this density distribution, with their masses

being reassigned so that mp is equal to .V,n divided by the number of inner particles,

where Ain is the mass interior to Ra as determined from E(r) in the equation above.

Using Vo = GMc/r for the inner disk, the tangential velocities Vi = kVo(ri) are

determined. The radial velocity is specified by assuming that the mass accretion rate

rh = 2rE.(r)Vr(r) is constant throughout the inner disk. Since the values of E, V, and

Vr at the boundary are determined at each time step, the inner disk must be reinitialized

at each time step as well.

The inner disk described above is specifically for 7 = 2; for other values of -y an

alternate model of the inner disk must be used. The second inner disk that is used is

one with an exponential density profile,

E(r) = Eoe-r/r.. (3.22)

The scale length is found from the linear regression used to estimate E(Ra), and the

central density is deduced from E(Ra) and the scale length, r,. Again, as a fixed number

of particles is used to model this inner disk, their individual masses are being changed

during the run, and as in the first inner disk described, they are placed on concentric

rings. The velocity profile for the inner disk is found from assuming radial hydrostatic

equilibrium and using Vo = v'GMc/r + Ar, the parameter A being estimated at Ra. As

in the previous disk, the radial velocity profile is determined by assuming a constant

accretion rate rh throughout the inner disk.

It has been stated above for both of the inner disks that the particles are distributed in

concentric rings. I now wish to describe more carefully the algorithm used to accomplish


this. In order to assure a continuous particle distribution at the boundary r = Ra,

the inner disks are built from the outside in. The first ring radius r, is the radius

at which the interparticle distance (mp/E(ri))1/2 is equal to 2(Ra rl), and is found

using the Newton-Raphson method and the appropriate density profile (equation 3.17

or 3.18). Particles placed in this and subsequent rings n, will describe the desired mass

distribution subject to the condition that the mass in each ring, Amn, is equal to the mass

in an annulus, centered on the ring radius rn, with a width Arn equal to the interparticle

distance. Hence, for the first ring, the mass contained within Arl = (mp/E(ri))1/2 is

equal to Ami = m(Ra) m(Ra An1). [In general m(r) will signify the mass within

radius r.] The number of particles to be placed in the ring is determined by the mass in

the annulus, Am,, and the particle mass, mp. However, since the number of particles

NV placed in a ring must be an integer, and Amr./mp is in general not an integer, a

small adjustment in r, and Arn must be made. For the first and subsequent rings, Nn is

found by rounding Amn/mp to the nearest integer value. Then the mass in the ring is

redetermined: Amn = Nnmp. For the first ring a new Art is redetermined by finding

the radius rM for which m(rM) = m(Ra) Ami, using the Newton-Raphson method

and the appropriate expression for m(rM), and then setting Art = Ra rM. Finally,

r, = Ra Ari/2. For the remaining rings the same procedure is used to redetermine

Arn from the new Amn, excepting that, instead of Ra, rL = rn-1 Arm_1/2 is used.


For clarity's sake, I restate the algorithm for finding the radii of the concentric rings:

1. Using the Newton-Raphson method, find r, which satisfies the condition

2(rL rn) = (mp/E(rn))1/2, where rL = Ra for n = 1, but rL = r,n-1 Arn_1/2


2. Arn = 2(rL r,)

3. Am, = m(rL) m(rL Ar1)

4. Find the number of particles in ring n, Nn, by rounding Amn/mp to the nearest

integer value.

5. Redetermine the mass in ring n: Amn = Nnamp.

6. Using the Newton-Raphson method, determine the radius rM for which m(rM) =

m(rL) Aman.

7. Set Arn = rL rM.

8. Finally, rn = rL Arn/2.

9. Proceed to next ring.

The building of the disk is discontinued when rM < Arn/2.

In order to use the above algorithm, the expressions for the mass with a given radius

for the two inner disks must be known. For the inner polytropic Keplerian disk:

m(r) = (1 k2) r 2 + 2(Ra), (3.23)

where k here is V/Vo and shouldn't be confused with the index. The expression for

m(r) for the exponential disk is given above in equation (3.7). For the polytropic disk

the above expression can be inverted to find r(m), the radius within which there is mass

mn. This will allow us to circumvent the use of the Newton-Raphson method in step 6)

above, and instead use the expression:

rM = 2 (3.24)
B + /1B2 + 4(7r(Ra) B/2Ra)m,

where B = rGMc(1 k2)/K, and mr = m(RL) Amk.



To gain insight into the general behaviour and evolution of protostellar disks, a suite

of models of accretion disks were simulated. Three physical parameters were varied:

the central object to disk mass ratio, Mc/MD, the ratio of specific heats, 7, and the

artificial viscosity parameter a. As I am particularly interested in the early stages of the

evolution of protostellar disks when the disk is self-gravitating, values of Mc/MD of 3

and 1 were chosen.

The second physical parameter that is adjusted, 7, affects the hydrodynamical char-

acter of the gas. In particular, as it appears as an exponent in the polytropic equation of

state, 7 determines the compressibility of the gas. While in three dimensions a value of

-7 < 4/3 is unstable to gravitational collapse, in two dimensions 7 < 3/2 is unstable to

collapse. This is one example of how gravity is more effective in two dimensions than

in three. Values of 7 = 2 and 5/3 were chosen for modeling; -7 = 2 corresponds to a

gas with two degrees of freedom, while -y = 5/3 to a gas with three degrees of free-

dom. The value -3 = 2 was chosen because previous theoretical and numerical work on

two-dimensional disks have used this value (Papaloizou and Savonije, 1991). However,

though convenient for solving the hydrodynamical equations, 7 = 2 results in a stiff equa-

tion of state. Therefore -y = 5/3 was also used in the simulations. This value of 7 allows

the gas to be more compressible, though it is still stable against gravitational collapse.


Of particular interest in this study is the effect of viscosity on the evolution of self-

gravitating accretion disks, and the parameter used to vary the effective viscosity in SPH

is the artificial viscosity parameter a. Values of a = 1, 0.5, 0.25 were chosen, which

were shown in Chapter 2 to correspond roughly to effective shear viscosity parameter

values of a, = .1, .05, and .02. On the basis of time scale arguments applied to

real astrophysical disks, a, is expected to be less than 0.1, where values between 0.04

and 0.002 are most commonly argued for. As mentioned in Chapter 1, these values are

much greater than those from molecular viscosity alone. Due to the as yet unresolved

question of what mechanism effects this anomalous viscosity, the value of av represents

the greatest unknown in the physics of accretion disks.

Varying the three parameters -1c/0MD, 7, and a as described above gives twelve

possible combinations, all of which are modeled. Table 1 summarizes the models

discussed in this chapter, as well as the time at which each simulation was terminated,

in dynamical times (see Chapter 3 for definition of units). All of the Mc/MD = 1 disks

terminate due to computational errors associated with adjusting the smoothing length

during the neighbor searching phase of a time step, and the finite size of the arrays in

which the lists of neighbors are stored. That is, the program has difficulty in adjusting the

smoothing length so that particles have more than ten, but more than sixty, neighbors.

This problem arises in disks in which an extreme density gradient borders a region

where the surface density, and hence also the number density of particles, is very low.

Such conditions arise in the disks with a mass ratio of one. In contrast, most of the

Mc/MD = 3 models ran to the specified time without difficulty. Their simulations


were terminated after they had run for a time comparable with the .1[c/MD = 1 disks.

Model A3 is the sole model with a mass ratio of three that ended due to the type of

difficulty mentioned above. In any event, the lengths of all the simulations were sufficient

to observe important properties of accretion disks, which are noted below and in the

following chapter; continuing the simulations much further, after a significant number of

particles already have been removed due to accretion, would be of questionable value.

Table 4-1: Models


Inner Disk





While the model parameters discussed above describe the physical state of the initial

conditions, another set of parameters which may affect the simulations is that of the

code parameters, which determine details of the numerical simulation method. These

parameters are listed in Table 2-1, with their nominal values. Ideally, if the values





chosen for the code parameters are sufficient to yield accurate simulations, the resulting

simulations will be insensitive to variations in the code parameters. To test whether

the simulations are accurate, several test models have been numerically evolved, each

with the value of a single code parameter changed from its nominal value, so as to

yield a more numerically accurate simulation. If the nominal value is sufficient, then

the test simulation will show little or no change from the model run with nominal code

parameters. A number of such tests are done with model D2 for the parameters 0, C, and

N. The comparison between the tests and model D2 are given at the end of this chapter.

The gravitational smoothing parameter, e, is not varied for testing, because its range

is restricted by requiring that it be smaller than the accretion radius, Ra, and that it must

be larger than the local inter-particle separation. This last condition must be satisfied

if N particles are to behave as a system with a much greater number of particles, or

even a continuous medium (N -+ oo), as is the case here. Because e is of single value,

independent of the local number density, the latter condition cannot be satisfied in the

outer, diffuse, portions of the disks. The local inter-particle separation in two dimensions

is n-1/2, with n representing the local number density. An average inter-particle distance

is defined as n-1/2 = V/7rRD/N, which is equal to 0.025 for RD = 1 and N = 5131.

The condition e > fi-1/2 is satisfied with e = 0.027, but e < n-1/2 for r > 0.587 when

N = 5131 and rs = 0.25. When N is increased to 10093 particles the condition fails

for r > 0.757.

The other restriction, that 2e < Ra, is imposed so that all particles experience

a Keplerian, unsoftened, gravitational potential from the central object. Therefore


increasing Z is possible if R, is also increased. However, this is not desired, as the

inner region around the central particle, which is not hydrodynamically modeled, should

be a small percentage of the disk.

Evolution of Disks

General Features

Since the evolution of the models is terminated at different times, it is more useful

to compare the evolution of the models at a time that all the models attain. In Table

4-2 are tabulated the total accreted mass, the central object to disk mass ratio, and the

angular momentum (AM) that has been transported beyond the initial outer radius of the

disk, after twenty eight dynamical times. The units of mass and angular momentum are

the dimensionless units described in the previous chapter.

To qualitatively show both the effect of the viscosity parameter a and the ratio of

specific heats 7, the models are shown after ten dynamical times in Figures 4-1 and 4-3,

and after thirty dynamical times in Figures 4-2 and 4-4. For both the Mc/MD = 1 and

3 models, increasing the viscosity effectively damps the amplitude of the modes. This

is even more noticeable at early times, where viscosity also suppresses the growth of

the modes. Of particular of interest is that the viscosity has suppressed the m= 1 mode

in the models having Mc/MD = 1. The mode is weakly visible in model B3 in Figure

4-4, but not visible in model D3, and a dominating m=1 mode remains absent in this

last model until the simulation is terminated at T=32. This later model represents the

one exception in the modal evolution of the Mc/MD = 1 disks, which otherwise become

dominated by a slowly growing m=l mode. This mode is of a tidal nature, resulting


from the gravitational interaction between the massive central object and the disk, and

has been predicted by others (Savonije et al., 1992).

Table 4-2: Models at T=28

mass (10-2)

dynamical times


AM (10-2)

The effects of varying the ratio of specific heats are more subtle. At early times,

when the modes are growing, the models with smaller -7 show slightly more density

contrast between arm and inter-arm regions. At the later times, such as those shown in

the Figures 4-3 and 4-4, this trend is not apparent in the a = 0.25 and 0.5 disks, and is

reversed in the a = 1 disks. Another effect that is more obvious in the a = 0.25 and

0.5 disks is that those with -y = 5/3 tend to show greater density variations along the

arms of the non-axisymmetric modes. At times this density variability will develop into

clumping in the outer portions of the arms, such as that seen in model D 1.



A more detailed description of the evolution of the modes, and of the mass accretion

observed in these disks are given in the next two sections. In the Mc/MD = 1 disks

another distinctive feature is seen, and that is the formation of small satellites in the outer

portions of disks. This is discussed in the following chapter.

Modal Evolution

The evolution of the simulated disks is typified by the presence of many transient

modes. In none of the Mc/MD = 3 disks does a single mode dominate the disk for more

than a few dynamical times. In all disks the growth of the non-axisymmetric modes

follows the general pattern that modes with higher mode number show faster growth

rates. However, as the amplitude in all modes increase, the modes of higher mode

number will peak and then decrease, while the modes of lower mode number continue to

grow. Hence, as the amplitude of the density fluctuations increases during the early part

of the evolution, different modes achieve temporary dominance. The maximum amplitude

is attained usually by an m=2 or m=3 mode. When the maximum power in a mode has

diminished it will usually not remain suppressed, but will later increase again to become

the mode with maximum power. In other words, the modes will initially grow at different

rates, but eventually the amplitude of all the modes seem to fluctuate at about the same

amplitude. The evolution of the maximum power of one model is shown in Figure 2-6.

The same behavior is initially seen in the Mc/MD = 1 disks, but a slowly growing

m=1 mode eventually dominates most of these disks, with one exception. Unlike the

other models, the dominating m=l mode is not a spiral wave, but instead is the result

of the central object and disk rotating about their center of mass. That is, the "central"


object, representing the star, is no longer centered on the mass distribution of the system.

Another characteristic of the modal evolution of these disks is that an m=2 mode will

dominate before this m=l mode becomes the dominant feature. Figure 4-5 shows a plot

of the maximum power reached by modes I through 4, between a radius of 0.1 and 1.0,

in the simulation of model B2.

It should be noted that characterizing the evolution by the maximum power present

in each mode may have limitations. Because the Fourier analysis was done on the

normalized density, the outer portions of the disk will primarily be represented, for this

part of the disk possesses the greatest density contrasts.

This transient behavior, as well as the number of modes present, does not allow the

frequencies to be unambiguously identified with a particular mode. The presence of many

initial transient modes is not surprising, as many modes will be excited by the perturbation

that has been seeded into the initial conditions. However, the persistence of this transient

behavior in a dissipative system is unexpected. In addition, the details of the evolution

of the modes for any given disk seem to be sensitive to initial conditions, in the sense

that a change in a computational parameter, such as the tolerance angle or the number of

particles, will cause the disk's evolution to diverge in detail from a corresponding model

that is otherwise the same. This chaotic behavior brings into question the practicality of

describing in detail the complex evolution of the modes that are observed, though general

features of the evolution have been noted. (This same sensitivity to initial conditions is

seen in meteorology, where dissipation also is present.)


This behavior can be explained in two ways; the fluctuations are either reflecting

the physical behavior of such a system, or they are due to the failure of the numerical

method. If the behavior is physical then two processes may be occurring either separately,

or in conjunction. One process is the interaction between modes, where energy is being

transferred between modes continuously. In the disks that are modeled here, interaction

between the modes is enhanced because of the proximity of the corotation of a mode m

with the outer Lindblad resonance of the m 1 mode. and the inner Lindblad resonance

of the m+l mode, in the frequency radius domain. In Figure 4-6 the location of the

resonances are shown. In such a system energy is easily transferred from one mode to

another. The other possible process is that the power in a mode is fluctuating due to

the presence of two modes with different frequencies, but the same mode number, which

would result in the oscillation of the amplitude of that mode.

Alternatively, the cause of the fluctuation in power may be due to a failure of the

SPH method. An example would be if the method could not describe a shock front

properly once it had developed a large amplitude. In this case the post-shock oscillations

could destroy a wave mode that has developed a shock front. There may be additional

problems in the outer portions of the disks where the number density of SPH particles can

become too low to describe regions of low density. It is in this part of the disk that the

highest normalized density amplitudes appear. Gaps in the particle distribution appear

between the arms of the non-axisymmetric modes; these gaps are low density regions

that are under-represented. To see if the transient character of the modes resulted from

such problems, models were run with 10093 particles. This increased the resolution of


the method, as the smoothing length h is inversely proportional to the square root of the

number density. While the under-represented regions of the disk were further out in the

disk than in the models with 5130 particles, the same transient behavior of the modes

was observed, though differing in detail. This suggests that the transient behaviour of

modes is an inherent effect found in such disks, rather than being a numerical effect.


The rate of accretion onto the central object is determined by the rate of angular

momentum transport effected in the disk by shear viscosity and the non-axisymmetric

modes. In most of the models simulated, the mass accretion rate typically begins at a

large value, and within 27r dynamical times approaches a constant rate, which is more

sensitive to the viscosity parameter than to the ratio of specific heats. In some cases a

constant accretion rate is immediately realized (models A3 and C3). The mass accretion

is shown in Figures 4-7 through 4-11 for all models. The constant accretion rates in the

models were measured by fitting a line to the appropriate section of the mass accretion

curves, as shown by one example in the latter figure. Not surprisingly, the models with

Mc/MD = 1 show higher constant accretion rates than the Mc/MD = 3 disks. Although

it is not clear why constant accretion rates should be an endemic feature of accretion

disks, they are seen to some degree in most of the simulated disks, with only models

Al and Cl showing a small continual decrease of their accretion rate with time. For

these two cases the accretion rates measured are the asymptotic accretion rate for most

of the simulations.


Interestingly, though a constant accretion rate is seen in almost all models, it is not

necessarily maintained throughout the run; in some cases the accretion rates change to

new constant values, sometimes with almost discontinuous abruptness. Also worthy of

note is the insensitivity of the initial accretion rates of models DI and D2 to viscosity

(a = 0.25 and 0.5), as the mass accretion of these models is nearly equal. It is after new

constant accretion rates are re-established that the rates of these two models diverge. This

behavior is also observed in models C 1 and C2, and therefore seems to be characteristic of

the Mc/MD = 1 models with lower values of a. With the exception of model D3, all of

the MC/MD = 1 models show at least one accretion rate shift, as well as the Mc/MD = 3

models A3 and C3, which both have a viscosity parameter equal to one. Two models

show a second accretion rate change: models A3 and Dl. For the Mc/MD = 1 models,

the first accretion rate change takes place at a transition point in the modal evolution of

the disk. The maximum power in the m=2 mode peaks preceding the change in accretion

rate, and the m= mode becomes the dominant mode in the disk. The second accretion

rate change seen in model Dl takes place when the power in all modes with m > 1

abruptly falls to 0.2, one third the value of the m= 1 mode.

The exception to this behavior seen in the Mc/MD = 1 models is model D3, which

does not develop a dominant m=l mode. In this model the accretion rate changes very

slowly until a constant accretion rate is established. Unlike the Mc/MD = 1 models, the

constant accretion rate changes observed in models A3 and C3 are not clearly matched

with identifiable events in their modal evolution. In Table 3 the accretion rates, and the

times that they are attained, are given for all of the models in dimensionless units. To


convert the accretion rates to physical units, multiply the rates by the conversion factor

(MT/TD), the total mass divided by the dynamical time in the desired units. For instance,

using MT = I1M. and the dynamical time TD = 159 years, the accretion rates are found

to range between 2.36 x 10-6M,/yr and 4.31 x lO-5MoV/yr.

Table 4-3: Accretion Rates (x l0-3)

Model ril t rm t2 rh3 t3
Al 1.00 NA
A2 0.813 8
A3 1.13 0 0.688 8 0.375 20
A4 0.670 10
B1 3.88 7 6.11 18
B2 4.00 7 5.00 20
B3 3.00 9 7.32 22.5
C1 0.813 NA
C2 0.531 7
C3 0.875 0 0.375 20.5
DI 4.25 6 6.86 18 4.63 27
D2 4.51 4 3.05 26
D3 2.13 17

More surprising is the dependence of the mass accretion on viscosity, for the mass

accretion increases as a decreases. This stated dependence is evident in Figures 4-7

through 4-10, though it is not established in the Mc/MD = 1 disks with lower viscosity

parameter values until the constant accretion rate changes. The inverse dependence of

mass accretion upon shear viscosity is reflected in the mass accretion rates, which are

plotted in Figures 4-12 and 4-13. For the Mc/MD = 3 models the trend is established


after the secondary accretion rates have been attained in the a = 1 disks, which are the

rates that coincide temporally with the mass accretion rates observed in the models with

smaller viscosity parameters. It is also the secondary accretion rates in the Mc/MD = 1

models that show the same trend with shear viscosity, with the exception of the secondary

accretion rate of model B3.

This counter intuitive relationship between viscosity and mass accretion occurs

due to the combined actions of the global non-axisymmetric modes and local viscous

processes, which are not independent of one another. When both are present the viscosity

damps the more efficient mechanism of angular momentum transport provided by the

non-axisymmetric modes. The numerical experiments show that lower effective shear

viscosity allows the non-axisymmetric modes to attain greater strength, which in turn

become more effective at transporting angular momentum. This damping action of the

viscosity is not only qualitatively apparent, but is also evidenced by the rms amplitudes

of the density fluctuations.

In all the models, the total rms amplitude of the normalized density initially increases

exponentially with time. Most models then undergo a period of linear growth until a

saturation level is reached by the model; model D3 attains a saturation level immediately

after the initial exponential growth. Two models, DI and D3, also show a second

saturation level being attained later in the simulation. In addition to this general behavior

of the normalized amplitudes, erratic fluctuations in its magnitude are apparent that take

place on a time scale of about two dynamical times. The fluctuations grow in amplitude

until the saturation level is reached, and are not obviously present until the period of


linear growth begins. The shear viscosity suppresses the rms amplitude of the density

fluctuations, causing lower saturation levels to be attained. Also, the viscosity generally

enhances the rate of the initial exponential growth.

Table 4-4: Normalized density amplitude growth rates and saturation levels.

Model Viscosity (a) Growth Rate Growth Rate Saturation
(at R = 0.55) (at R = 1) levels)
Al 0.25 0.7 1.7 0.7
A2 0.5 1.3 1.8 0.55
A3 1 1.7 2.1 0.5
A4 0.5 1.3 1.8 0.5
BI 0.25 1.1 1.5 1.1
B2 0.5 0.9 1.5 0.9
B3 1 1.4 2.1 0.8
C1 0.25 0.7 2.0 0.9
C2 0.5 0.9 2.1 0.8
C3 1 1.5 2.6 Not attained
Dl 0.25 1.0 1.8 0.65
D2 0.5 1.0 2.2 1.2 (0.7)
D3 1 1.2 3.0 0.3 (0.45)

The growth rates and saturation levels of the rms amplitude of the normalized density

fluctuations also vary with radius. The saturation level of the rms amplitude increases

with radius, as do the growth rates. Table 4 shows the growth rates for the various

models at the radial distances of 0.55 and 1.0 from the center of mass of the system. The

saturation level given corresponds to the smaller radius of 0.55.


The method of treating the central region and accretion in the hydrodynamic code

was intended to let the properties of the disk drive accretion, rather than be determined

by the imposed inner boundary condition around the central object. To test the degree to

which I am successful, models were run with two different inner disks (models A2 and

A4), and another with a larger Ra (model A5). Both types of inner disks, an exponential

disk and a massless Keplerian disk, provide pressure support at the inner boundary by

maintaining a continuous surface density and velocity field across the boundary (see

Chapter 3 for further details). Models A2 and A4 are identical except for the treatment

of the inner region, and the mass accretion in the two models is nearly the same (see

Figure 4-11). However, since the massless Keplerian disk results in the density gradient

being discontinuous at the boundary, the inner exponential disk was preferred. To what

degree the size of Ra alone influences accretion is tested by model A5, which is identical

to model A4 excepting the size of the inner region. The value of Ra for model A5,

0.0635, deviates from the nominal value of 0.0544. As can be seen in Figure 4-11, the

size of the inner accretion region does not significantly affect the mass accretion.

The initialization of an encounter model is described in chapter 3, which includes an

encountering particle with a mass of 0.5, and a disk that is otherwise equivalent to model

D2. In the resulting simulation it is found that the encounter takes place approximately

after ten dynamical times. The time sequence of this encounter is shown in the following

chapter, in Figures 5-9 through 5-12. Beside robbing a significant fraction of mass from

the disk, the encountering particle also greatly enhances the mass accretion onto the

central object. This is shown in Figure 4-14.


Tests of Code Parameters

As discussed at the beginning of this chapter, models with differing values of code

parameters, but equivalent model parameters, should be compared to confirm that the

choice of code parameters yield consistent results. Since it was impractical to run such a

series of tests for all combinations of the physical parameters, a single model was chosen:

model D2. The test models are designated TI, T2, and T3. In model TI the tolerance

angle, 0, is given a value of 0.5, resulting in a gravitational approximation closer to a

direct summation of the gravitational terms due to .V particles. In model T2 the Courant

parameter, C, is equal to 0.2, resulting in smaller integration time steps. In model T3

the number of particles used in the simulation is increased to 10093 particles. The mass

accretion of these three test models are shown with the mass accretion of model D2 in

Figure 4-15.

Model TI, with a smaller tolerance angle than that of D2, is nearly equivalent to

model D2, indicating that the gravitational forces are being calculated to a sufficient

degree of accuracy with 0 = 0.7. Model T2 initially has the same mass accretion as D2,

but then deviates after twenty dynamical times. This may indicate that the hydrodynamic

calculations require smaller time steps to yield sufficiently accurate results.

Model T3, which possesses 10093 particles, differs most dramatically from model

D2. However, these two models are not physically equivalent for several reasons. By

increasing the number of particles by nearly a factor of two the local smoothing length is

decreased by approximately a factor of 2-1/2. This effectively changes the local viscosity

of the disk. An additional difference in model D2 and T3 is also seen when the modal


evolution of model T3 is inspected: a m= 1 mode does not eventually dominate the disk

of model T3 as it does model D2. This doubtless accounts for the smaller mass accretion

rate in model T3 in spite of the fact that, with a smaller smoothing lengths, the effective

local shear viscosity is smaller than in model D2. As mentioned above, weaker shear

viscosity typically leads to larger accretion rates. The reason that model T3 does not

develop a dominating m= mode is most likely due to the fact that its initial density

perturbation also differs from that of model D2. This is because the density perturbation

was initially seeded into the disk by displacing the particles from the positions required

to describe a smooth exponential disk. These displacements are random in direction and

Gaussian in magnitude, with a standard deviation that is a fraction of the local interparticle

spacing. Because model T3 has a smaller interparticle spacing than model D2 the initial

noise in density in the two models are not the same. Hence, as models T3 and D2 are

not equivalent physically, these two models can not be directly compared.

\,* ^*^-
-l -1 0 I

-" -I 0 1 2

0' ;, I ++ 2

L y^ -*.'~.* **! -.**

1 + .

. . .
-a -i 0 1 2

+ .. ."

-2 -1 0 1 2

-2 -1 0 1 21


Figure 4-1: Particle distribution of Mc/MD = 3 models at Time = 10.0 dynamical times.


;** ,. .- >

i .. ...^ i..'.. i .... i:

-* -I o I -
S .

: .

-1 0 I <


.' : :' '" "' *-

.. .*Jt ','.

0,- I 0 *

-z -I 0 1 2

^ .. ^.. r

-2 -1 0 I 2
]- *^^ iii ^7

K**' w

". .-OL a


Figure 4-2: Particle distribution of Mc/MD = 3 models at Time = 29.5 dynamical times.

-, : ,W I

- -I 0 I

M,/MD= 1

Figure 4-3: Particle distribution of Mc/MD = 1 models at Time = 10.0 dynamical times.


2, 0 ".

-2 0 I 2

- -

-2 -1 1 a ll I

^ .* ,. '

<:,, o
i ..* i .. *. i
- -I 0 I 2

-z -i o I 2

S .

-l 0 I I I "o


-2 -I 0 I 2
"!-.*it" 0

^^ ,.^ ,,-

' i -- i-.. ---i* i .

;': -" ,* ,ill '.

-I' -'I .. + -, i '


Figure 4-4: Particle distribution of Mc/MD = 1 models at Time = 30.0 dynamical times.






0 10 20 30

Figure 4-5: Maximum power in modes I through 4 for model B2.


0 I I I ,
0 0.2 0.4 0.6 0.8

Figure 4-6: Resonances in an accretion disk with Mc/MD = 3. Solid lines
correspond to the corotation resonance, the dotted lines to the Inner
Lindblad resonance, and the dashed lines to the Outer Lindblad resonance.

0 10 20 30 40

Figure 4-7: Mass accretion, as a fraction of initial disk mass, for models Al through A3.

0.5 .






0 10 20 30 40

Figure 4-8: Mass accretion, as a fraction of initial disk mass, for B models.


Mass accretion, as a fraction of initial disk mass, for C models.

,0.3 / .- -



0 10 20 30 40

Figure 4-10: Mass accretion, as a fraction of initial disk mass, for D models.

Figure 4-9:





0 10 20 30 40

Figure 4-11: Mass accretion in models A2, A4, and A5, as a fraction of the initial disk mass.

0 0.2 0.4 0.6

0.8 1 1.2

Figure 4-12: Constant mass accretion rates for Mc/MD = 3 accretion disks.

0.01 1 1 i *- i 1 1

0.8 1 1.2

Figure 4-13: Constant mass accretion rates for Mc/MD = 1 accretion disks.


0.0015 r







0 0.2

0.4 0.6






0 10 20 30 40

Figure 4-14: Mass accretion of the encounter model (solid line)
compared to the mass accretion of model D2 (dotted line).

Figure 4-15: Mass accretion, as a fraction of initial disk
mass, for model D2 is shown with three test models.


The greatest difference between disks with higher and lower star to disk mass ratios

is that all disks that develop a dominating m= 1 mode also form satellites. This includes

all disks with Mc/MD = 1, with the exception of model D3. By way of example, the

later evolutionary sequence of model D2 is shown in Figure 5-1 and Figure 5-2. Figure

5-3 shows the final configuration of the disk, with the path of the satellite shown since

its formation. The satellite forms from a clump of gas that is recognizable as a distinct

structure at time = 25. At this time it is part of a spiral arm of a m=1 mode that extends

to large radii. The mass in the initial clump of gas is 0.013MT, and at time = 39 the

satellite has a mass of 0.028MT. Also, as the mass of the satellite increases, and the

combined mass of the star and the disk interior to the companion decreases, their mean

separation decreases. A closer view of the satellite is shown in Figure 5-4, which also

shows its radial density profile.

Other models form multiple satellites, which can interact with one another. Table 2

provides an overview of the formation and evolution of the satellites, giving the initial

separation Ri between the mass and the central object disk center of mass, the initial

mass Mi, the time of formation Ti, the final separation Rf, and the final mass of the

satellite, Mf. The final time, Tf, corresponds to either the end of the simulation, or to

when the satellite was destroyed.

Table 5-1: Satellites

Model Satellite Ri Mi Ti Rf Mf Tf
(x 10-2) (x 10-2)
BI 1 1.96 0.448 31.5 3.25 0.536 38.8
2 1.57 0.409 22.5 4.72 0.692 38.8
3 2.01 0.370 23.0 4.57 0.419 38.8
B2 1 2.23 0.546 24.0 3.72 1.277 38.5
2 1.68 0.975 24.0 0.81 4.045 37.5*
B3 1 1.33 1.374 21.0 1.48 4.776 31.0
2 1.41 0.760 24.0 2.26 1.72 31.0
DI 1 1.21 0.770 18.5 2.31 1.920 33.0
2 1.17 1.150 17.5 0.98 2.174 33.0
3 1.30 0.741 22.0 1.65 0.955 33.0
4 1.17 1.053 19.0 0.67 1.803 24.5*
D2 1 1.50 1.257 25.0 1.24 2.797 39.0

* Indicated satellites are reabsorbed into the disk.

Orbital Evolution

In the absence of close encounters with other satellites, the orbital evolution can be

understood by considering the simpler system of two point masses in circular motion

about a common center of mass. In such a two body system the separation of the two

masses is determined by the masses and the orbital angular momentum, given explicitly by

L2 m
a G (mG m),


where L is the orbital angular momentum about the center of mass, and m = m1 + m2

is the total mass. The above equation shows how the equivalent separation between


two bodies can decrease, assuming they remain on circular orbits. If the total mass

and orbital angular momentum remain constant, and mass is transferred from the more

massive body (the primary), mi, to its less massive companion (the secondary), m2,

then the separation will decrease as 1/(mim2)'. The separation will also decrease if

the orbital angular momentum decreases. A change in the orbital angular momentum

can occur by a variety of mechanisms: a) external torques on the system, such as an

outer satellite, b) internal torques, which can transfer orbital angular momentum into,

or out of, the spin angular momentum associated with an extended asymmetrical mass

distribution about one of the bodies, c) mass loss from the system, which will carry away

orbital angular momentum with it, and conversely, d) mass accretion onto the system.

While mass loss from the system decreases the orbital angular momentum, it may cause

the separation to increase. From the above equation, and assuming that the mass loss

process does not alter the velocities of the two masses, and that the secondary's mass

remains constant, the condition for increasing the separation with mass loss from the

primary is m1/m2 > 1(1 + VTT/) 2.56.

To understand the orbital evolution of the formed satellites, the above model is applied

as an approximation. The mass of each of the satellites, and the mass within the minimum

separation between the satellite and central object, is found. These masses are taken to

be the masses of the primary and secondary respectively. Hence, the primary mass m,

consists of the central object and the mass in the disk, as well as any inner satellites, that

lie within the minimum separation. In addition, the orbital angular momentum associated

with these two masses, with respect to their center of mass, is found. Using the above


formula (equation 5.1), the equivalent separation, a, of a system in circular motion, with

the same orbital angular momentum and masses, was found. Often the formed satellites

exhibit substantial eccentricity, which results in the actual separation oscillating about the

equivalent separation. All of the above mechanisms for altering the separation between

the disk-star and the satellite are observed to some degree. Mass is commonly lost

from the disk as mass is transported outward with angular momentum. As the primary

mass of the disk-star far exceeds the mass of the satellite, this mass loss results in the

separation increasing, if the mass of the satellite remains constant. Interaction with the

non-axisymmetric modes of the disk, which act as a reservoir of angular momentum,

is also important. Though this model is an approximation it gives sufficient insight to

understand the orbital evolution of the satellites.

In the case of model D2, where only a single satellite is present, the decrease in the

separation is primarily due to a change in the mass ratio of the primary and the satellite. In

model B2 (Figure 5-5) two satellites form almost simultaneously at the same polar angle.

The outer satellite then spirals outward, while the inner satellite spirals inward. Initially

this can be understood as resulting from the interaction of the two satellites. However,

while the outer satellite's orbital angular momentum rises sharply soon after formation,

a comparable drop in the inner satellite's orbital angular momentum is not observed; the

orbital angular momentum of the inner satellite initially remains constant. Therefore any

angular momentum transferred from the inner to the outer satellite must be offset by an

equal transfer of angular momentum from the disk to the inner satellite, or else angular

momentum is being transferred from the disk directly to the outer satellite. Later the


angular momentum of the inner satellite increases dramatically, due to interaction with

the non-axisymmetric modes in the disk. However, this increase of the inner satellite's

orbital angular momentum is not enough to counter the effects of mass transfer from the

disk to the satellite, and it continues to spiral inward to eventually be reabsorbed into

the disk. At the same time the outer satellite continues to slowly gain orbital angular

momentum, even while its own mass, and the mass interior to it remains constant. The

source of this orbital angular momentum must be from the disk and satellite interior to

it, as no external torques are present.

In model B3 two satellites also form (Figure 5-6). However, due to a different

configuration of the satellites, they do not strongly interact. The innermost satellite

remains at nearly a constant separation for four dynamical times, then gently increases

to its maximum peak value at T=28.5, and then diminishes to its final value. Though the

separation does not vary greatly, the mass and orbital angular momentum show abrupt

increases, making "steps" to higher values. These steps are correlated with collisions of

the satellite with arms of a mode that extends out to the orbital radius of the satellite and

overtakes it. Between these collisions the separation between the satellite and primary

continues to decrease due to mass loss from the primary. In contrast to the inner satellite,

the outer satellite's separation grows continuously after its formation. Its mass and orbital

angular momentum also remain constant for four dynamical times, after which both

increase as it encounters material moving outward from the disk. The continual change

in the separation is again due to the decrease of the primary's mass.


In model Bl three satellites form (Figure 5-7), and all move outward to greater

radii. Soon after their formation, the three satellites' masses remain constant. The orbital

angular momentum of the first two satellites increase nearly linearly with time during

the simulation, while the outermost satellite's orbital angular momentum remains nearly

constant during its trajectory. Again, the increase in separation is due largely to the

primary mass, that corresponds to each satellite, decreasing in magnitude. However, the

net increase in orbital angular momentum in this system of bodies, despite the loss of

mass to the primary, is clear evidence of angular momentum being transferred from the

disk, which has a strong m=l mode, to the satellites.

Model DI, which I have left for last, is the most complicated, as four satellites

are formed (Figure 5-8). The first two satellites are clearly self-gravitating and distinct

structures. A third is much more diffuse, but remains as a recognizably distinct structure

for eleven dynamical times, perhaps with the aid of the first two satellites, which are at

larger radii. The fourth satellite is short lived, encountering the first satellite soon after

its own formation, and is consequently reabsorbed into the disk. Until just previous to

this collision the first satellite's mass and angular momentum are constant with time, but

both increase abruptly during the collision. Afterward, the orbital angular momentum

of the first satellite increases linearly with time, presumably interacting further with the

fourth satellite, as well as with the m=l mode of the disk. The mass of the second

satellite remains constant with time soon after its formation, while the orbital angular

momentum initially increases, then levels, and gradually decreases. The separation also

increases, then decreases. This is due not only to ellipticity but also the loss of orbital

Full Text
xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd