UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
NUMERICAL SIMULATION OF VISCOUS ACCRETION DISKS By RONALD EUGENE DRIMMEL A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1995 For my parents, for the poor, and to Sacred Truth and Beauty. ACKNOWLEDGMENTS 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 Nbody 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. TABLE OF CONTENTS ACKNOWLEDGMENTS ............ LIST OF TABLES ................ . . . iii . vii LIST OF FIGURES .......................... ABSTRA CT .......................... ... .. ... ... .... xi CHAPTERS 1. INTRODUCTION .................................... 1 2. METHODS ........ .... ................ Development of the Hydrodynamic Code ........ Nbody 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 ................ TwoDimensional 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 3. INITIALIZATION OF DISKS ............ Twodimensional Disks ............... M aclaurin Disk ................. Exponential Disks ............... Inner Disks ....... ............. 4. SIMULATION OF TWODIMENSIONAL DISKS Param eters ...................... Evolution of Disks ................. General Features ................ Modal Evolution ................ Accretion .................... Tests of Code Parameters ............. . . . 56 . . . 56 . . 60 . . 60 . . 62 . . 65 . . . 7 1 5. FORMATION OF SATELLITES ...... Orbital Evolution .............. Formation Conditions ............ Encounter Model .............. 6. SUMMARY AND CONCLUSIONS . 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 BIOGRAPHICAL SKETCH . 281 LIST OF TABLES Code Parameters ............... M odels . . . Models at T=28 dynamical times ...... Accretion Rates (x 103) .......... Normalized density amplitude growth rates Satellites ........... ........ ........and saturation levels... and saturation levels. . . . Table 21: Table 41: Table 42: Table 43: Table 44: Table 51: LIST OF FIGURES Figure 21: Figure 22: Figure 23: Figure 24: Figure 25: Figure 26: Figure 27: Illustration of a twodimensional 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 twodimensional 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 31: Initial distribution of particles for the Maclaurin Disk. ... 41 Figure 32: The initial Toomre Q parameter (solid line) and tangential velocity (dotted line) with respect to radius for a Mc/MD = 3 disk. .. 46 Figure 33: Initial distribution of particles in an exponential disk. . 48 Figure 34: Figure 41: 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 42: Particle distribution of Mc/MD = 3 models at Time = 29.5 dynamical tim es. . .. . . 74 Figure 43: Particle distribution of Mc/MD = 1 models at Time = 10.0 dynamical times. ......... ..................... 75 Figure 44: Particle distribution of Mc/MD = 1 models at Time = 30.0 dynamical tim es. . .. . . . 76 Figure 45: Maximum power in modes 1 through 4 for model B2. ... 77 Figure 46: 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 47: Mass accretion, as a fraction of initial disk mass, for models Al through A 3. . . . . 79 Figure 48: Mass accretion, as a fraction of initial disk mass, for B models. .. 79 Figure 49: Mass accretion, as a fraction of initial disk mass, for C models. .. 80 Figure 410: Mass accretion, as a fraction of initial disk mass, for D models. 80 Figure 411: Mass accretion in models A2, A4, and A5, as a fraction of the initial disk m ass. . . . . ... 81 Figure 412: Constant mass accretion rates for Mc/MD = 3 accretion disks 82 Figure 413: Constant mass accretion rates for Mc/MD = 1 accretion disks 82 Figure 414: Mass accretion of the encounter model (solid line) compared to the mass accretion of model D2 (dotted line). . ... 83 Figure 415: Mass accretion, as a fraction of initial disk mass, for model D2 is shown with three test models. ........................ 83 Figure 51: Evolutionary sequence of model D2 showing formation of satellite.. 96 ix Figure 52: Evolutionary sequence of model D2 continued. . ... 97 Figure 53: Final configuration of model D2, showing the position of the satellite at earlier times. ............. .................. 98 Figure 54: Radial density profile and particle distribution of satellite D21. .99 Figure 55: Final configuration of model B2, showing the position of the satellites at earlier tim es .. ... .. .. .. .. .. .. .. .. 100 Figure 56: Final configuration of model B3, showing the position of the satellites at earlier tim es .. .. . . 101 Figure 57: Final configuration of model Bl, showing the position of the satellites at earlier tim es .. .. .. .. .. .. .. ... .. ... .. 102 Figure 58: 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 59: Encounter model ................................ 104 Figure 510: Encounter model continued. . . ... 105 Figure 511: Encounter model continued. . . ... 106 Figure 512: 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 NUMERICAL SIMULATION OF VISCOUS ACCRETION DISKS By RONALD EUGENE DRIMMEL August, 1995 Chairman: James H. Hunter, Jr. Major Department: Astronomy A two and threedimensional 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 selfgravitating 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 nonaxisymmetric 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. CHAPTER 1 INTRODUCTION 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 selfgravity 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 I 2 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, disklike 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 "flattop" 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 JCMTCSO 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 selfgravitating 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 3 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 (LyndenBell & 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 magnetorotational 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 selfgravity of the disk is neglected, and only pressure supports the vertical structure. In contrast, the theory of selfgravitating 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 selfgravitating accretion disks began with Nbody 4 modeling (Cassen et al. 1981; Tomley et al., 1991). Papaloizou and Savonije (1991) analytically investigate the instabilities that exist in twodimensional selfgravitating 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 selfgravitating threedimensional disks been investigated through numerical simulation. The stability of threedimensional 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 threedimensional 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 5 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 selfgravitating 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. CHAPTER 2 METHODS 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. Nbody Seed Code The code developed to evolve accretion disks is itself evolved from a version of Hernquist's Nbody 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 21). 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 21: Illustration of a twodimensional 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 N1 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, SN A(2.1) 1 and the absolute average acceleration, _ arect (2.2) i He define the mean error as 1 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 N (ba oai)/arec (2.5) j=1 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 Nbody 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; 10 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 Nbody 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. SPH Using this Nbody code as a skeleton, a twodimensional 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 threedimensional version of this hydrodynamics code was then developed and dubbed T3DSPH. 11 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 n(r) = 6(r rj), (2.9) j=1 then the integral expression for the estimate (f) (equation 2.3) can be written as a sum: N (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 N 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: N (pA)i = WmAjWi. (2.13) J 13 This convention is considered a "golden rule of SPH" implementation by Monaghan (1992). 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 estimated. 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 14 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). 15 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 Nbody 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) where ( P) smooth P d[ = P( v. v), 2V" )i'" (2.19) and (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: dri dt p VP + asc VPi, dt p, 17 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) where 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 18 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 nonadiabatic viscous heating, and the final terms represent other nonadiabatic 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 Courantlike 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 ' 20 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 Nbody codes. Real gravitating systems often possess orders of magnitude more gravitating particles than can be modeled directlythe computational requirements are simply 21 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 Nbody code, is often as many orders of magnitudes as the number of particles. When using an Nbody 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: mimjrij 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.), (2.35) 22 from Poisson's equation. Since the operator (V) = Irr2 in spherical coordinates we have 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) 0 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. 23 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. Accretion During the hydrodynamic simulation an inner region, within a radius Ra of the central particle, is treated semianalytically 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 nonaxisymmetric 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 24 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 (340), 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. 25 Summary of Code Parameters A list of the parameters which control the performance of T2DSPH is listed for convenience in Table 21. The first parameter, N, determines how well a continuous Table 21: 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 26 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. Tests 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 twodimensional 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 selfgravity, is a simulation of a Maclaurin disk, the twodimensional 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 22 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 R Figure 22: Surface density and velocity profile of the Maclaurin disk after 26 dynamical times. TwoDimensional Shocks Twodimensional adiabatic shocks, without selfgravity, 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 twodimensional 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 28 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 postshock 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/ nonadiabatic 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 postshock 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 where b = 2vo + Q. 2 L?'o The solution above is given in the postshock reference frame, with vo and v2 representing the preshock 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~~00 a 0 I I Figure 23: Surface density of three twodimensional 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. I I I (2.45) (2.46) "( 30 In Figure 23 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 31 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 byproduct 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 selfgravity. 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. 32 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 r v = (Er )' IVrEr(r 2)'dr. (2.48) 0 For a Keplerian potential f22 = GM/r3, which upon substitution results in r v = ( Er 1/2) I / vrdr. (2.49) 0 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 33 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 24). 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 34 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 25). 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 26). 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 27). 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.8 0.6 0.4 0.2 0 2 4 6 8 10 Mode INT=0.05, BOT=0.01 TIME=21.0055 Figure 24: 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  0.2 10 15 20 25 30 Time INT=0.05, BOT=0.05 Figure 25: 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 Time Figure 26: Maximum of the power spectrum, in modes = 1 4, of the normalized density in model A2. 38 15 I I I 10 5 ^ ^ 0 0 0.2 0.4 0.6 0.8 Radius BOT=13.2012=99.9% cont. Time = 20.0 to 30.0 Figure 27: The average Lombe periodogram of normalized density temporal fluctuations along Cartesian axes for model A2. CHAPTER 3 INITIALIZATION OF DISKS 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. 40 Twodimensional Disks Maclaurin Disk A twodimensional disk has already been presented in Chapter 2 as one of the tests to the SPH code: the Maclaurin disk. This is a twodimensional 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/ 41 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 31). 1 0.5 h > 01 1 0.5 0 0.5 1 x Figure 31: 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 42 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 31, 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= (1k2) RD (3.5) 91MD 43 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 44 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 eRDtl+ ~. (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 NewtonRaphson method, with the expression for the mass that lies interior to a given radius for the exponential disk: m(r) = Eo27rr [1 er/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 rn1)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 n there are Nr = E Nj particles. This iterative procedure is continued until Nr > Ntot. If j=1 Nr > Ntot then Ntot is set equal to Nr, the particle mass mrn is reinitialized, and the disk 45 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 r function m(r) = 27 f E(u)udu be specified, it is not known whether the iteration will 0 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 twodimensional 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 46 evaluated for each ring n with the expression 1 (V2(rn+) V'2(r)) V (r) n= 2 (3.11) rk ((n+1 rn1) 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 equation: K (Qmin7rG)2  (3.13) / r37y 7mink Kkk For all the twodimensional exponential disks, Qmin is set to 1.15. In Figure 32 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 32: The initial Toomre Q parameter (solid line) and tangential velocity (dotted line) with respect to radius for a Mc/MD = 3 disk. 47 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 33, while the initial radial profile of the density is shown in Figure 34. 0.5  > 0 0.5 1 1 0.5 0 0.5 1 X Figure 33: 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 " 2 0 0.2 0.4 0.6 0.8 1 1.2 R Figure 34: 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 onedimensional 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) 50 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 NewtonRaphson method. To find the associated perturbed velocities, the velocities of the onedimensional 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 expression V2() = V02 7 K S 1 p(r)cos(my)(1 + ( )cot[ r(rr }. (3.17) 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 51 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) with 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 = WFE. 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 52 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) = Eoer/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 53 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 NewtonRaphson 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 NewtonRaphson 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 = rn1 Arm_1/2 is used. 54 For clarity's sake, I restate the algorithm for finding the radii of the concentric rings: 1. Using the NewtonRaphson method, find r, which satisfies the condition 2(rL rn) = (mp/E(rn))1/2, where rL = Ra for n = 1, but rL = r,n1 Arn_1/2 otherwise. 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 NewtonRaphson 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 55 mn. This will allow us to circumvent the use of the NewtonRaphson method in step 6) above, and instead use the expression: 2mr rM = 2 (3.24) B + /1B2 + 4(7r(Ra) B/2Ra)m, where B = rGMc(1 k2)/K, and mr = m(RL) Amk. CHAPTER 4 SIMULATION OF TWODIMENSIONAL DISKS Parameters 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 selfgravitating, 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 twodimensional 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. 57 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 58 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 41: Models a 0.25 0.5 1 0.5 0.25 0.5 1 0.25 0.5 1 0.25 0.5 1 Inner Disk Keplerian it Exponential Exponential Exponential Exponential 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 21, with their nominal values. Ideally, if the values Mc/MD 3 Model: Al A2 A3 A4 BI B2 B3 Cl C2 C3 DI D2 D3 Tf 33.10 32.11 33.01 24.01 38.81 38.51 31.01 32.61 29.97 29.55 33.01 39.00 32.63 59 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 interparticle 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 interparticle separation in two dimensions is n1/2, with n representing the local number density. An average interparticle distance is defined as n1/2 = V/7rRD/N, which is equal to 0.025 for RD = 1 and N = 5131. The condition e > fi1/2 is satisfied with e = 0.027, but e < n1/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 60 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 42 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 41 and 43, and after thirty dynamical times in Figures 42 and 44. 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 44, 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 61 from the gravitational interaction between the massive central object and the disk, and has been predicted by others (Savonije et al., 1992). Table 42: Models at T=28 Accreted mass (102) 3.504 2.851 2.042 15.26 14.23 11.69 2.856 1.808 2.125 16.34 13.47 7.466 dynamical times MIc/IMD 3.790 3.645 3.477 1.985 1.895 1.697 3.646 3.430 3.493 2.084 1.833 1.422 Transported AM (102) 2.501 2.923 2.562 10.68 9.877 11.67 2.844 2.808 2.440 14.15 9.461 8.323 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 interarm regions. At the later times, such as those shown in the Figures 43 and 44, 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 nonaxisymmetric 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. Models Al A2 A3 BI B2 B3 Cl C2 C3 Dl D2 D3 a 0.25 0.5 1.0 0.25 0.5 1.0 0.25 0.5 1.0 0.25 0.5 1.0 62 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 nonaxisymmetric 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 26. 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" 63 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 45 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.) 64 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 46 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 postshock 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 nonaxisymmetric modes; these gaps are low density regions that are underrepresented. To see if the transient character of the modes resulted from such problems, models were run with 10093 particles. This increased the resolution of 65 the method, as the smoothing length h is inversely proportional to the square root of the number density. While the underrepresented 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. Accretion 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 nonaxisymmetric 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 47 through 411 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. 66 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 reestablished 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 67 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 106M,/yr and 4.31 x lO5MoV/yr. Table 43: Accretion Rates (x l03) 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 47 through 410, 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 412 and 413. For the Mc/MD = 3 models the trend is established 68 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 nonaxisymmetric 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 nonaxisymmetric modes. The numerical experiments show that lower effective shear viscosity allows the nonaxisymmetric 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 69 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 44: 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 (<0.6) 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. 70 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 411). 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 411, 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 59 through 512. 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 414. 71 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 415. 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 21/2. This effectively changes the local viscosity of the disk. An additional difference in model D2 and T3 is also seen when the modal 72 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 M,/MD=3 Figure 41: Particle distribution of Mc/MD = 3 models at Time = 10.0 dynamical times. 7=2 ;** ,. . > i .. ...^ i..'.. i .... i: * I o I  S . : . 1 0 I < 7=5/3 .' : :' '" "' * .. .*Jt ','. 0, I 0 * z I 0 1 2 I J6 ^ .. ^.. r 2 1 0 I 2 ] *^^ iii ^7 K**' w ". .OL a Mc/M,=3 Figure 42: Particle distribution of Mc/MD = 3 models at Time = 29.5 dynamical times. , : ,W I  I 0 I M,/MD= 1 Figure 43: Particle distribution of Mc/MD = 1 models at Time = 10.0 dynamical times. y=5/3 I 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 7=5/3 2 I 0 I 2 "!.*it" 0 ^^ ,.^ ,, ' i  i.. i* i . ;': " ,* ,ill '. I' 'I .. + , i ' M,/MD= Figure 44: Particle distribution of Mc/MD = 1 models at Time = 30.0 dynamical times. 0.8 0 0.6 0.4 0.2 0 10 20 30 Time Figure 45: Maximum power in modes I through 4 for model B2. x 0 I I I , 0 0.2 0.4 0.6 0.8 radius Figure 46: 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 Time Figure 47: Mass accretion, as a fraction of initial disk mass, for models Al through A3. 0.5 . a0.25 0.4 0.3 0.2 0.1 0 10 20 30 40 Time Figure 48: Mass accretion, as a fraction of initial disk mass, for B models. Time Mass accretion, as a fraction of initial disk mass, for C models. ,0.3 / .  0.2 0.1 0 10 20 30 40 Time Figure 410: Mass accretion, as a fraction of initial disk mass, for D models. Figure 49: 0.2 0.15 0.1 0.05 0 10 20 30 40 Time Figure 411: Mass accretion in models A2, A4, and A5, as a fraction of the initial disk mass. 0 0.2 0.4 0.6 G( 0.8 1 1.2 Figure 412: Constant mass accretion rates for Mc/MD = 3 accretion disks. 0.01 1 1 i * i 1 1 0.8 1 1.2 Figure 413: Constant mass accretion rates for Mc/MD = 1 accretion disks. 0.002 0.0015 r A U 0.008 0.006 0 A A A U U A 0 0.2 0.4 0.6 a I I I I I I 0.4 0.3 0.2 01 0 10 20 30 40 Tine Figure 414: Mass accretion of the encounter model (solid line) compared to the mass accretion of model D2 (dotted line). Figure 415: Mass accretion, as a fraction of initial disk mass, for model D2 is shown with three test models. CHAPTER 5 FORMATION OF SATELLITES 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 51 and Figure 52. Figure 53 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 54, 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 51: Satellites Model Satellite Ri Mi Ti Rf Mf Tf (x 102) (x 102) 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), (5.1) 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 86 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 87 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 diskstar 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 diskstar 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 nonaxisymmetric 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 55) 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 88 angular momentum of the inner satellite increases dramatically, due to interaction with the nonaxisymmetric 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 56). 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. 89 In model Bl three satellites form (Figure 57), 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 58). The first two satellites are clearly selfgravitating 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 UTF8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd INGEST IEID EST0SV8GC_QOZJBT INGEST_TIME 20120220T20:59:35Z PACKAGE AA00009020_00001 AGREEMENT_INFO ACCOUNT UF PROJECT UFDC FILES 