Equilibrium Initialization and Stability of Three-Dimensional Gas Disks
Hsiang-Hsu Wang, Ralf S. Klessen, Cornelis P. Dullemond, Frank C. van den Bosch, Burkhard Fuchs
aa r X i v : . [ a s t r o - ph . GA ] A p r Mon. Not. R. Astron. Soc. , 1– ?? (2009) Printed 11 November 2018 (MN L A TEX style file v1.4)
Equilibrium Initialization and Stability ofThree-Dimensional Gas Disks
Hsiang-Hsu Wang , ⋆ , Ralf S. Klessen , , Cornelis P. Dullemond , Frank C. vanden Bosch , Burkhard Fuchs Max-Planck-Institut f¨ur Astronomie, K¨onigstuhl 17, D-69117 Heidelberg, Germany Zentrum f¨ur Astronomie der Universit¨at Heidelberg, Institut f¨ur Theoretische Astrophysik, Albert-Ueberle-Strasse 2, 69120 Heidelberg, Germany Zentrum f¨ur Astronomie der Universit¨at Heidelberg, Astronomisches Rechen-Institut, M¨onchhofstr. 12-14, 69120 Heidelberg, Germany Department of Physics and Astronomy, University of Utah, 115 South 1400 East, Salt Lake City, UT 84112, U.S.A. Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, Menlo Park, CA 94025, U.S.A.Accepted by MNRAS 2010 April 28
ABSTRACT
We present a new systematic way of setting up galactic gas disks based on the as-sumption of detailed hydrodynamic equilibrium. To do this, we need to specify thedensity distribution and the velocity field which supports the disk. We first show thatthe required circular velocity has no dependence on the height above or below themidplane so long as the gas pressure is a function of density only. The assumptionof disks being very thin enables us to decouple the vertical structure from the ra-dial direction. Based on that, the equation of hydrostatic equilibrium together withthe reduced Poisson equation leads to two sets of second-order non-linear differentialequation, which are easily integrated to set-up a stable disk. We call one approach‘density method’ and the other one ‘potential method’. Gas disks in detailed balanceare especially suitable for investigating the onset of the gravitational instability. Werevisit the question of global, axisymmetric instability using fully three-dimensionaldisk simulations. The impact of disk thickness on the disk instability and the forma-tion of spontaneously induced spirals is studied systematically with or without thepresence of the stellar potential. In our models, the numerical results show that thethreshold value for disk instability is shifted from unity to 0.69 for self-gravitatingthick disks and to 0.75 for combined stellar and gas thick disks. The simulations alsoshow that self-induced spirals occur in the correct regions and with the right numbersas predicted by the analytic theory.
Key words: galaxy: disc – galaxies: evolution – galaxies: structure – methods: nu-merical
The stability of gas disks plays an important role in govern-ing the structure of disk galaxies and in regulating their starformation rate. Although important insights can be obtainedusing perturbation theory (Toomre 1964, Lin & Shu 1964,Rafikov 2001), the onset of stability and its impact on thestar formation and evolution of gas disks is best studied us-ing hydrodynamical simulations. These can follow the non-linear behavior of the system, which cannot be addressedby linear analysis. With the recent advances in computing ⋆ E-mail: [email protected] power and the development of new numerical techniques,we are now in a good position to treat a three-dimensional,isolated galaxy self-consistently.However, in order for a stability analysis to be mean-ingful and reliable, it is of paramount importance that onecan specify equilibrium initial conditions. After all, if theinitial disk is not in equilibrium, its relaxation during thefirst time-steps of the simulation may trigger instabilitiesthat are of little relevance for our understanding of the sta-bility of disk galaxies. Unfortunately, no analytical solutionis known for the density, velocity field and temperature of athree-dimensional gas disk in hydrostatic equilibrium in theexternal potential of a dark matter halo and/or a stellar disk.Consequently, previous hydrodynamical simulations have ei- c (cid:13) Wang et al. ther started from non-equilibrium initial conditions, or haveresorted to iterative techniques to set-up the initial condi-tions, at the cost of having little control over the resultingequilibrium configuration. In this paper we present a newmethod that allows one to compute the density and veloc-ity structure of a realistic, isothermal, three-dimensional gasdisk in hydrostatic equilibrium in an abritrary external po-tential.Hydrostatic equilibrium implies a balance betweengravity and pressure. Gravity includes the self-gravity ofthe disk plus that of external components (i.e. dark mat-ter halo, bulge, stellar disk, etc), while the pressure is givenby an equation of state p = p ( ρ g , T ), with p being the gaspressure, ρ g the gas density and T the temperature. Thechallenge is to find a ρ g , T and the velocity field, ~v , suchthat the system is self-consistent (i.e., obeys the Poissonequation) and in hydrostatic equilibrium.In the case of an isothermal, axisymmetric, perfectlyself-gravitating disk (i.e., no external potential), the equi-librium disk has a sech distribution (Spitzer 1942) in thevertical direction, with a scale-height that is proportional to p c /ρ g ( R, z = 0), where c s is the sound speed. Here, cylin-drical coordinates, ( R, φ, z ), are used to describe the den-sity field. This immediately shows that since ρ g ( R, z = 0) istypically a decreasing function of radius, one generally ex-pects the scale-height to be a function of R . In particular,in the case of a globally isothermal disk, the sound speed c ∝ T is constant in space, giving rise to a flaring disk, i.e.,the scale-height increases with increasing R (Narayan & Jog2002, hereafter NJ02; Agertz et al. 2009). Alternatively, ifwe want to initialize a disk with a constant scale-height, aradial temperature gradient needs to be introduced. Tasker& Bryan (2006) initialize their disks to be isothermal and tohave a constant scale-height. As indicated above, this can-not be an equilibrium configuration. Consequently, the diskis expected to experience an unavoidable relaxation processwhich makes the initialization not well-controlled and mightpotentially contaminate the physics, e.g., star formation, gasdynamics etc., of interest. Agertz et al. (2009) set-up theirisothermal disks based on the local total surface density ofgas plus dark matter. Although the scale-height of their ini-tial disk changes with radius, the local total surface densityis not defined in a mathematical way and therefore elusive.In addition, their surface density does not follow an expo-nential profile.An important assumption underlying Spitzer’s analysisis that the radial variation in the potential is negligible com-pared to that in the vertical direction. This assumption issupported by observation that disks typically have verticalscale-heights that are an order of magnitude smaller thantheir radial scale-length (van der Kruit & Searle 1981a,b).A well studied example is the Milky Way, whose scale-heightis less than 200 pc for the cold gas (Jackson & Kellman 1974;Lockman 1984; Sanders et al. 1984; Wouterloot et al. 1990;see also Narayan & Jog 2002) and roughly 300 pc for thestars in the Solar neighborhood (Binney & Tremaine 2008,Kent, Dame & Fazio 1991), compared to a radial scale-lengthof ∼ . × ×
64) on a distorted grid structure in the radial,the azimuthal and the vertical directions. Unlike the live par-ticles, these particles are simply used as markers for massdistribution. Second, they compute the joint total potentialand the resulting force field numerically with a hierarchicalmultipole expansion based on a tree code. Third, given thepotential just computed, integrating Eq. (2) for a given mid-plane volume density, ρ g ( R, z = 0), gives the vertical struc-ture of density. Fourth, adjust the midplane volume densityto fulfill Eq. (24). Repeat the procedure between the secondstep and the fourth until the result converges.Although this description is quite general and flexible,for several reasons, this is not commonly used in the grid-based codes which are featured with adaptive-mesh refine-ment (AMR). The first and also the most fundamental oneis that the grid structure is normally unknown before we ac-tually initialize the disk. Except the uniform-grid initializa-tion, the grid structure is automatically generated based onthe criterion for refinement. Second, for a fully parallelizedcode, the initial data is distributed over different processorsand memory storages. This means that the data exchangebetween processors is necessary in order to fully computethe joint total potential. The situation becomes even moretechnically challenging when initializing with AMR. Third,The vertical structure of the gas disk depends only on thevertical potential difference (see Eq. (7) and Eq. (9) below).A description of the equatorial potential is enough for speci-fying the velocity field (See Eq. (13), Eq. (29) and the resultsshown in Sec. 3). In general, given the density distributioncomputed by the methods introduced in Section 2.2 togetherwith the conclusion in Section 2.1, we are allowed to acquirethe exact velocity field by Eq. (A.17) in Casertano (1983).Fully solving the Poisson equation becomes not necessary.Fourth, initializing a disk over distributed memories allowus to deal with a larger data set which cannot be fully con-tained in a single memory storage.We propose a simple but very effective way of initializ-ing a three-dimensional gas disk. This method can be easilyincorporated into any existing code based on either a La-grangian or Eulerian approach. No data exchange betweenprocessors is needed. Vertical density profile is obtained self-consistently without solving the full Possion equation. Weimplement these ideas with the adaptive mesh refinementmagnetohydrodynamics code RAMSES (Teyssier 2002) andapply our concepts to probe the onset of the disk instability.We modify the dispersion relation for the infinitesimally thindisk (Lin & Shu 1964) to be able to treat thick disks (Gol-dreich & Lynden-Bell 1965; Kim & Ostriker 2002a, 2006;Shetty & Ostriker 2006; Lisker & Fuchs 2009). The thresh-old value Q th is then obtained semi-analytically. Previousstudies on this subject are either focused on a small patchof a galaxy (2D/3D: Kim & Ostriker 2002a) or are globallytwo-dimensional but with the reduction of gravity includedin the governing equations (Shetty & Ostriker 2006). In thispaper, we revisit the subject as a test of our fully three- c (cid:13) , 1– ?? quilibrium Initialization and Stability of Three-Dimensional Gas Disks dimensional isolated galaxy models. Models with or withoutstellar potential are investigated.Galactic disks are comprised of stars and gas. Both com-ponents are coupled to each other via the Poisson equation.Since the stellar disk dominates the mass budget within theluminous disk, its presence has great impact on the scale-height of the gas disk as described in NJ02. A balancedinitial condition depends not only on the correct verticalstructure but also on the correct rotation velocity. To spec-ify the rotation velocity needed, the mass enclosed within acertain radius must be under control. Although it is com-mon practice to specify the functional form of the volumedensities of 3D disks, we show that because of the flaringdisk this typically results in a surface density profile thatcontains a central ‘hole’ (Agertz et al. 2009). This problemcan be trivially avoided by specifying the desired surfacedensity profile instead. We show that the corresponding vol-ume density can easily be obtained using a simple iterativescheme. The surface density of the total gas (HI + H ) fromobservation (Leroy et al. 2008) typically follows an expo-tential profile in disk galaxies. This profile gives an analyticdescription of the total mass enclosed within a radius as wellas a reasonable approximation for velocity field as shown byEq. (29) below (Binney & Tremaine 2008).Describing the stellar disk with a fixed background po-tential is at best an approximation to reality. The interactionbetween live stellar disk and gas can potentially destabilizethe system (Rafikov 2001; Li, Mac Low & Klessen 2005a,2005b, 2006; Kim & Ostriker 2007). After all, the gas is coldcompared to the stellar disk and has highly non-linear re-sponse to the asymmetric stellar potential. The gravitationalinterplay between the collisionless stars and dissipative gasis important for a number of key questions in galactic dy-namics. For example, what is the physical origin of grand de-sign spirals? Or what initiates and regulates the formationof stars? Having access to well-controlled initial and envi-ronmental conditions is a prerequisite to discovering theircauses.This paper is organized as follows. The ideas of initial-izing a gas disk are outlined in Section 2. Details of the sim-ulation parameters and test runs are described in Section 3.Axisymmetric instability of the disk is revisited in Section4. The self-induced spirals due to swing amplification willbe discussed in Section 5. A brief summary and the possibleextension of this work is put in Section 6. In this Section, we develop the required relations and equa-tions to immerse a 3D gas disk in a preexisting static poten-tial. Assuming that the gas disk and the preexisting poten-tial share the same symmetry axis, cylindrical coordinates,(
R, φ, z ), are adopted to formulate the dynamics of the sys-tem. Axial-symmetry enables us to discard the terms de-scribing the variation in azimuthal direction, i.e., ∂/∂φ = 0.A gas disk which is in detailed balance should fulfill thefollowing equations:1 ρ g ∂p∂R + ∂ Φ ∂R = V R , (1)1 ρ g ∂p∂z + ∂ Φ ∂z = 0 , (2) with ρ g , p , V rot and Φ being the volume density of the gas,the gas pressure, the azimuthal rotation velocity (“rotationvelocity” in short) and the joint total potential. Equation(1) describes that the gravitational pull in radial directionis counterbalanced by the centrifugal force and the pressuregradient. Equation (2) states that hydrostatic equilibriumalong the symmetry axis, the z -direction, is achieved by thebalance between vertical pull of the gravity and the pressuregradient in z .To make the system self-consistent, the Poisson equa-tion must be involved: ∇ Φ = 4 πG ( ρ g + ρ DM + ρ s ) , (3)with G , ρ DM and ρ s being the gravitational constant, and thevolume density of dark matter and stars. The total potentialis comprised of the contributions from the dark matter halo,the stellar disk and the self-gravity of the gas, i.e., Φ =Φ DM + Φ s + Φ g . In addition, the ideal gas law provides thelink between the gas density, the gas temperature and thegas pressure: p = ρ g ( γ − e ( T ) , (4)where γ represents the ratio of the heat capacities (adiabaticindex), e the specific internal energy and T the gas temper-ature. In the case of an ideal gas, the specific internal energydepends only on temperature, and is given by e = 1 γ − k B Tµm p , (5)with k B being the Boltzmann’s constant, µ the atomicweight and m p the mass of a proton. However, to close theset of equations, we should either invoke the energy equationor an equation of state (EoS), which will be used to evolvethe system.A disk which is in hydrodynamic equilibrium shouldstay in its original state if we evolve the disk with the sameequation of state which is used to set-up the disk. Thenumerical results throughout this paper are based on theisothermal equation of state, i.e., p = c ρ g , (6)with c s being the sound speed, a temporal and spatial con-stant. Equations (1) to (6) then form the basis of our discus-sion. In this paper, all the disks are in detailed equilibriumwith the isothermal EoS. If those disks are adopted to evolvewith a cooling function or a polytropic EoS, we can makesure any change in temperature or dynamics is purely causedby a cooling or a heating source.For a polytropic gas, p = Kρ Γg . with Γ and K beingconstant, integrating Eq. (2) gives: ρ g ( R, z ) = ρ g ( R, z = 0) (cid:20) − Γ − c ( R, z = 0) Φ z ( R, z ) (cid:21) − , (7)where Φ z ( R, z ) = Φ(
R, z ) − Φ( R, z = 0) defines the verti-cal potential difference. We have used the fact that c ≡ ∂p/∂ρ g = K Γ ρ Γ − when approaching Eq. (7). Note thatgiven Γ = 1, the internal energy has the following relation: e ( T ) = Kρ Γ − γ − . (8)Combining Eq. (5) and Eq. (8) gives the temperature field c (cid:13) , 1– ?? Wang et al. as a function of position if the gas disk is initialized with anon-isothermal EoS. As a special case, when Γ →
1, Eq. (7)then converges to a form for the isothermal gas: ρ g ( R, z ) = ρ g ( R, z = 0) exp (cid:18) − Φ z ( R, z ) c (cid:19) . (9)As we can see from Eq. (7) and Eq. (9), the vertical structureof gas disk depends on the gas properties in the midplaneand the vertical potential difference.To fully characterize a gas disk which is in detailed bal-ance, we need to specify the velocity, the density and thetemperature at every location in the beginning of the sim-ulation. In the following sub-sections we study the generalproperties of the velocity and density field, which allows usto devise a simple, but effective method to initialize a 3Dgas disk in hydrostatic equilibrium. In this sub-section, we treat the azimuthal rotation velocityas generally as possible. To make the notation concise, wedrop the subscript of gas density, ρ g , and restore the sub-script after this sub-section. Without further assumption,integrating Eq. (2) from 0 to z gives: Z z ρ ∂p∂z d z = − Φ z ( R, z ) . (10)By integrating Eq. (10) in parts, we have: p ( R, z ) ρ ( R, z ) = p ( R ) ρ ( R ) (cid:12)(cid:12)(cid:12)(cid:12) z =0 − Z z pρ ∂ρ∂z d z − Φ z ( R, z ) . (11)Inserting Eq.(11) into (1), we get (see Appendix A for fur-ther details): V ( R, z ) R = 1 ρ ∂p∂R + ∂ Φ( R, z ) ∂R = 1 ρ ∂p∂R (cid:12)(cid:12)(cid:12)(cid:12) z =0 + ∂ Φ( R ) ∂R (cid:12)(cid:12)(cid:12)(cid:12) z =0 − Z z (cid:26)(cid:16) ∂ρ∂z (cid:17) ∂∂R (cid:18) pρ (cid:19) − (cid:16) ∂ρ∂R (cid:17) ∂∂z (cid:18) pρ (cid:19)(cid:27) d z. (12)Equation (12) shows that the rotation velocity is indepen-dent of height above or below the midplane so long as theintegral vanishes. It is evident that gas with a barotropicequation of state, i.e., p ( ρ g ), fulfills this requirement. In ad-dition, for an initially constant temperature ( T is everywherethe same in the beginning) disk, the initial pressure is afunction of volume density only and therefore the integralbecomes zero. In this case, equation (12) can be simplifiedfurther: V ( R, z ) =
R ∂ Φ ∂R (cid:12)(cid:12)(cid:12) z =0 + ( γ − e ∂ ln ρ∂ ln R (cid:12)(cid:12)(cid:12) z =0 . (13)Equation (13) states that the process of specifying the initialvelocity in 3D comes down to knowing the rotation velocityin the equatorial midplane. From now on, to avoid confusion, we restore the subscript forthe gas density. To proceed further, we consider the gas layer to be a very thin structure embedded in a static potentialcontributed by the background spherical dark matter andthe stellar disk. Because the gas disk is observationally thinwe neglect the radial variation compared to the vertical one(i.e., | ( ∂/∂R ( R∂ Φ g /∂R )) /R | ≪ | ∂ Φ g /∂z | ). In Appendix Ewe show that this is a valid assumption for realistic gas disks.For an axisymmetric thin disk, the Poisson equation thenreduces to (Binney & Tremaine 2008):d Φ g d z = 4 πGρ g . (14)with Φ g being the potentials contributed by the gas. In thefollowing, we focus only on disks with initially constant tem-perature, i.e., the rotation velocity required for equilibriumhas no dependence on the height above or below the mid-plane. Differentiation Eq. (2) with respect to z and insertingEq. (14) leads to a second-order non-linear differential equa-tion:d p d z − ρ g d ρ g d z d p d z + ρ g (4 πGρ g + d Φ s d z + d Φ DM d z ) = 0 , (15)with Φ s and Φ DM being the potentials contributed by thestellar disk and the dark matter, respectively. So far, Eq.(15) is still general with respect to any kind of equation ofstate. However, a single equation with two unknowns p and ρ g is not solvable. To continue with Eq. (15), in this sub-section, we assume that the gas is barotropic, i.e., p ( ρ g ).Given density distributions of stars, the dark matter andthe boundary conditions in the midplane: ρ g ( R,
0) = ρ ( R ) and d ρ g d z = 0 , (16)equation (15) can be solved by numerical integration, e.g.,using the Runge-Kutta method. For a single-component,self-gravitating, locally isothermal disk ( c s ( R ) can be a func-tion of radius), Eq. (15) has an exact solution: ρ g ( R, z ) = ρ ( R )sech ( z/h ) , (17)with ρ ( R ) being the gas volume density in the midplane, h = p c / πGρ the scale-height and c s the local isothermalsound speed. According to Eq. (17) and since the midplanevolume density, ρ ( R ), generally decreases with radius, tokeep the scale-height a constant, the sound speed and there-fore the temperature must be a function of radius.Equation (15) is the simplified version of the formuladerived by NJ02 (see also Kim et al. 2002a), where they in-vestigated the vertical structure in a gravitationally coupled,multi-component galactic disk. It is important to noticethat all calculations can be done locally without the needof exchanging information between processors and thereforegreatly reduces the complexity of coding.In principle Eq. (15) allows one to compute the den-sity of the gas such that the disk initially is in hydrostaticequilibrium. The actual implementation using Eq. (15) doesnot guarantee the positivity of the density. In particular,at large radii ρ g ( R, z ) is typically close to zero, and smallerrors due to the numerical integration often yield negativedensities. This problem is especially relevant when using theadaptive-mesh refinement techniques. c (cid:13) , 1– ?? quilibrium Initialization and Stability of Three-Dimensional Gas Disks Initializing a gas disk with AMR usually starts with thecoarsest grid. A natural selection of the integration step isthe cell size. Then a problem immediately rises when solvingEq. (15) to specify the volume density. Supposing that thecell size is much larger than the scale-height of the gas disk,the errors introduced by the coarse integration may lead tonegative densities on the outskirt of the computation do-main. One might think the integration can be done by usingeither adaptive integration intervals or simply a fixed inte-gration interval which is much smaller than the cell size.However, the improvements only guarantee the convergenceof the solution not the positivity. Nevertheless, because ofthe generality of Eq. (15), density method is still valuablefor semi-analytic research.
In this sub-section, we develop another route for specifyingthe density distribution. We stress that the following deriva-tion is only applicable to initially isothermal disks. With thisconstraint, integrating Eq. (2) gives: ρ g ( R, z ) = ρ ( R ) exp (cid:18) − Φ z ( R, z )( γ − e (cid:19) . (18)Combining Eq. (14) and Eq. (18), a second-order non-linearequation for the vertical potential difference of gas is ob-tained:d Φ g ,z d z = 4 πGρ ( R ) exp (cid:18) − Φ z ( R, z )( γ − e (cid:19) . (19)Given the analytic forms of Φ DM and Φ s the only unknown isthe potential difference of gas, Φ g ,z = Φ g ( R, z ) − Φ g ( R, z =0). Similar to the density method, given the boundary condi-tions ρ ( R ), Φ( R, z = 0) and dΦ(
R, z = 0) / d z = 0, numeri-cal integration can be applied to solve Eq. (19). By insertingΦ z obtained by integrating Eq. (19) into Eq. (18), the den-sity distribution is acquired. Notice that what really mattersto us is the potential difference, not the absolute value. Thismeans the value of Φ g ( R, z = 0) can be an arbitrary con-stant, although we do know the values of Φ DM ( R, z = 0)and Φ s ( R, z = 0).The merit of this formulation is evident, the occurrenceof negative density is avoided by Eq. (18). Tiny errors in thepotential difference will not do any harm to the positivityof the gas density. Numerical experiments show that in nor-mal cases in which both the density method and potentialmethod work, the solutions are consistent.At a given radius, R , solving Eq. (19) only provides uswith information about the potential difference, Φ z ( R, z ).However, a useful byproduct of the potential method is thatit is possible to acquire a good approximation of the totalpotential by:Φ g ( R, z ) = Φ g ( R, z = 0) + Φ g ,z ( R, z ) , (20)as long as we know the potential in the midplane, Φ g ( R, z =0). Equation (20) is an approximation because the use ofEq. (19) is based on the reduced Poisson equation Eq. (14)in which the variation in radial direction is ignored. Thegradient of the potential Φ g ( R, z = 0) determines the ve-locity field required while the vertical potential differenceΦ g ,z ( R, z ) gives the vertical structure of the disk. In princi-ple, the radial force, which is associated with Φ g ( R, z = 0), in the equatorial plane for an axially symmetric density dis-tribution can be evaluated precisely by the equation (A.17)in Casertano (1983). This allows us to obtain the total po-tential without fully solving the Poisson equation. In prac-tice, if the initialization is performed with multi-node clus-ters, each node only keeps part of the information about thedensity distribution, data exchange with AMR itself is tech-nically challenging. In Section 3, for an exponential disk, thenumerical results will show that the use of Eq. (29) is a goodapproximation for most of our interests. The correspondingΦ g ( R, z = 0) associated with Eq. (29) can be found in thebook by Binney & Tremaine (2008), Eq. (1.164a).Equation (20) is useful, because involving the total po-tential into the formulation is an important step for self-consistently building up the combined disks comprised ofa live stellar disk and a gas disk. Extension to the workof Shu (1969), Kuijken & Dubinski (1995, hereafter, KD95)develop a self-consistent disk-bulge-halo model for galaxies.The distribution function built by Eq. (6) in KD95 involvesthe potential differences Φ z and Φ( R, − Φ( R c , R c the radius of the guiding center. The potential method pre-sented here can be naturally incorporated into the frame-work of KD95. Therefore, in this paper, all the disks areinitialized by the potential method. Some studies have assumed that the midplane density of a3D gas disk has an exponential form (Tasker 2006, Agertz2009). However, as we now demonstrate, in general this re-sults in a surface density distribution that peaks at a specificnon-zero radius, giving rise to a ring-like feature. We assumea gas disk with the popular sech vertical profile: ρ g ( R, z ) = ρ c exp( − R/R d )sech (cid:18) zh ( R ) (cid:19) , (21)with ρ c being the central volume density, R d the disk scale-length and h ( R ) the scale-height as a function of radius. Thesurface density then reads:Σ( R ) = Z ∞−∞ ρ g ( R, z )d z = 2 ρ c exp( − R/R d ) h ( R ) . (22)Based on Eq. (22), we measure the scale-height of a diskat certain radius by h ( R ) = Σ( R ) / (2 ρ ( R )). The extrema ofthe surface density can be evaluated by taking the derivativeto Eq. (22):dΣ( R )d R = 2 ρ c exp( − R/R d ) (cid:18) d h ( R )d R − h ( R ) R d (cid:19) = 0 . (23)We Suppose that the disk is linearly flaring, i.e., h ( R ) = h + R/R h , with h being the minimum scale-height of the diskand R h a factor controlling the degree of flaring. The peakof the surface density then locates at R peak = R d − h R h .Whenever the R peak is positive, we get a ring in surfacedensity. However, a ring in the surface gas density is notcommonly seen in a real disk galaxy. An exponential pro-file in the total gas is prevalent in disk galaxies (Leroy etal. 2008).In order to avoid this feature, it is advantageous to spec-ify the actual surface density of the disk, rather than its c (cid:13) , 1– ?? Wang et al.
Table 1.
Models’ Parameters.Run T (K) M s (M ⊙ ) FigureGas0 4 × - (1),(2),(3)Gas1 2 × - (4),(5)Gas2 1 × - (4),(5),(8)Gas3 9 × - (4),(5)Gas4 8 × - (4),(5)GasStar1 7 × × (6),(7),(8)GasStar2 6 × × (6),(7)GasStar3 5 × × (6),(7)GasStar4 4 × × (6),(7)*All disks have a gas mass of 10 M ⊙ . midplane density. In the case of the exponential profile, thesurface density reads:Σ( R ) = Σ exp ( − R/R d ) = Z ∞−∞ ρ g ( R, z )d z, (24)with Σ being the surface density in the galactic centre.Combining Eq. (24) and Eq. (18), the volume density in themidplane can be expressed as: ρ ( R ) = Σ exp ( − R/R d ) R ∞−∞ exp ( − Φ z / [( γ − e ])d z . (25)It shows that the correct volume density in the midplanefor the desired surface density profile can be obtained iter-atively. Given a initial guess for ρ ( R ), Φ z is evaluated viaEq. (19) and also the integral appears in Eq. (25). One needsto iterate between Eq. (19) and Eq. (25). However, depend-ing on the quality of the initial guess, convergence can bereached very fast. For instance, with the initial guess being ρ ( R ) = Σ exp( − R/R d ), a six-time iteration already givesus a reasonable exponential disk.We pursue the exponential disk for several reasons. Oneis simply because it is commonly seen in disk galaxies. An-other is that we have a better control of the total mass. Aswe can see, if we specify the midplane volume density insteadof the surface density, we do not exactly know the total massuntil we finish the integration. Without knowing the totalmass in advance, evaluating the circular velocity contributedby the self-gravity will not be a trivial task. Nevertheless, inprinciple, any profile of the surface density can be achievedsimply by the process introduced in this sub-section. In this Section, we test the ideas outlined in the previous Sec-tion. We implement the method in the AMR-code RAMSES(Teyssier 2002). RAMSES uses grid-based Riemann-solversfor the magneto-hydrodynamics (MHD) and particle-mesh(PM) technique for the collisionless physics. It has a fullyparallelized Poisson solver with periodic boundary condi-tions, which we use for this paper. Gas disks which are ini-tialized isothermally with an exponential surface density ofa scale-length of 3.5 kpc and a total mass of 10 M ⊙ areembedded in a static potential. An isothermal equation ofstate is used to evolve the disks throughout this paper.The tests are mainly divided into two groups, one group is evolved with a static stellar potential (models with theprefix GasStar), the other without (models with the prefixGas). Gas1 to Gas4 are M33-like gas-rich galaxies, whileGasStar1 to GasStar4 are more similar to the Milky-Way.The main parameters of the models are listed in Table 1.The size of the computational domain is 250 kpc on a side.Up to 12 levels of refinement are used for those runs withoutstellar potential, and 13 levels for the other group, i.e., thecorresponding highest spatial resolutions are about 60 pcand 30 pc, respectively.The volume density of the halo is described by the NFWprofile (Navarro, Frenk & White 1997): ρ DM ( r ) = M πf ( c ) r cxr (1 + x ) , (26)with the Virial mass M = 10 M ⊙ , x = rc/r , concen-tration parameter c = 12, distance r = √ R + z , Virialradius r = 213 kpc and f ( c ) = ln(1 + c ) − c/ (1 + c ). TheVirial radius ( r ) is a radius within which the averagedmatter density is 200 times the critical density.The density distribution of the stellar disk reads(Miyamoto & Nagai 1975, Binney & Tremaine 2008): ρ s ( r ) = (cid:18) b M s π (cid:19) aR + ( a + 3 √ z + b )( a + √ z + b ) (cid:2) R + ( a + √ z + b ) (cid:3) / ( z + b ) / , (27)with M s = 4 × M ⊙ being the mass of the stellar disk, a = 3 . b = 0 . V rot , is decomposed into fourcomponents: V = V + V + V + V , (28)where V DM , V s , V g are the circular velocities correspondingto the dark matter halo, the stellar disk and the gas disk, and V p denotes the contribution due to the pressure gradient.In this paper, we have the analytic form for V DM and V s .For the contribution from the gas disk and pressure gradient,we take the approximation for an infinitesimally thin diskwith exponential surface density as described in Eq. (24).We set: V ( R ) = 4 πG Σ R d y [ I ( y ) K ( y ) − I ( y ) K ( y )] (29) V ( R ) = ( γ − e ∂ ln ρ∂ ln R (cid:12)(cid:12)(cid:12) z =0 , (30)with y = R/ (2 R d ), I , K , I and K being the modifiedBessel functions of the first and second kinds of zeroth/first-order, respectively. Equation (30) derives from the secondterm of Eq. (13). However, contribution from pressure gra-dient in the midplane can only be evaluated after the gasdisk is set up. Note that for an exponential disk, surfaceand volume densities decrease with radius and hence V isnegative. To demonstrate that the disk built by the potential methoddescribed in Section 2 is in detailed equilibrium, we startwith a stable equilibrium disk in model Gas0. In this test, thestellar disk is deliberately removed. Without the dynamical c (cid:13) , 1– ?? quilibrium Initialization and Stability of Three-Dimensional Gas Disks V r o t [ k m / s ] halogaspressuretotal Q R [Kpc] R [Kpc] (a) (b)
Figure 1.
Model Gas0: (a) The total rotation velocity (solid) and contributions from dark matter halo (dashed), gas (dotted) and gaspressure (dash-dotted). Note that we plot the absolute value of the pressure gradient to have positive values for the direct comparison. Itshould be in opposite sense to the gravity. In this model contributions from the gas self-gravity and the pressure gradient is not negligible.(b) The Q value of model Gas0 as a function of radius as defined by Eq. (32). The Q is well above the threshold value Q th = 1, thus thedisk is expected to be stable. No structure should develop with time. R [Kpc] Σ [ M O • / kp c ] R [Kpc] V r o t [ k m / s ] R [Kpc] s ca l e - h e i gh t [ K p c ] t=0 Gyrt=1.6 Gyr (a) (b) Figure 2.
Model Gas0: The evolution of 1.6 Gyrs of a stable disk. (a) The evolution of the surface density (upper panel) and therotation curve (lower panel) at t = 0 Gyr (solid) and t = 1 . t = 0 Gyr (solid) and t = 1 . support from the stellar disk, the self-gravity of the gas playsthe dominant role in determining the vertical structure ofthe disk and provides a not negligible contribution to therotation velocity.Figure 1a decomposes the rotation curve into the dif-ferent contributions by the halo, the gas and the pressuregradient. Note that the forces of the self-gravity and thepressure gradient are in opposite sense, the self-gravity pullsmatter inwards while the pressure gradient pushes outwards.In this figure, V p is shown in its absolute value. If we ignorethe pressure gradient, the disk would rotate too fast andgradually drift outward. Figure 1b shows the conventionalToomre’s Q defined by: Q = c s κπG Σ g . (31)with κ being the epicyclic frequency. It shows that the Q iswell above Q th = 1, the threshold value for stability, at all radii. The disk is hot enough to keep the disk stable and nostructure should develop with time.We let the disk evolve for 1.6 Gyrs (four orbits for thegas at 10 kpc) and check how well the disk properties arekept. Figure 2a presents the evolution of the surface den-sity and the rotation curve. The solid lines represent theinitial states and the diamond symbols are the status afteran evolution of 1.6 Gyrs. The surface density is obtainedby projecting along the symmetry axis and the rotationcurve is evaluated by the mass-weighted circular velocity,¯ v rot ( R ) = R ρ g ( R, z ) v rot ( R, z )d z/ Σ g ( R ). Although a smallamount of mass accretes onto the very central part of thedisk, overall the surface density and the rotation curve arekept very well. Mass accretion into the center seems unavoid-able for a Cartesian-grid code mainly because too small anumber of cubic cells is used to mimic the circular motion c (cid:13) , 1– ?? Wang et al. (a) (b)
Figure 3.
Model Gas0: The size of the images is 20 kpc ×
20 kpc. The evolution of the surface density at (a) t = 0 Gyr and (b) at t = 1 . in the centre. This accretion will be eventually halted by thepressure gradient built by the accumulating material.Figure 2b shows the evolution of the scale-height. Thesolid line represents the initial state and the dotted linethe evolution after 1.6 Gyrs. Upon closer inspection wefind that the disk undergoes a very small amount of mass-redistribution in the radial direction, which we believe tobe a consequence of our two approximations when initializ-ing the disk. One is from the reduced Poisson equation, Eq.(14), and the other is from the use of Eq. (29). Equation(29) overestimates the circular velocity needed to supportthe disk. The thickness of the disk reduces the potential inthe midplane by a few percent (see appendix B). Figure 3shows the snapshots of the face-on surface density at t = 0(Fig. 3a) and at t = 1 . The question of disk stability has been investigated for morethan four decades since the pioneering works by Toomre(1964) for collisionless stars and Goldreich & Lynden-Bell(1965) for gas sheets. Understanding the origin and evolu-tion of disk structure is challenging. If the disk is stablelike our model Gas0, no structures can form. On the otherhand, if the disk is highly unstable, the surface density willquickly fragment and develop a clumpy and chaotic-lookingappearance. There will be no well-organized structures. Thestriking spiral appearance of many nearby disk galaxies in-dicates that those disks are marginally stable.For an infinitesimally thin disk, the instability thresh-old is at Q th = 1 (Toomre 1964). The first theoretical workto include the finite thickness of a self-gravitating gas disk is that by Goldreich & Lynden-Bell (1965). Some authors haveinvestigated the stability of finite thickness gas disks in nu-merical simulations (both in 2D and 3D) using local patcheswithin a shearing box (Kim & Ostriker 2006, 2002a; Gam-mie 2001). This technique, in 2D, has also been used by Kim& Ostriker (2007) to investigate the interaction between thegas disk and a live stellar disk. Shetty & Ostriker (2006) usedglobal 2D simulations in which they incorporated the effectof finite disk thickness by diluting the gravitational force.For 3D global disk calculations, see Li, Mac Low & Klessen(2005a, 2005b, 2006), who investigate the relation betweendisk instability and star formation rate. These studies allagree that although the inclusion of the thickness does nothave a qualitative impact on the disk instability, it does shiftthe threshold value of instability quantitatively. In addition,accounting for disk thickness may have a large impact onthe evolution of a disk, such as the development of spurs orthe wiggle instability (Kim & Ostriker 2002b, 2006).In this paper, armed with a well-balanced gas disk, werevisit the axisymmetric instability of disks in 3D globalfashion. We first derive the reduction factor F which re-flects the reduction of the gravity due to the finite thicknessof the disk. Then the corresponding instability threshold Q th ( R, T ) derived from a semi-analytic calculation is com-pared with the numerical results. In the final sub-section,we also explore the impact of the presence of a static stellarpotential on the axisymmetric instability.
The Fourier component of the perturbed gravitational po-tential, Φ k , of an infinitesimally thin disk is given by:Φ k = − πG Σ k | k | e ikx − k | z | , (32)where k represents the wave number of the Fourier compo-nents and x = R − R being the radial deviation for an ax-isymmetric perturbation. Supposing that a 3D disk is piledup by a stack of infinitesimally thin gas layers, we approx-imate the effect of the disk thickness by superimposing thecontribution from every razor-thin layer: c (cid:13) , 1– ?? quilibrium Initialization and Stability of Three-Dimensional Gas Disks Φ k ( z ) = − πG Σ k e ikx | k | Z ∞−∞ e − k | z − h | sech ( h/h z )2 h z d h, (33)with h z being the scale-height of the disk. In Eq. (33), wemodel the vertical structure of the gas disk by a sech func-tion. This is valid especially for the inner part of diskswhere the vertical structure is mainly determined by theself-gravity of the gas. See also the Fig. D1 in Appendix D.Equation (33) leads to the Fourier potential in the midplane:Φ k ( z = 0) = − πG Σ k e ikx | k | F ( k, h z ) , (34)with F ( k, h z ) being the reduction factor described by (seeAppendix C): F ( k, h z ) = 1 − kh z h H (cid:16) kh z (cid:17) − H (cid:16) kh z − (cid:17)i , (35)with H being the harmonic number defined by: H ( α ) = Z − y α − y d y. (36)The Lin-Shu (1964) dispersion relation for the axisymmetricperturbation is then modified to: ω = κ − πG Σ | k | F ( k, h z ) + c k . (37)The dispersion relation states that on small scales ( k → ∞ )the disk is stabilized by gas pressure, i.e., the term c k .Large scales ( k →
0) are regulated by global shear, i.e., the κ term. The instability however happens at intermediatewavelengths, much smaller than the disk size but still largerthan the thickness of the disk. In this region, neither globalshear nor gas pressure can resist the gravitational collapse.The reduction factor, 0 < F ≤
1, softens the effect of self-gravity and makes the disk more stable.Given a certain radius R and temperature T , we ob-tain the threshold value Q th ( T, R ) by probing the maximumvalue along the neutral curve defined by setting ω = 0 inEq. (37) and calculating the epicyclic frequency, κ , fromthe rotation curve. Similar to the conventional Toomre cri-terion for the stability of an infinitesimally thin disk, Q th is a threshold curve for thick disks. Above Q th the disk isstable and otherwise unstable. Since the Q th is a functionof both temperature and radius, it is convenient to definethe critical value Q crit , which is the value of Q th for which Q th ( T, R ) /Q ( R ) = 1, and the corresponding critical tem-perature T crit .The solid lines shown in Fig. 4 represent the thresholdvalue Q th as a function of radius. Each plot corresponds to adisk of different temperature. The dash-dotted lines are theactual Q values defined by Eq. (31) of the different models.From these figures, the most unstable radius is around R = 2kpc. The corresponding surface densities after an evolutionof 750 Myrs are shown in Fig. 5. The gas at the most unsta-ble region has revolved for more than four orbital periodsaround the disk center.These figures shows that the prediction of Q crit and thenumerical results match quite well. The Q value of Gas1is well above the solid line and shows a featureless surfacedensity. As shown in Gas2 and Gas3, with the decrease intemperature, the Q th curves shift up and the disks’ Q curvescome down. As a consequence, the disk starts to developmulti-armed structure, which is very likely caused by swing amplification, as discussed in Section 5. And finally in Gas4,the curves Q and Q th intersect. The disk fragments andstarts to behave chaotically. A more detail calculation showsthat the two curves just touch each other at a temperature T crit = 8 . × K with the maximum threshold Q crit =0 . Q crit = 0 .
676 of Goldreich & Lynden-Bell’s (1965) analysis but away from the numerical result, Q crit = 0 . Q crit is model dependent. Different models of thedark matter, the stellar disk and even the EoS will all affectthe resulting value of Q crit . The inclusion of a static stellar disk alters two importantfactors which influence the stability of the disk. One is therotation curve and the other is the thickness of the gasdisk. By changing the rotation curve, the epicyclic frequency, κ , changes accordingly. Supposing a flat rotation curve de-scribed by Ω = V /R , the epicyclic frequency κ then reads: κ = 2Ω = 2 V R , (38)with V being the rotation velocity. The presence of a stel-lar disk tends to stabilize the gas disk via increasing V .However, by increasing the gravitational pull in the verticaldirection, the gas disk becomes thinner and therefore moresusceptible to gravitational collapse. In Section 4.1, we havealready seen that the scale-height, which is governed by thetemperature of the disk, is a very sensitive factor for thedisk stability. GasStar1 to GasStar4 are designed to explorethe competition between the two opponents.From Fig. 6, we first notice that, compared to Fig. 4, thethreshold value, Q crit , is boosted from 0.693 to 0.75 due tothe decrease in scale-height. This makes the disk more proneto gravitational instability. On the other hand, the changeof the rotation curve drastically shifts the dash-dotted curveupwards. Instability only sets in once the temperature of thegas disk drops below T crit ∼ Q (dash-dotted line)is approaching the threshold Q th (the solid line). The spi-rals are sheared, become tighter and tighter and enhanced.Once the density reaches the supercritical point, instabilitysets in. c (cid:13) , 1– ?? Wang et al. ( c) Gas3 (d) Gas4 R [Kpc] Q Q th Q 1 2 3 40.20.40.60.811.21.4
R [Kpc] Q Q th Q1 2 3 40.20.40.60.811.21.4
R [Kpc] Q Q th Q 1 2 3 40.20.40.60.811.21.4
R [kpc] Q Q th Q (a) Gas1 (b) Gas2 Figure 4.
Plots (a) to (d) correspond to models from Gas1 to Gas4, respectively. In each plot, curves of the disk Q (dash-dotted) andthe threshold value Q th (solid) are put together to probe the onset of axisymmetric instability. Q th ( R ) is a obtained by probing themaximum value along the neutral curve for a given radius. Information of the disk thickness has been encapsulated in the reductionfactor defined by Eq. (36). When the two curves meet, we expect the disk fragments very fast. This figure shows that the most unstableregion is about the radius R = 2 kpc. The fact that the Q th curves are well below unity shows the impact of the disk thickness on thedisk stability. An interesting feature which is hard to ignore in Fig. 5and Fig. 7 is that the marginally stable disks are sponta-neously developing multi-arm spiral structures. We have al-ready seen in Section 3 and 4 that the effect of the diskthickness is to shift the range of the marginally stable re-gion downwards and therefore to stabilize the disk. As wesystematically lower the temperature to probe the onset ofinstability, runs with as well as without stellar potential areexperiencing swing amplification.Hohl (1971) found that disks which are marginally sta-ble to axisymmetric perturbation are prone to develop alarge-scale bar structure. This finding initiated both numer-ical (Zang & Hohl 1978; Sellwood 1981, 1985; Fuchs & vonLinden 1998; Sellwood & Moore 1999) and theoretical stud-ies (Kalnajs 1978; Sawamura 1988; Vauterin & Dejonghe1996; Pichon & Cannon 1997; Evans & Read 1998; Fuchs2001) of marginally stable disks. Goldreich & Lynden-Bell(1965) and Toomre (1981) pointed out that self-gravitating,differentially rotating disks are able to amplify spiral wavesby shearing a leading wave into a trailing one. Three key in- gredients, self-gravity, shearing and epicyclic motions workharmonically to make the phenomenon now coined with thename ‘swing amplification’ happen.Three necessary conditions need to be fulfilled in orderto facilitate the swing amplification (Toomre 1981; Fuchs2001; Fuchs & von Linden 1998; Binney & Tremaine 2008).First, the disk must be marginally stable, i.e., for an in-finitesimally thin disk, 1 < Q <
2, as defined by Eq. (31).Second, the parameter X = k crit R/m = k crit /k y (Toomre1981; Binney & Tremaine 2008), with m being the num-ber of arms and k crit = κ / (2 πG Σ g ) the critical wave num-ber, has to be of order unity, i.e., somewhere between 1 and3 (Goldreich & Lynden-Bell 1965; Julian & Toomre 1966;Toomre 1981). Third, there must be a mechanism that isable to induce leading arms in the system either explicitlyby hand (Toomre 1981) or implicitly by random fluctuationinduced by numerical noise (Toomre 1990; Sellwood & Carl-berg 1984; Fuchs 2001). We notice that most of these worksmentioned above are for live stellar disks not directly for thegas disk. But since the amplification principles are the same,the results are still applicable to pure gas disks.As shown in Fig. 8a and 8b, GasStar1 gets more armsthan Gas2 does. To be more quantitatively, Fig. 8c and c (cid:13) , 1– ?? quilibrium Initialization and Stability of Three-Dimensional Gas Disks (c) Gas3 (d) Gas4(a) Ga s1 (b) Gas2 Figure 5.
Images (a) to (d) correspond to models Gas1 to Gas4, respectively. They show the face-on surface density at t = 750 Myr.The size of the images are 20 kpc ×
20 kpc. The gas at the most unstable radius has orbited around the center for more than fourtimes. (a) Since the disk Q is well above the threshold value Q th , the disk is featureless. In models Gas2(b) and Gas3(c) the disk Q isapproaching Q th around R = 2 kpc, both disks are developing self-induced spirals due to swing amplification. (d) The disk fragmentsvery fast once Q and Q th intersect.
8d show the Fourier components as a function of radius.They are obtained by doing Fourier transform to (Σ g ( R, φ ) − Σ g ( R )) / Σ g ( R ), where Σ g ( R ) is the averaged surface densityof a given radius. Note that the dominating modes tendsto be multiples of m = 4. This is a consequence of usinga Cartesian grid, for which m = 4 is the natural mode.However, the dominating mode is determined by physics.The dominating mode of Gas2 is m = 8 while in GasStar1 m = 12.As is apparent from Eq. (38), including a stellar diskcauses an increase in k crit . Consequently, a larger value of m is required in order to satisfy 1 < X <
3. From the imageshown in Fig. 8a and the relation, k crit ∝ κ , to keep X a constant, the number of spiral arms in GasStar1 can becrudely estimated as m ≃
15. More precisely, the number ofspiral arms, m , is predicted by (Toomre 1981; Athanassoula,Bosma & Papaioannou 1987; Fuchs 2001, 2008): m = 2 πRλ max , (39) with λ max being defined by: λ max = λ crit χ ( A/ Ω) , (40)where λ crit = 2 π/k crit . The coefficient χ is a function ofrotation curve (Fuchs 2001), as measured by Oort’s constant A . We employ Eq. (39) to analytically estimate the numberof arms and compare the predictions with the images shownin Fig. 8. For Gas2, spirals appear between 2 kpc and 5kpc. Within this radial range, the most unstable wavelengthranges from 2.0 to 3.6 kpc. The corresponding predictionfor m ranges from 6 to 9, while the simulation reveals aspiral pattern with 8-fold symmetry. For GasStar1, spiralsare prominent between 3 and 4 kpc, while the correspondingmost unstable wavelength ranges from 1.4 to 2.0 kpc. Thetwelve arms developing in GasStar1 should be compared to apredicted m ranging from 13 to 14. Hence, overall the trendsin the simulations are in reasonable agreement with our pre-dictions. Note that the spatial resolution in both simulations c (cid:13) , 1– ?? Wang et al. (c) Ga sStar
R [Kpc] Q Q th Q 1 2 3 40.20.40.60.811.21.4
R [Kpc] Q Q th Q1 2 3 40.20.40.60.811.21.4
R [kpc] Q Q th Q 1 2 3 40.20.40.60.811.21.4
R [kpc] Q Q th Q Figure 6.
Plots (a) to (d) correspond to models GasStar1 to GasStar4, respectively.: The Q (dash-dotted) and Q th (solid) curves ofthe gas disks of different temperatures. The presence the stellar potential stabilizes the disks through changing the rotation curve anddestabilizing the disk by increasing local gravitational force. The effect of disk thickness is included via the reduction factor Eq. (36).We need to lower the temperature down to T= 7 × K in order to probe the onset of axisymmetric instability. Overall, the presenceof the stellar potential stabilizes the disk. ranges from 60 pc to 120 pc, indicating that the most un-stable wavelengths are well-resolved.The observed small deviations can be explained as fol-lows. First, the formulation used to predict the numberof arms is precise only for stellar disks. However, Toomre(1981) has shown the strikingly similar behavior of gaseousdisks (Goldreich & Lynden-Bell 1965) and stellar disks (Ju-lian & Toomre 1966). Therefore, we have confidence thatEq. (39) is still applicable to gaseous disks. Second, thenumber of arms has to be an integer, a number of fractiongiven by Eq. (39) has no physical meaning. Third, the usageof a Cartesian grid introduces the multiples of the natural m = 4 mode, which manifests itself in the Fourier transformof the surface density. Fourth, swing amplification picks upthe dominating mode. It takes some time to fully developthe dominating mode. All these factors combined determinethe number of spiral seen in our simulations. It is impor-tant to realize that the most unstable radius according tothe axisymmetric instability criterion might not be the mosteffective site for swing amplification, since the shear playsan important role in this process.Without any external pumping source, spiral waves pro-duced by swing amplification should be a transient phe- nomenon. Similar to material spirals, swing amplified spi-ral waves will experience azimuthal shearing which reducestheir pitch angle until they become too tightly wound to beidentified. As an example, in the Gas2 simulation, the spiralarm that appears around R = 2 kpc initially has a pitchangle of 90 ◦ and should be sheared to less than 1 ◦ within2.2 Gyr. On the contrary, we find that the spontaneously in-duced spirals seen in Gas2 can last for more than 3 Gyr andstill keep the pitch angle relatively open. This result sug-gests at least one mechanism keeps replenishing noise intothe disk, leaving the physics to pick up the dominating modeand sustain the waves. This noise can be caused by numericsor preexisting waves. In this paper we have developed a simple and effectivemethod to compute the three-dimensional density and ve-locity structure of an isothermal gas disk in hydrodynamicequilibrium in the presence of an arbitrary external potential(i.e., dark matter halo and/or stellar disk). This is ideallysuited to set-up the initial conditions of a three-dimensionalgas disk in equilibrium in hydrodynamical simulations. We c (cid:13) , 1– ?? quilibrium Initialization and Stability of Three-Dimensional Gas Disks (c) Ga sStar Figure 7.
Images (a) to (d) correspond to models GasStar1 to GasStar4. They show the face-on surface density at t = 250 Myr. Thesize of the images are 20 kpc ×
20 kpc. The gas at the most unstable radius has orbited around the center about two and half times.Spirals seen in model GasStar2(b) and GasStar3(c) are due to swing amplification. In (d) the disk fragments very fast mainly due toboth the axisymmetric instability and swing amplification. first notice that as long as the gas is barotropic or has con-stant temperature at t = 0, the circular velocity needed tosupport the self-consistent disk is independent of the heightabove or below the midplane. This feature greatly simpli-fies the process of specifying the initial velocity field. All weneed to know is the rotation velocity in the midplane.To specify the density distribution self-consistently, thehydrostatic equation coupled with the reduced Poisson equa-tion is adopted to develop the vertical structure of the gas.Two sets of second-order non-linear differential equationsare found. One is directly associated with the gas densitycalled the density method, the other associated with the gaspotential called the potential method. In a simulation in-volving a huge dynamic range (using AMR techniques), thepotential method is shown to be numerically more stable.A simple local iteration can be performed to gain a bettercontrol on the shape and the mass of disks. These ideas aresimple enough to be incorporated into any existing code,and most importantly they are very effective.With gas disks that are in detailed balance, we are ableto systematically investigate the axisymmetric stability of a fully three-dimensional disk for the first time. We probe theonset of instability both semi-analytically and numerically.Simulations without stellar disk show that the thickness ofthe gas disk, which is governed by the temperature of thedisk, has a huge impact on the disk stability. The reduc-tion of the gravity decreases the threshold value by around30 percent in our models. As we gradually lower the gastemperature, the threshold Q th shifts up, the disk Q shiftsdown, and the system starts to develop multi-arm structurevia swing amplification. The onset of the instability in simu-lations matches the theoretical prediction very well as shownin Fig. 4 and Fig. 5. The disk fragments as the two curves, Q and Q th , come very close to each other.The influence of the stellar disk is less obvious. Its pres-ence has a stabilizing effect on the gas disk through changingthe rotation curve and a destabilizing one through the in-crease of the local gravitational force. The simulation resultsshow that overall the presence of the stellar disk tends to sta-bilize the gas disk. But this conclusion comes with a caveat.The interaction between live stars and gas might be impor-tant. A live stellar disk itself can be unstable or marginally c (cid:13) , 1– ?? Wang et al.
Figure 8.
The image size in (a) and (b) is 20 kpc ×
20 kpc. (a) The surface density of Gas2 at t = 750 Myr. (b) The surface density ofGasStar1 at t = 500 Myr. In both cases, the inner parts of the gas disks, which have been evolved for about four orbital times, developingspiral structure. Contour plots (c) and (d) are the Fourier maps of (a) and (d), respectively. In (c) and (d), the horizontal axis representsradius, the vertical axis is the number of arms, m , obtained by Fourier analysis. The color represents the intensity of each Fourier mode,the redder the stronger. stable. Perturbations from the interstellar medium can trig-ger instabilities in the stellar disk. Since stars dominate themass budget in Milky-Way type galaxies (more than 90 per-cent), and because gas is highly responsive and dissipative,the interplay between both components is one of the mostinteresting subjects in galactic dynamics. Tackling this prob-lem needs elaborate initial conditions for the live stellar diskor the combined disk. We stress that the potential methoddeveloped in this paper is compatible with the formulationin KD95. This makes the self-consistent combined disk anatural direction for future work.Marginally stable disks are susceptible to the processof swing amplification, a prevalent mechanism that triggersself-induced spirals. Simulations Gas2 and GasStar1 showthe spirals are prominent in the regions in which the gascan respond to swing amplification. Semi-analytic result re-lates the most vulnerable wavelength in azimuthal direction, λ max, to the number of arms. Numerically, The naturalmode of a Cartesian grid together with the swing amplifica-tion determine the dominating mode of the spiral structure.Our numerical results with or without stellar disk shows thecorrect characteristics of the swing amplification. It happensin marginally stable disks and the number of arm fits reason-ably well to the analytic prediction. In the run of GasStar2, swing amplification eventually leads to disk fragmentationonce the density becomes supercritical to the gravitationalinstability. However, in a sub-critical case like Gas2, the spi-ral structure can survive more than 3 Gyrs without frag-menting the disk, suggesting at least one mechansim is sus-taining the waves. The number of arms suggests a charac-teristic wavelength relating to the upper limit of the massof giant molecular clouds (Escala 2008). ACKNOWLEDGMENTS
R.S.K. acknowledges financial support from the German
Bundesministerium f¨ur Bildung und Forschung via the AS-TRONET project STAR FORMAT (grant 05A09VHA) andfrom the
Deutsche Forschungsgemeinschaft (DFG) undergrants no. KL 1358/1, KL 1358/4, KL 1359/5, KL 1359/10,and KL 1359/11. R.S.K. furthermore thanks for subsidiesfrom a Frontier grant of Heidelberg University sponsored bythe German Excellence Initiative and for support from the
Landesstiftung Baden-W¨urttemberg via their program Inter-national Collaboration II (grant P-LS-SPII/18 ). R.S.K. alsothanks the KIPAC at Stanford University and the Depart-ment of Astronomy and Astrophysics at the University of c (cid:13) , 1– ?? quilibrium Initialization and Stability of Three-Dimensional Gas Disks California at Santa Cruz for their warm hospitality during asabbatical stay in spring 2010. H.-H.W. acknowledges the fi-nancial support from the International Max Planck ResearchSchool (IMPRS) Heidelberg and the technical support withRAMSES from Dr. R. Teyssier. Authors acknowledge theanonymous referee’s comments and discussions that makethis paper better. Numerical simulations were performed onthe PIA cluster of the Max-Planck-Institut f¨ur Astronomieat the Rechenzentrum in Garching.
REFERENCES
Agertz O., Lake G., Teyssier R., Moore B., Mayer L., RomeoA.B., 2009, MNRAS, 392, 294Athanassoula E., Bosma A., Papaioannou S. 1987, A&A, 179, 23Binney J., Tremaine S. 2008, Galactic Dynamics (Princeton:Princeton Univ. Press, 2nd Ed.)Casertano S., 1983, MNRAS, 203, 735Escala A., Larson R.B., 2008, ApJ, 685, L31Evans N.W., Read J.C.A., 1998, MNRAS, 300, 106Fuchs B., 2001, A&A, 368, 107Fuchs B., 2008, Invited contribution to Galactic and Stellar Dy-namics in the era of high resolution surveys, Strassburg,March 16-20 [arXiv:0810.3503]Fuchs B., von Linden S., 1998, MNRAS, 294, 513Gammie C.F., 2001, ApJ, 553, 174Goldreich P., Lynden-Bell D., 1965, MNRAS, 130, 125Gradshteyn I. S., Ryzhik I. M., 1965, Table of Integrals, Series,and Products, 4th Ed. (Academic Press, New York)Hohl F. 1971, ApJ, 168, 343Jackson, P.D., Kellman, S.A., 1974, ApJ, 190, 53Julian W.H., Toomre A., 1966, ApJ, 146, 810 (JT)Kalnajs A.J., 1978, In IAU Symposium 77, Structure and Prop-erties of Nearby Galaxies, ed. Berkhuijsen E.M., WielebinskiR. (Dordrecht: Reidel), 133Kent, S. M., Dame, T. M., Fazio, G. 1991, ApJ, 378, 131Kim W.T., Ostriker E.C., 2007, ApJ, 660, 1232Kim W.T., Ostriker E.C., 2006, ApJ, 646, 213Kim W.T., Ostriker E.C., Stone J.M., 2002a, ApJ, 581, 1080Kim W.T., Ostriker E.C., 2002b, ApJ, 570, 132Kuijken K., Dubinski J., 1995, MNRAS, 277, 1341Leroy A.K., Walter F., Brink E., Bigiel F., de Blok W.J.G.,Modore B., Thornley M.D. 2008, AJ, 136, 2782Li X., Mac Low M.-M, Klessen R.S., 2005a, ApJ, 620, 19Li X., Mac Low M.-M, Klessen R.S., 2005b, ApJ, 626, 823Li X., Mac Low M.-M, Klessen R.S., 2006, ApJ, 639, 879Lin C.C, Shu F.H., 1964, ApJ, 140, 646Lisker T., Fuchs B., 2009, A&A, 501, 429Lockman, F. J. 1984, ApJ, 283, 90Miyamoto M., Nagai R., 1975, Astron. Soc. Japan, 27, 533Narayan C.A., Jog C.J., 2002, A&A., 394, 89Navarro J., Frenk C., White S.D.M., 1997, ApJ, 490, 493Pichon C., Cannon R.C., 1997, MNRAS, 291, 616Rafikov R. R., 2001, MNRAS, 323, 445Sanders, D. B., Solomon, P. M., Scoville, N. Z. 1984, ApJ, 276,182Sawamura M., 1988, PASJ, 40, 279Sellwood J.A., Moore E.M., 1999, ApJ, 510, 125Sellwood J.A., 1985, MNRAS, 217, 127Sellwood J.A., Carlberg R.G., 1984, ApJ, 282, 61Sellwood J.A., 1981, A&A, 99, 362Shetty R., Ostriker E.C., 2006, ApJ, 647, 997Shu F.H., 1969, 158, 505Spizter L., JR., 1942, ApJ, 95, 329Springel V., Matteo T.D., Hernquist L., 2005, MNRAS, 361, 776Tasker E., Bryan G., 2006, ApJ, 641, 878 Teyssier R., 2002, A&A, 385, 337Toomre A., 1990, in Dynamics and Interactions of Galaxies, ed.R. Wielen (Springer, Berlin), 292Toomre A., 1981, in Structure and Evolution of Normal Galaxies,ed. Fall, S.M., Lynden-Bell, D. (Cambridge: CUP), 111Toomre A., 1964, ApJ, 139, 1217van der Kruit, P. C., Searle, L. 1981a, A&A, 95, 105van der Kruit, P. C., Searle, L. 1981b, A&A, 95, 116Vauterin P., Dejonghe H., 1996, A&A, 313, 465Wouterloot, J. G. A., Brand, J., Burton, W. B., Kwee, K. K.1990, A&A, 230, 21Zang T.A., Hohl F. 1978, ApJ, 226, 521
APPENDIX A: THE DERIVATION OFROTATION VELOCITY
Equation (11) can be re-written as p ( R, z ) = ρ ( R, z ) p ( R ) ρ ( R ) (cid:12)(cid:12)(cid:12)(cid:12) z =0 − ρ ( R, z ) Z z pρ ∂ρ∂z d z − ρ ( R, z )[Φ(
R, z ) − Φ( R, z = 0)] , (A1)where we have replaced Φ z = Φ( R, z ) − Φ( R, z = 0). Insert-ing Eq. (A1) in Eq. (1) involves a partial derivative to theintegral, let us prepare it first: ∂∂R (cid:18)Z z pρ ∂ρ∂z d z (cid:19) = Z z ∂∂R (cid:18) pρ ∂ρ∂z (cid:19) d z = Z z (cid:26) ∂ ∂R∂z (cid:18) pρ ρ (cid:19) − ∂∂R (cid:20) ρ ∂∂z (cid:18) pρ (cid:19)(cid:21)(cid:27) d z = Z z (cid:26) ∂∂z (cid:20) pρ ∂ρ∂R + ρ ∂∂R (cid:18) pρ (cid:19)(cid:21) − ∂∂R (cid:20) ρ ∂∂z (cid:18) pρ (cid:19)(cid:21)(cid:27) d z = p ( R, z ) ρ ( R, z ) ∂ ln ρ ( R, z ) ∂R − p ( R, z = 0) ρ ( R, z = 0) ∂ ln ρ ( R, z = 0) ∂R + Z z (cid:26)(cid:16) ∂ρ∂z (cid:17) ∂∂R (cid:18) pρ (cid:19) − (cid:16) ∂ρ∂R (cid:17) ∂∂z (cid:18) pρ (cid:19)(cid:27) d z (A2)With Eq. (A2), the first term of Eq. (1) then becomes:1 ρ ( R, z ) ∂p ( R, z ) ∂R = ∂∂R (cid:18) p ( R, z = 0) ρ ( R, z = 0) (cid:19) − (cid:20) ∂ Φ( R, z ) ∂R − ∂ Φ( R, z = 0) ∂R (cid:21) + p ( R, z = 0) ρ ( R, z = 0) ∂ ln ρ ( R, z = 0) ∂R − Z z (cid:26)(cid:16) ∂ρ∂z (cid:17) ∂∂R (cid:18) pρ (cid:19) − (cid:16) ∂ρ∂R (cid:17) ∂∂z (cid:18) pρ (cid:19)(cid:27) d z + ∂ ln ρ ( R, z ) ∂R (cid:26) p ( R, z = 0) ρ ( R, z = 0) − p ( R, z ) ρ ( R, z ) − Φ( R, z ) + Φ(
R, z = 0) − Z z pρ ∂ρ∂z d z (cid:27) (A3)Equation (11) says that the term in the big brace shouldvanish. And therefore, Eq. (1) reduces to1 ρ ∂p∂R + ∂ Φ ∂R = ∂∂R (cid:20) p ( R, z = 0) ρ ( R, z = 0) (cid:21) + ∂ Φ( R, z = 0) ∂R + p ( R, z = 0) ρ ( R, z = 0) ∂ ln ρ ( R, z = 0) ∂R c (cid:13) , 1– ?? Wang et al. − Z z (cid:26)(cid:16) ∂ρ∂z (cid:17) ∂∂R (cid:18) pρ (cid:19) − (cid:16) ∂ρ∂R (cid:17) ∂∂z (cid:18) pρ (cid:19)(cid:27) d z (A4)For the barotropic gas, i.e., p ( ρ ), the integrand of the integralvanishes: (cid:16) ∂ρ∂z (cid:17) ∂∂R (cid:18) pρ (cid:19) − (cid:16) ∂ρ∂R (cid:17) ∂∂z (cid:18) pρ (cid:19) = (cid:16) ∂ρ∂z (cid:17) ∂∂ρ (cid:18) pρ (cid:19) ∂ρ∂R − (cid:16) ∂ρ∂R (cid:17) ∂∂ρ (cid:18) pρ (cid:19) ∂ρ∂z = 0(A5)For the cases of initially constant temperature, the specificinternal energy, e , is a constant and therefore the pressureis a function of density only, the integrand vanishes. APPENDIX B: THE EFFECT OF THE DISKTHICKNESS ON THE MIDPLANE POTENTIAL
For an axisymmetrically and infinitesimally thin disk, thepotential can be evaluated by the following relation (Binney& Tremaine, 2008):Φ(
R, z ) = Z ∞ d kS ( k ) J ( kR ) e − k | z | , (B1)where J is the Bessel function of the first kind of order zeroand S is the Hankel transform of − πG Σ defined by: S ( k ) = − πG Z ∞ d R ′ R ′ J ( kR ′ )Σ ( R ′ ) (B2)With Eq. (B1) and Eq. (B2), we can superimpose the poten-tial contributed by each gas layer. For the sake of simplicity,we assume that the volume density has the double exponen-tial profile: ρ ( R, z ) = Σ e − R/R d e − z/h z h z , (B3)with h z being the scale-height of the gaseous disk, Eq. (B2)then becomes (Gradshteyn & Ryzhik 1965, hereafter GR65,6.623-2): S ( k, z ) = − πG Σ e − z/h z h z ∆ z Z ∞ d R ′ R ′ J ( kR ′ ) e − R ′ /R d = − πG Σ e − z/h z h z ∆ z ξ ( ξ + k ) / , (B4)with ξ = 1 /R d . ∆ z represents the infinitesimal thickness in-troduced to keep the dimension correct. The potential whichtakes into account the thickness of the disk then reads:Φ( R, z ) = − πG Σ Z ∞ d k ξ ( ξ + k ) / J ( kR ) × Z ∞−∞ e − k | z − h | e − h/h z h z d h. (B5)Evaluating the potential at the midplane, z = 0, yields:Φ( R, z = 0) = − πG Σ Z ∞ d k ξ ( ξ + k ) / J ( kR ) 11 + kh z . (B6)Given the finite scale-height, the integral can be eval-uated numerically and compared with the result of the in-finitesimally thin disk. F z , d m / F z , ga s r = 2kpcr = 5kpcr = 8kpc Figure D1.
The force ratio F z, DM /F z, gas at R=2, 5 and 8 kpc.It shows that the vertical structure of the inner disk is determinedmostly by the self-gravity of gas. APPENDIX C: THE DERIVATION OF THEREDUCTION FACTOR
To derive the reduction factor F defined by Eq. (35) we needto evaluate the integral of the form: Z ∞−∞ e − k | h | sech ( ah )d h = 2 Z ∞ e − kh sech ( ah )d h = 2 a n − k a h H (cid:16) k a (cid:17) − H (cid:16) k a − (cid:17)io . (C1)The last line can be reached by looking up the formulae3.541, 8.370, 8.361-7 listed in the integral table (GR65) andthe definition Eq. (36). In the last line, we have employedthe recursive relation (8.365-1 GR65): H ( α ) = H ( α −
1) + 1 α . (C2)The asymptotic behavior of the harmonic number reads(8.367-2, 8.367-13 GR65): H ( α ) = ln α + γ + 12 α − − α − + 1120 α − + O ( α − ) , (C3)with γ = 0 . α ≥
1. We employ the recursive relation (C2) to eval-uate H ( α ) for − < α < APPENDIX D: THE VERTICAL FORCE RATIO
The vertical force ratio measures the impact of the halo forceon the vertical structure. The simplified Poisson equation forisothermally self-gravitating gas disk reads: ∂ Φ g ∂z = 4 πGρ ( R )sech (cid:16) zh z (cid:17) , (D1)where h z being a measure of the scale-height. Parameter h z can be related to the volume density in the midplane, ρ ( R )by: h z = r c πGρ . (D2) c (cid:13) , 1– ?? quilibrium Initialization and Stability of Three-Dimensional Gas Disks radius [Kpc] he i gh t [ K p c ] Figure E1.
Contour map of ǫ . The black lines represent thescale-height of the gas disk. The corresponding vertical force for the gas then becomes: F z, gas = − ∂ Φ ∂z = − πGh z ρ ( R )tanh( z/h z ) . (D3)For a NFW halo, the vertical force can be written downdirectly: F z, DM = GM f ( c ) (cid:16) cr (cid:17) x/ (1 + x ) − ln(1 + x ) x z √ R + z , (D4)with x = c √ R + z /r . Figure D1 then shows the forceratio F z, DM /F z, gas as a function of vertical height | z | at dif-ferent radii. Comparing to the rotation curve shown in theleft panel of Fig. 1, although the dynamics is still dictatedby the potential of the dark halo, the vertical structure ofthe gaseous disk is mainly determined by the self-gravityof the gas component. However, it raises another issue, thepresence of the stellar disk will dominate both the dynam-ics and the vertical structure of the gas and will affect thestability of the gas component via changing the thickness ofthe gaseous disk and the rotation curve. APPENDIX E: VALIDITY CHECK OF THEREDUCED POISSON EQUATION FOR THEGAS DISK
Throughout this paper we have assumed that the radial po-tential gradients of the disk are negligible compared to thevertical gradients, such that the Poisson equation reducesto Eq. (14). We now test this assumption by computing theratio ǫ ≡ | R ∂∂R ( R ∂ Φ g ∂R ) / ∂ Φ g ∂z | , (E1)with Φ g the gravitational potential of the gas disk. For arealistic, analytical disk model, our assumption will be validas long as ǫ ≪ g ( R, z ) = − GM g q R + ( a + p z + b ( R )) . (E2) Here a is a constant that controls the scale-length of thedisk and b ( R ), which we take to be a function of radius,modulates the scale-height of the disk. In the limit b → a = 3 . b ( R ) = − . × − R + 1 . × − R + 0 . . (E3)Using the Poisson equation to solve (numerically) forthe corresponding density distribution yields the radial-dependent scale-height shown as the solid black lines inFig. E1, and which is comparable to that of the Gas0 disk.The contours in Fig. E1 are defined by constant values of ǫ .These show that our assumption that ǫ ≪ ∼ ǫ ≪ ǫ that are even smaller. Based on theseresults, and based on the absence of significant disk thick-ening in our simulations, we are confident that Eq. (14) issufficiently accurate for all realistic gas disks. c (cid:13) , 1–, 1–