On the conservation of the vertical action on galactic disks
DDraft version May 19, 2016
Preprint typeset using L A TEX style emulateapj v. 5/2/11
ON THE CONSERVATION OF THE VERTICAL ACTION IN GALACTIC DISKS
Carlos Vera-Ciro , Elena D’Onghia Draft version May 19, 2016
ABSTRACTWe employ high-resolution N -body simulations of isolated spiral galaxy models, from low-amplitude,multi-armed galaxies to Milky Way-like disks, to estimate the vertical action of ensembles of starsin an axisymmetrical potential. In the multi-armed galaxy the low-amplitude arms represent tinyperturbations of the potential, hence the vertical action for a set of stars is conserved, although afterseveral orbital periods of revolution the conservation degrades significantly. For a Milky Way-likegalaxy with vigorous spiral activity and the formation of a bar, our results show that the potentialis far from steady, implying that the action is not a constant of motion. Furthermore, because of thepresence of high-amplitude arms and the bar, considerable in-plane and vertical heating occurs thatforces stars to deviate from near-circular orbits, reducing the degree at which the actions are conservedfor individual stars, in agreement with previous results, but also for ensembles of stars. If confirmed,this result has several implications, including the assertion that the thick disk of our Galaxy formsby radial migration of stars, under the assumption of the conservation of the action describing thevertical motion of stars. Subject headings: galaxies: kinematics and dynamics - Galaxy: disk - Galaxy: evolution - stars:kinematics and dynamics INTRODUCTION
A stellar system can be fully modeled by a phase-spacedensity distribution function (DF), whose evolution is de-scribed by solutions of the Boltzmann equation (Binney& Tremaine 2008). Is also a well-known result that forsteady-state potentials, the DF describing the galaxy canbe represented as a function of the integrals of motionof the system (Jeans 1919). There are several choicesone can make to use then as coordinates of the DF, e.g.energy and angular momentum (Eddington 1915). In-tegrated orbits in axisymmetric potentials indicate theexistence of a non-classical isolating integral of motion(Ollongren 1962), which is closely related to the actionsof the system (Binney & Spergel 1984).Besides being a natural choice to describe the evolutionin the phase-space of collisionless stellar systems, the ac-tions are also adiabatic invariants. This is a very appeal-ing property of spiral galaxies, where the high frequencyof the vertical motion, compared to in-plane evolution,makes the action describing the z motion a prime candi-date for being a conserved quantity, even in the presenceof spiral arms. This assumption has been tested by Sol-way et al. (2012), who showed that for isolated galaxiesthe vertical action is not a constant of motion of individ-ual stars, and is only conserved in average for samples ofstars, with an intrinsic dispersion of ∼ Department of Astronomy, University of Wisconsin, 2535Sterling Hall, 475 N. Charter Street, Madison, WI 53076, USA.e-mail:[email protected] Alfred P. Sloan Fellow ble increase of random motion, causing a process termedradial migration. Following studies confirmed that thisprocess occurs in 3D simulations of galactic disks eitherin isolation or in a cosmological context (Roˇskar et al.2008; Minchev et al. 2011; Grand et al. 2012; Loebmanet al. 2011; Kubryk et al. 2014; Halle et al. 2015; Grandet al. 2015).If radial mixing in galactic disks occurs, this processhas interesting implications for the spread in the age-metallicity relation observed in the solar neighborhood(Edvardsson et al. 1993). If stellar migration never oc-curred, the stars in the solar vicinity would have metallic-ities ranging from zero to the present-day value (Serenelliet al. 2009) with a significant correlation between the ageof the stars and their metallicities. However this correla-tion has never been found in the current data, raising thepossibility that the stars migrated from their birth radii(Wielen et al. 1996). Because the radii of stars at thebirth depends on their specific angular momentum, stel-lar migration occurs if there is a change of the angularmomentum of the stars. Changes in the angular mo-mentum for individual stars can be induced by scatter-ing with giant molecular clouds (Spitzer & Schwarzschild1953; Mihalas & Binney 1981) or transient spiral struc-tures (Jenkins & Binney 1990; Jenkins 1992), interfer-ence of coherent long-lived modes and density waves witha bar (Roˇskar et al. 2012), or by interactions with satel-lite galaxies (Kazantzidis et al. 2008; Villalobos & Helmi2008; Bird et al. 2012).In this paper we aim to study in detail the effect ofmigration on isolated disks, by using N -body simula-tions for disk galaxies evolving in isolation, to achievehigh-resolution; with the aim of better determining thecontribution of the spiral waves to the heating of the stel-lar disk, to the radial migration of stars in the disk, andultimately to the formation of a thick disk.This manuscript is laid out as follows. In Section 2 weintroduce the N -body simulations used in our work and a r X i v : . [ a s t r o - ph . GA ] M a y Vera-Ciro & D’Onghiadescribe the method to measure the vertical action. InSection 3 we discuss to what extent this quantity is con-served during the evolution of the disk, and in Section 4we show how this can be used to study the formation of athick disk. Finally in Section 5 we draw our conclusions. NUMERICAL PRELIMINARIES
Simulations
For this study we employed N -body simulations oftwo galaxy models labeled as disk 1 and disk 2 . Eachgalaxy consisted of a static dark matter halo and a liverotationally supported stellar disk of 5 × particles.Because the spiral morphology in stellar disks dependsin detail on the mass distribution of the galaxy andon the self-gravity of the disk, we selected the struc-tural parameters of the galaxies such that one galaxydeveloped a low-amplitude multi-armed disk ( disk 1 ),whereas the other developed a spiral morphology moresimilar to a Milky Way (MW)-like galaxy ( disk 2 ). Thiswas achieved by choosing models with very different diskmass fractions within 2 . M disk /M total =0 .
271 for disk 1 and 0 .
478 for disk 2 , respectively. Byvarying the mass distribution of the disk within 2 . λ crit = 4 π G Σ /κ , where Σ is thedisk surface mass density and κ the epicycle frequency,and lead to a different spiral morphology (Toomre 1981;Sellwood & Carlberg 1984; Carlberg & Freedman 1985;D’Onghia 2015).The parameters describing each galaxy componentwere assumed to be independent and both models wereconstructed similarly to the approach described in (Hern-quist 1993; Springel et al. 2005). The dark matter massdistribution was modeled assuming the Hernquist profile(Hernquist 1990): ρ halo ( r ) = M halo πR r/R halo )(1 + r/R halo ) , (1)with M halo being the total halo mass and R halo the scaleradius. The effect of the dark halo is modeled assuming astatic potential, this is of course a limitation of our modelsince the system never completely relaxes, however, thecomplication of adding an alive component comes withthe price of including further scattering to the particlesof the disk. The stellar disk is modeled in the initialconditions as a thin exponential surface density profileof scale-length R disk and total mass M disk . The verticalmass distribution of the stars in the disk is specified bythe profile of an isothermal sheet with a radially constantscale height z disk : ρ disk ( R, z ) = M disk πR z disk e − R/R disk sech zz disk . (2)Model disk 2 included a static potential for the bulgedescribed by the Hernquist model with mass M bulge andscale-length R bulge . We aim to build this galaxy modelwith a flat rotation curve consistent with the observa-tions of the MW within two optical radii ( R (cid:46)
14 kpc)(Reid et al. 2014; Bovy et al. 2009). We noted that thisrequirement is satisfied by adopting a Hernquist profilefor the dark halo with a large scale radius, and a large
TABLE 1Structural properties of the galaxy models disk 1 disk 2 M disk (10 M (cid:12) ) 1.905 4.000 R disk (kpc) 3.130 2.500 z disk /R disk M halo (10 M (cid:12) ) 0.933 9.500 R halo (kpc) 29.775 130.0 M disk /M total (2 . R disk ) 0.271 0.479 M bulge (10 M (cid:12) ) . . . 1.400 R bulge (kpc) . . . 0.350 total mass. However, the dynamical range of interest ofthis study is a few optical radii, thus the mass distri-bution outside this radius does not affect our analysis.Table 1 summarizes the structural parameters adoptedfor each galaxy model in this study.Fig. 1 shows the azimuthally averaged radial profilesfor the total circular velocity V circ (stars and dark mat-ter), the Toomre parameter Q Toomre of the stellar disk,the disk surface density Σ and its vertical velocity dis-persion σ z . The profiles are displayed at four differenttimes, spanning 6 Gyr of evolution. In what follows weuse the notation t n to refer to the time of integrationin our simulations in Gyr, with t thus labeling the theinitial time, and t corresponding to 6 Gyr.We note that in the multi-armed galaxy ( disk 1 model), Q Toomre increases with time, as a consequenceof the increase of the radial velocity dispersion. This isexpected, since in-plane heating occurs as a consequenceof the spiral structure activity, while the vertical struc-ture remains unchanged.This is, however, not the case for our disk 2 run. Abar instability grows at the center of the galaxy after ∼ ∼ disk 1 model grows a large numberof low-amplitude spiral arms that self-perpetuate for atleast 6 Gyr (D’Onghia et al. 2013). After 1 Gyr disk 2 already developed four high-amplitude spiral arms at twoscale-lengths ( ∼ ∼ disk 1 model) has been recently studied by Vera-Ciroet al. (2014). As first shown in Sellwood & Binney (2002),the authors confirmed that spiral activity in a low-massdisc causes stars to diffuse radially, with stars near coro-tation of a spiral pattern exchanging angular momentumand moving to a new radius without adding random mo-tion (see also Daniel & Wyse 2015). Indeed, the mech-n the conservation of the vertical action 3 Fig. 1.—
Radial profiles of total circular velocity, V cir , Toomreparameter, Q Toomre , surface mass density, Σ, and vertical velocitydispersion, σ z , for the two studied galaxy models: disk 1 (left) and disk 2 (right). t n refers to the time of integration in our simula-tions in Gyr, with t being the the initial time and t correspondingto 6 Gyr. anism of scattering of stars at corotation predicts a lackof significant disk heating and leaves the surface densityprofile unchanged. These two features are reproduced inour disk 1 model as shown in Fig. 2 (left panels). Inthe MW-like galaxy model ( disk 2 ) the bar formationcauses greater angular momentum changes than thoseof a multi-armed spiral structure (Minchev & Famaey2010). Furthermore, because the bar persists after itsformation, the associated angular momentum changesmight differ from those caused by recursive spiral arms.Therefore the surface density profile of the disk is ex-pected to change. Some disk heating is also expected tooccur with time. Measuring the vertical action
The action-angle variables can be used to describe theevolution of orbits in static potentials. The actions, j ,are constants of motion in the unperturbed field, or whenthe potential changes slowly so that they are adiabaticinvariants (Arnold 1978; Goldstein et al. 2002). The ver-tical action is defined as: j z ≡ π (cid:90) d z d v z = 12 π (cid:73) d z v z , (3)where z and v z describe the position and velocity alongthe orbit, respectively. There are a number of ways toestimate of the vertical action of a star particle in a Fig. 2.—
Time sequence of density projections of the face-onviews of stellar disks presented in this study. The disk 1 simula-tion develops a large number of low-amplitude spiral structure thatperpetuate for at least 6 Gyr (left panels); disk 2 develops approx-imately a four-fold rotational symmetry and becomes bar-unstableafter ∼ nearly axisymmetric potential. Solway et al. (2012), forinstance, calculated the quadrature in Eq. (3) by mea-suring the area enclosed by the orbit in the z − v z plane.We adopted here the prescription as described in Ap-pendix A. This method allows us to identify and rejectstars trapped in resonances, a case where the adiabaticinvariance breaks (Pfenniger 1984; Pfenniger & Friedli1991).In Fig. 3 we show an example for a random particle inour disk 1 simulation at two different times. The greenpoints show the surface of section of the integrated orbitand the red line is the connected loop used to calculateEq. (3). The blue line represents the zero-velocity curve,and gray points are other orbits, shown here to depictthe structure of the phase-space in our runs. The areadefined by the red-loop in Fig. 3 is the vertical action j z ,for this particular example, the vertical action measuredat t is a factor of ∼
80 larger than it is at time t , thisis clear evidence that the value of j z is not conserved forindividual orbits for slow-varying, low-amplitude pertur-bations, such as the ones developed in our disk 1 run. A CONSERVED QUANTITY?
Our previous calculations demonstrated that the ver-tical action j z is not a conserved quantity for individ- Vera-Ciro & D’Onghia Fig. 3.—
Surface of section for a star particle measured at twodifferent times. In the left panel the simulation has already devel-oped spiral arms, the right panel is the surface of section after 2Gyr. In this specific example, R gui ( t ) ≈ R gui ( t ) indicates thatthe particle has not migrated significantly. However, its verticalaction has changed by a factor ∼
80 ( δj z = 4 . Fig. 4.—
Probability distribution of fractional changes of thevertical action j z , as a function of changes in guiding radius R gui for particles selected at 1 Gyr with R gui = 5 . ± .
25 kpc. Thedistribution of fractional changes in j z does not depend on howefficiently particles migrate. The median and 1 σ equivalent disper-sion are represented with solid lines in the left panels. The rightpanels show the marginalized distribution of changes in the verticalaction. ual stars. This result agrees with findings of previousstudies (Solway et al. 2012). We investigated this prop-erty in our disk 1 and disk 2 simulations. After 1 Gyrof evolution an ensemble of stars with guiding radius: R gui = 5 . ± .
25 kpc was selected and the relativechanges in the vertical action were computed at a latertime (See Appendix A). This radius is close to two diskscale-lengths in both galaxy models. The vertical action changed significantly over a period of time as shown inFig. 3. Thus, its change from time t to a later time t (cid:48) > t was quantified in the following way: δj z ≡ ln j z ( t (cid:48) ) j z ( t ) . (4)In addition, we noted that over the same period of timea significant radial migration occurs, with a consequentchange of the guiding radii. Thus, a similar definitionwas also applied to δR gui . Fig. 4 (top panel) displays theoutcome of our study for the multi-armed galaxy disk1 over a period of 2 Gyr. In agreement with previousstudies, our findings suggest the vertical action is a con-served quantity on average for an ensemble of stars. Notethat the modest fractional changes of δj z in our simula-tion match the values reported by Solway et al. (2012) Furthermore, Fig. 4 indicates that any variation in thevertical action, δj z , is independent of the radial migra-tion, measured by changes in guiding radii, δR gui . Next,we computed the median of the fractional change in thevertical action µ ( δj z ) and the 1 σ equivalent dispersionaround the median σ ( δj z ). We obtained µ ( δj z ) = 0,with standard deviation σ ( δj z ) = 0 .
28 (displayed as thesolid black line in Fig. 4), also in agreement with thevalues reported in Solway et al. (2012).The same analysis performed on the more massive disk2 produces different results as shown in Fig. 4 (bottompanel). The median of fractional changes in the verticalaction is µ ( δj z ) = 0 .
09 with a dispersion σ ( δj z ) = 0 . µ ( δj z ) and its dispersion σ ( δj z ) of fractional changes in the vertical action are dis-played in Fig. 5 for the two galaxy models. Ensembles ofstars are selected at after 1 Gyr of evolution ( t = t ) withguiding radii R gui = 3 , , µ ( δj z ) and dis-persion σ ( δj z ) of fractional changes in the vertical actionagainst the corresponding variation in the velocity dis-persion of the vertical component.To test for the origin of the changes in the average ac-tion we note that the vertical energy of samples of stars ln x ≈ x −
1, for x ∼ n the conservation of the vertical action 5 Fig. 5.—
Left panels: evolution of average changes in the verticalaction (top panels) and its 1 σ equivalent dispersion (bottom pan-els) as a function of time. Right panels: the same quantities as inthe left panels are plotted against the changes in vertical velocitydispersion for ensembles of stars selected at different radii. proportionally depends on the vertical velocity dispersion (cid:104) j z (cid:105) ∼ σ z (cf. Eq. (7)), and so if the average vertical ac-tion increases as a consequence of the vertical heating inthe disk, then δj z = 2 δσ z . The multiplying factor in thisapproximation is the vertical epicycle frequency, whichremains mostly constant in time (See Fig. 1). Fig. 5shows that this is not the case, and that the increasein dispersion of the changes in the vertical action is notmerely determined by the vertical heating of the disk. IMPLICATIONS FOR THE VERTICAL STRUCTURE
What does the conservation of the vertical action im-ply for the structure and evolution of the stellar disk? Ithas been argued that under the assumption of the conser-vation of the vertical action j z stars migrating outwardsin the disk tend to move far from the mid-plane andpopulate the thick disk (Roˇskar et al. 2008; Sch¨onrich &Binney 2009). We address this question here.The epicycle approximation regulates the motion ofstars on near-circular orbits. In the approximation thevertical motion of a star is also decoupled from the in-plane motion. The motion perpendicular to the plane isthen described by the Hamiltonian (Binney & Tremaine2008): H z = 12 v z + 12 ν z . (5)For an ensemble of star particles at radius R , the av-erage vertical energy is (cid:104) H z (cid:105) = 12 σ z + 12 ν (cid:104) z (cid:105) , (6)where σ z is the vertical velocity dispersion of the ensem-ble. Note that σ z in Eq. (6) does not necessarily refer Fig. 6.—
Top: average vertical action as a function of radius, thesolid lines are the values calculated using the method described inthe previous section (Eq. (3)) whereas the dashed lines correspondto the epicycle approximation (Eq. (6)). The difference betweenthese two values is ∼ z -coordinates (solid), and by solvingEq. (7). to the value of velocity dispersion of all stars with homeradii at R . Indeed, Eq. (6) holds also for any popula-tion of stars whose guiding radii have changed over time(e.g. stars migrating from inner or outer locations withinthe disk) and that are used to calculate it. Thus, in theepicycle approximation, the vertical action is the follow-ing j z = H z /ν (Binney & Tremaine 2008): (cid:104) j z, epi (cid:105) = (cid:104) H z (cid:105) ν = 12 ν ( σ z + ν α h ) , (7)where, α h = (cid:104) z (cid:105) , (cid:104) z (cid:105) is the second moment of thevertical component, and α = π/ √
12. With this nor-malization, the parameter h is equivalent to the diskscale height z disk , h ≡ z disk , when the second momentis calculated for all stars sampled from the sech z/z disk model used for the vertical density of our simulations (cf.Eq. (2)).Once the average vertical action and vertical velocitydispersion are computed for an ensemble of stars, Eq. (7)is solved for the thickness h of that population of stars.Fig. 6 displays the average value of the vertical action j z as a function of the galactic radius at three differenttimes for both galaxy models (top panels). The solidlines show the average vertical action calculated for a setof stars using the algorithm described in Section 2. Thedashed lines illustrate the predicted vertical action usingthe epicycle approximation for that set of stars. Clearly,the epicycle approximation predicts values of the verticalaction 15% higher than the estimates obtained with thenumerical procedure. A similar trend is observed in thenumerical experiments described in Solway et al. (2012). Vera-Ciro & D’Onghia Fig. 7.—
Scale height h profile for particles in circular orbitsthat migrate from R gui ( t ) = 5 . ± .
25 kpc (left) and R gui ( t ) =8 . ± .
40 (right) as a function of their final guiding radii for thethe multi-arm simulation (top) and the MW-like model (bottom).The scale height profile of the whole disk at t is shown by the solidblue line. In green we show the thickness profile of the migratorsat t and in red the thickness at t . The orange dashed line isthe result of using Eq. (7) and assuming that the vertical action isconserved (cid:104) j z, epi ( t ) (cid:105) = (cid:104) j z, epi ( t ) (cid:105) . Then, Eq. (7) is used to compute the thickness of thedisk h at each radius, and the values obtained for theboth galaxy models are shown in Fig. 6 (bottom pan-els). A similar discrepancy is reflected in the estimate ofthe scale height at different radii. The values derived inthe epicycle approximation are 5% larger than the valuesobtained by evaluating Eq. (3) as previously described.In both simulations, after the disks evolved for 1Gyr ( t ), we selected a set of stars with guiding radii R gui = 5 . ± .
25 kpc (left panels in Fig. 7) and R gui = 8 . ± .
40 kpc (right panels), we follow thesestars for two billion years longer ( t ) and calculated theirscale height h profile as a function of their final positionin the disk. This would show if a subsample of migratingstars can change their thickness h as they migrate.We only select stars that are in nearly circular orbits L z /L circ > .
97, where L z is the vertical angular mo-mentum and L circ is the angular momentum that a starhas, should it move in a circular orbit with the sameenergy. By adding this constraint at both times t and t we ensure that the sample contains stars that migratewhile preserving their circular orbits (Sellwood & Binney2002).For this set of particles the average vertical action inthe epicycle approximation is also j z is computed atboth times t and t . The thickness at both times isalso included, in general stars that move toward the cen-ter tend to reduce their scale height, whereas stars thatmigrate outwards increase it. For comparison we alsoshow in blue the scale height of the disk as a functionof radius, calculated with the normalized RMS value of the z -coordinates of all particles in the disk. This valuematches the one reported as z disk in table 1 (blue solidline).The bias reported in (Vera-Ciro et al. 2014) is observedfor both disks. This explains the fact that for stars thatmigrate outwards, the slight increase in their verticalscale height is not enough to create a thicker componentat their final location. We also included in Fig. 7 thepredicted thickness h at final t under the assumptionthat the action is conserved. It can clearly be seen thatfor our disk 1 no thick distribution of stars is either ob-served or predicted even when the average vertical actionis clearly conserved.In the MW-like galaxy, disk 2 , a bar grows in the in-ner parts of the disk, in qualitative agreement with sim-ulations of Brunetti et al. (2011). As a consequence ofthe bar formation and the vigorous spiral activity thereare not only significant changes of angular momentum,but also a considerable in-plane heating that forces starsin this simulation to deviate from a near-circular orbits.There are good reasons to believe that in this case the en-ergy of vertical motion is clearly not decoupled from thatof the radial part of the motion. Indeed, the epicycle ap-proximation is a poor description of the motion of mostparticles, for which the radial and vertical oscillations areneither harmonic nor are the energies of the two oscilla-tions decoupled. Furthermore, the actions are constantsof motion under the conditions of slow changes of thepotential. In the multi-armed galaxy the low-amplitudearms represent tiny perturbations of the potential, hencethe vertical action might still be conserved, although af-ter several orbital periods of revolution the vertical ac-tion is not conserved at the same level of accuracy. How-ever, in the MW-like galaxy the high-amplitude spiralstructure and the bar formation seem to compromise thevalidity of these conditions. SUMMARY AND CONCLUSIONS
We have presented quantitative estimates of the verti-cal action of ensembles of stars in the presence of spiralactivity. The vertical action of star particles has beenevaluated in isolated disks of galaxies with different spi-ral morphologies. The two galaxy models consist of alow-mass multi-armed disk ( disk 1 ) and a MW-like diskthat develops a bar after ∼ . disk2 ).For the multi-armed disk galaxy, disk 1 , we find thatthe vertical action for an ensemble of stars is approxi-mately a conserved quantity, with a dispersion of frac-tional changes around its median of 20%, in agreementwith the results reported in the literature. However, astime progresses the action is conserved with less accu-racy, as indicated by the increase of the standard devi-ation up to the value of 50% after 6 Gyr of evolution.This also holds for a subset of stars that radially mi-grates by resonant scattering at corotation. These re-sults have implications for the vertical structure of thedisk when stars radially migrate to the outer part of thedisk. In this model stars migrating outwards are a heav-ily biased subset of stars with preferentially low verticalvelocity dispersion. Thus, extreme migrators outwardsare characterized by having a have small amplitude oftheir vertical excursion, independent of the conservationof the vertical action.n the conservation of the vertical action 7In the MW-like galaxy model, disk 2 , the bar for-mation causes large angular momentum changes, withsignificant in-plane and vertical heating. Thus, in thepresence of the formation of the bar and combined withsignificant disk heating, our results show that the verti-cal action cannot be assumed a constant of motion. In-deed, in this case the fractional changes in the estimatesof the action for a set of stars progressively increase bytime, with a standard deviation that approaches 80% af-ter 6 Gyr of disk evolution. If confirmed, this result castsdoubts on the possibility of predicting the vertical struc-ture of the disk for a set of extreme migrators outwards,by assuming that their vertical action is conserved. This is funded by NSF Grant No AST-1211258 andATP NASA Grant No NNX144AP53G. ED gratefullyacknowledges the support of the Alfred P. Sloan Founda-tion. We express our appreciation toward the Aspen Cen-ter for Physics for their hospitality, funded by the NSFunder Grant No. PHYS-1066293. Simulations have beenrun on the High Performance Computing cluster pro-vided by the Advanced Computing Infrastructure (ACI)and Center for High Throughput Computing (CHTC)at the University of Wisconsin. We thank Daniel Pfen-niger, Jerry Sellwood, Martin Weinberg, Victor Debat-tista, Rob Grand and Jonathan Bird for the very helpfulcomments and constructive discussions. REFERENCESArnold, V. I. 1978, Mathematical methods of classical mechanics(Springer)Athanassoula, E., & Mart´ınez-Valpuesta, I. 2008, in AstronomicalSociety of the Pacific Conference Series, Vol. 390, PathwaysThrough an Eclectic Universe, ed. J. H. Knapen, T. J.Mahoney, & A. Vazdekis, 454Binney, J., & Spergel, D. 1984, MNRAS, 206, 159Binney, J., & Tremaine, S. 2008, Galactic Dynamics: SecondEdition (Princeton University Press)Bird, J. C., Kazantzidis, S., & Weinberg, D. H. 2012, MNRAS,420, 913Bovy, J., Hogg, D. W., & Rix, H.-W. 2009, ApJ, 704, 1704Brunetti, M., Chiappini, C., & Pfenniger, D. 2011, A&A, 534, A75Candy, J., & Rozmus, W. 1991, Journal of ComputationalPhysics, 92, 230Carlberg, R. G., & Freedman, W. L. 1985, ApJ, 298, 486Daniel, K. J., & Wyse, R. F. G. 2015, MNRAS, 447, 3576D’Onghia, E. 2015, ApJ, 808, L8D’Onghia, E., Vogelsberger, M., & Hernquist, L. 2013, ApJ, 766,34Dormand, J., & Prince, P. 1980, Journal of Computational andApplied Mathematics, 6, 19Eddington, A. S. 1915, MNRAS, 76, 37Edvardsson, B., Andersen, J., Gustafsson, B., et al. 1993, A&A,275, 101Goldstein, H., Poole, C., & Safko, J. 2002, Classical mechanics(Addison-Wesley)Grand, R. J. J., Kawata, D., & Cropper, M. 2012, MNRAS, 421,1529—. 2015, MNRAS, 447, 4018Halle, A., Di Matteo, P., Haywood, M., & Combes, F. 2015,A&A, 578, A58Hernquist, L. 1990, ApJ, 356, 359—. 1993, ApJS, 86, 389Jeans, J. H. 1919, Problems of cosmogony and stellar dynamics(Cambridge University Press)Jenkins, A. 1992, MNRAS, 257, 620Jenkins, A., & Binney, J. 1990, MNRAS, 245, 305Kazantzidis, S., Bullock, J. S., Zentner, A. R., Kravtsov, A. V., &Moustakas, L. A. 2008, ApJ, 688, 254Kubryk, M., Prantzos, N., & Athanassoula, E. 2014, preprint,arXiv:1412.0585 Loebman, S. R., Roˇskar, R., Debattista, V. P., et al. 2011, ApJ,737, 8Martinez-Valpuesta, I., Shlosman, I., & Heller, C. 2006, ApJ, 637,214Merritt, D., & Sellwood, J. A. 1994, ApJ, 425, 551Mihalas, D., & Binney, J. 1981, Galactic astronomy: Structureand kinematics /2nd edition/ (W H Freeman & Co (Sd))Minchev, I., & Famaey, B. 2010, ApJ, 722, 112Minchev, I., Famaey, B., Combes, F., et al. 2011, A&A, 527, A147Ollongren, A. 1962, Bull. Astron. Inst. Netherlands, 16, 241Pettitt, A. R., Dobbs, C. L., Acreman, D. M., & Bate, M. R.2015, MNRAS, 449, 3911Pfenniger, D. 1984, A&A, 134, 373Pfenniger, D., & Friedli, D. 1991, A&A, 252, 75Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783,130Roˇskar, R., Debattista, V. P., Quinn, T. R., Stinson, G. S., &Wadsley, J. 2008, ApJ, 684, L79Roˇskar, R., Debattista, V. P., Quinn, T. R., & Wadsley, J. 2012,MNRAS, 426, 2089Sch¨onrich, R., & Binney, J. 2009, MNRAS, 399, 1145Sellwood, J. A., & Binney, J. J. 2002, MNRAS, 336, 785Sellwood, J. A., & Carlberg, R. G. 1984, ApJ, 282, 61Serenelli, A. M., Basu, S., Ferguson, J. W., & Asplund, M. 2009,ApJ, 705, L123Solway, M., Sellwood, J. A., & Sch¨onrich, R. 2012, MNRAS, 422,1363Spitzer, Jr., L., & Schwarzschild, M. 1953, ApJ, 118, 106Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 361,776Toomre, A. 1981, in Structure and Evolution of Normal Galaxies,ed. S. M. Fall & D. Lynden-Bell, 111–136Vera-Ciro, C., D’Onghia, E., Navarro, J., & Abadi, M. 2014, ApJ,794, 173Villalobos, ´A., & Helmi, A. 2008, MNRAS, 391, 1806Wielen, R., Fuchs, B., & Dettbarn, C. 1996, A&A, 314, 438Yurin, D., & Springel, V. 2015, MNRAS, 452, 2343APPENDIX A. MEASURING THE VERTICAL ACTION
The vertical action j z is defined as: j z ≡ π (cid:90) d z d v z = 12 π (cid:73) d z v z , (A1)where z and v z describe the position and velocity along the orbit, respectively. If the potential is nearly axisymmetricthis quadrature can be calculated using the algorithm described in (Solway et al. 2012). We adopted here a similarprescription as explained below: Vera-Ciro & D’Onghia Fig. A1.—
Surface of section for a star particle measured at two different times. In the left panel the simulation has already developedspiral arms, the right panel is the surface of section after 2 Gyr. In this specific example, R gui ( t ) ≈ R gui ( t ) indicates that the particle hasnot migrated significantly. However, the particle moves away from a resonant orbit, a situation that is easily identified with our scheme.
1. At a given time t the gravitational potential Φ, its radial gradient ∂ Φ /∂R , and vertical gradient ∂ Φ /∂z arecomputed on a polar grid ( R, θ, z ) centered at the minimum of the (fixed) halo potential. The grid dimensionis N R × N θ × N z = 128 × × (cid:15) Φ < R <
30 kpc(logarithmically spaced), − π < θ < π and − < z < (cid:15) Φ = 50 pc is the softening length. Notethat these ranges span six times the scale-length and twenty times the scale height. The potential is numericallycalculated at each node of the grid using a tree algorithm that includes up to the quadrupole contribution. Wedo not impose a symmetry about the galactic mid-plane, in fact, modest asymmetries that can develop over timeare naturally captured by the grid.2. At each point ( R, z ) the potential is calculated by averaging over the N θ = 128 azimuthal points:Φ( R, z ) ≡ N θ N θ (cid:88) j =1 Φ( R, θ j , z ) , (A2)Equivalent expressions are derived for ∂ Φ /∂R and ∂ Φ /∂z . This approximation relies on the fact that the potentialis nearly axisymmetric, this however may not be the case after if a bar grows at the center of the galaxy.Once the potential is computed at t time, for each star particle the following quantities are evaluated: the verticalangular momentum L z , the guiding radius R gui as the solution to the equation: L z = R gui V circ ( R gui ) , (A3)and the vertical action j z , with the procedure described below.First, the position and velocity vectors are transformed to cylindrical coordinates ( R, z ). This point is integrated inthe fixed potential defined by the grid for a time t max = 10 T z , where T z = 2 π/ν is the period of vertical oscillations,and ν is the vertical epicycle frequency (Binney & Tremaine 2008) ν = ∂ Φ( R, z ) ∂z (cid:12)(cid:12)(cid:12)(cid:12) z =0 . (A4)The vertical frequency can be easily calculated using five-point stencil derivatives applied to the ∂ Φ /∂z grid.The integration is performed using a symplectic integrator of order 4 (Candy & Rozmus 1991). With this choice,typical deviations in the energy are always smaller than 10 − and show no secular growth over the integration time t max . This procedure yields a set of discrete points along the orbit of the form { t i , R i , z i , v R,i , v z,i } i .During the integration of each star, every time the condition R i < R gui < R i +1 , v R,i > { R ( t ) , z ( t ) , v R ( t ) , v z ( t ) } to { t ( R ) , z ( R ) , v R ( R ) , v z ( R ) } which is numerically integrated in therange R i ≤ R ≤ R gui . With this change of coordinates, the symplectic structure clearly breaks, therefore we have touse a different (non-conservative) integrator. In this case we use the explicit embedded Runge-Kutta Prince-Dormand(8, 9) method (Dormand & Prince 1980) implemented in the Gsl library. After the exact intersection with the R = R gui plane is found, the integration is resumed from R i +1 with our symplectic method.In Fig. A1 we show an example for a random particle in our disk 1 simulation. The green points show the setobtained with the procedure described above. The blue line represents the zero-velocity curve, and gray points areother orbits, shown here to depict the structure of the phase-space in our runs. At a first glance, the set in the left n the conservation of the vertical action 9panel resembles a simple closed structure in the ( z, v z )-plane, however, a careful inspection reveals a more complexstructure. Indeed, the surface of section of this particle is composed of five resonant islands. At later times, however,the orbit traces a single closed structure (right panel). This transition is non-adiabatic and as such we would like todevise a way to identify these cases.The task of calculating the vertical action for a given star is then reduced to measuring the area enclosed by the setof points in the ( z, v z )-plane found above as the star moves in the galaxy. To achieve this, the points in the surface ofsection were connected with the following procedure:1. In the ( z, v z )-plane a random point p ≡ ( z , v z, ) is selected and its closest neighbor p is found. For simplicity,both coordinates ( z, v z ) are rescaled in the interval ( − , N points to p , { q k } Nk =1 that have not been already included in a previous set are then selected. Foreach point q k the angles between the segments p p and p q k are defined.3. The point with minimum angle weighted by the distance is found and defined as the the new p . We noted thatin our implementation the weights ∼ | p q k | perform reasonably well.4. Steps 2 and 3 are repeated until p is the next point to be included in the cycle.5. Step 1 is then repeated until all points have been included in our analysis.This algorithm then yields a set of unconnected loops in the surface of section defined by R = R gui and v R >
0. Inthe cases where only one loop is present the action is just the area of such cycle and can be calculated using the popularformula based on the Green theorem. Otherwise, the particle is trapped in a resonance and it is no longer consideredin our results. This implementation uses N = 16 and rejects any cycle that has fewer points than this threshold. Inaddition, orbits for which the distance between the last and first point is larger than five times the average distancebetween the other points in the cycle are rejected. This choice avoids cases in which the algorithm ends before closinga loop, a circumstance that occurs for stars whose orbit takes them close to the zero-velocity curve. If the orbit iscomposed of several loops (e.g. left panel in Fig. 3) and the procedure fails to calculate the area for one of those loops,then the orbit is rejected altogether.Fig. A1 shows the results of the implementation of our algorithm to the particles of our simulations, as illustratedby the red line of the inset of the left panel. In our numerical experiments the procedure successfully computes thevertical action for ∼
85% of the particles located at a distance of R (cid:38) ∼∼