{'fraction_non_alphanumeric': 0.03140368548912721, 'fraction_numerical': 0.01632399123066896, 'mean_word_length': 4.393034030995366, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the conservation of the vertical action on galactic disks

We 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 stars in an axisymmetric potential. In the multi-armed galaxy the low-amplitude arms represent tiny perturbations of the potential, hence the vertical action for a set of stars is conserved, although after several orbital periods of revolution the conservation degrades significantly. For a Milky Way-like galaxy with vigorous spiral activity and the formation of a bar, our results show that the potential is far from steady, implying that the action is not a constant of motion. Furthermore, because of the presence of high-amplitude arms and the bar, considerable in-plane and vertical heating occurs that forces stars to deviate from near-circular orbits, reducing the degree at which the actions are conserved for 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 forms by radial migration of stars, under the assumption of the conservation of the action describing the vertical motion of stars.


INTRODUCTION
A stellar system can be fully modeled by a phase-space density distribution function (DF), whose evolution is described by solutions of the Boltzmann equation (Binney & Tremaine 2008). Is also a well-known result that for steady-state potentials, the DF describing the galaxy can be represented as a function of the integrals of motion of the system (Jeans 1919). There are several choices one can make to use then as coordinates of the DF, e.g. energy and angular momentum (Eddington 1915). Integrated orbits in axisymmetric potentials indicate the existence of a non-classical isolating integral of motion (Ollongren 1962), which is closely related to the actions of the system (Binney & Spergel 1984).
Besides being a natural choice to describe the evolution in the phase-space of collisionless stellar systems, the actions are also adiabatic invariants. This is a very appealing property of spiral galaxies, where the high frequency of the vertical motion, compared to in-plane evolution, makes the action describing the z motion a prime candidate for being a conserved quantity, even in the presence of spiral arms. This assumption has been tested by Solway et al. (2012), who showed that for isolated galaxies the vertical action is not a constant of motion of individual stars, and is only conserved in average for samples of stars, with an intrinsic dispersion of ∼ 20%.
It is known that recurring spiral arm-activity in a galactic disk causes stars to diffuse radially over time in a manner that largely preserves the overall structure of the disk. Sellwood & Binney (2002) used bi-dimensional stellar disks to show that stars near the corotation radius of a transient spiral pattern change angular momentum and move to new radii inwards or outwards with negligi-ble increase of random motion, causing a process termed radial migration. Following studies confirmed that this process occurs in 3D simulations of galactic disks either in isolation or in a cosmological context (Roškar et al. 2008;Minchev et al. 2011;Grand et al. 2012;Loebman et al. 2011;Kubryk et al. 2014;Halle et al. 2015;Grand et al. 2015).
If radial mixing in galactic disks occurs, this process has interesting implications for the spread in the agemetallicity relation observed in the solar neighborhood (Edvardsson et al. 1993). If stellar migration never occurred, the stars in the solar vicinity would have metallicities ranging from zero to the present-day value (Serenelli et al. 2009) with a significant correlation between the age of the stars and their metallicities. However this correlation has never been found in the current data, raising the possibility that the stars migrated from their birth radii (Wielen et al. 1996). Because the radii of stars at the birth depends on their specific angular momentum, stellar migration occurs if there is a change of the angular momentum of the stars. Changes in the angular momentum for individual stars can be induced by scattering with giant molecular clouds (Spitzer & Schwarzschild 1953;Mihalas & Binney 1981) or transient spiral structures (Jenkins & Binney 1990;Jenkins 1992), interference of coherent long-lived modes and density waves with a bar (Roškar et al. 2012), or by interactions with satellite galaxies (Kazantzidis et al. 2008;Villalobos & Helmi 2008;Bird et al. 2012).
In this paper we aim to study in detail the effect of migration on isolated disks, by using N -body simulations for disk galaxies evolving in isolation, to achieve high-resolution; with the aim of better determining the contribution of the spiral waves to the heating of the stellar disk, to the radial migration of stars in the disk, and ultimately to the formation of a thick disk.
This manuscript is laid out as follows. In Section 2 we introduce the N -body simulations used in our work and describe the method to measure the vertical action. In Section 3 we discuss to what extent this quantity is conserved during the evolution of the disk, and in Section 4 we show how this can be used to study the formation of a thick disk. Finally in Section 5 we draw our conclusions.

Simulations
For this study we employed N -body simulations of two galaxy models labeled as disk 1 and disk 2. Each galaxy consisted of a static dark matter halo and a live rotationally supported stellar disk of 5 × 10 6 particles. Because the spiral morphology in stellar disks depends in detail on the mass distribution of the galaxy and on the self-gravity of the disk, we selected the structural parameters of the galaxies such that one galaxy developed a low-amplitude multi-armed disk (disk 1), whereas the other developed a spiral morphology more similar to a Milky Way (MW)-like galaxy (disk 2). This was achieved by choosing models with very different disk mass fractions within 2.2 scale-lengths: M disk /M total = 0.271 for disk 1 and 0.478 for disk 2, respectively. By varying the mass distribution of the disk within 2.2 scalelengths, these models were designed to change the critical length scale parameter λ crit = 4π 2 GΣ/κ 2 , where Σ is the disk 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 component were assumed to be independent and both models were constructed similarly to the approach described in (Hernquist 1993;Springel et al. 2005). The dark matter mass distribution was modeled assuming the Hernquist profile (Hernquist 1990): with M halo being the total halo mass and R halo the scale radius. The effect of the dark halo is modeled assuming a static potential, this is of course a limitation of our model since the system never completely relaxes, however, the complication of adding an alive component comes with the price of including further scattering to the particles of the disk. The stellar disk is modeled in the initial conditions as a thin exponential surface density profile of scale-length R disk and total mass M disk . The vertical mass distribution of the stars in the disk is specified by the profile of an isothermal sheet with a radially constant scale height z disk : Model disk 2 included a static potential for the bulge described by the Hernquist model with mass M bulge and scale-length R bulge . We aim to build this galaxy model with a flat rotation curve consistent with the observations of the MW within two optical radii (R 14 kpc) (Reid et al. 2014;Bovy et al. 2009). We noted that this requirement is satisfied by adopting a Hernquist profile for the dark halo with a large scale radius, and a large total mass. However, the dynamical range of interest of this study is a few optical radii, thus the mass distribution outside this radius does not affect our analysis. Table 1 summarizes the structural parameters adopted for each galaxy model in this study. Fig. 1 shows the azimuthally averaged radial profiles for the total circular velocity V circ (stars and dark matter), the Toomre parameter Q Toomre of the stellar disk, the disk surface density Σ and its vertical velocity dispersion σ z . The profiles are displayed at four different times, spanning 6 Gyr of evolution. In what follows we use the notation t n to refer to the time of integration in our simulations in Gyr, with t 0 thus labeling the the initial time, and t 6 corresponding to 6 Gyr.
We note that in the multi-armed galaxy (disk 1 model), Q Toomre increases with time, as a consequence of the increase of the radial velocity dispersion. This is expected, since in-plane heating occurs as a consequence of the spiral structure activity, while the vertical structure remains unchanged. This is, however, not the case for our disk 2 run. A bar instability grows at the center of the galaxy after ∼ 3 Gyr. Despite the flat rotation curve, the surface mass density of the disk changes considerably with time and the vertical velocity dispersion increases by ∼ 15% in the shown radial range.
In Fig. 2 a sequence in time of the face-on views of the stellar disk of both models of galaxies are displayed. The colored bar indicates the values for the surface density of the disk. As expected in the context of swing amplification theory, the two models show very different spiral morphologies. Our disk 1 model grows a large number of low-amplitude spiral arms that self-perpetuate for at least 6 Gyr (D'Onghia et al. 2013). After 1 Gyr disk 2 already developed four high-amplitude spiral arms at two scale-lengths (∼ 5 kpc), in remarkable agreement with the number of arms predicted for a disk with 50% disk fraction (D'Onghia 2015; Pettitt et al. 2015). Later in time this disk model develops a bar (Merritt & Sellwood 1994;Martinez-Valpuesta et al. 2006;Athanassoula & Martínez-Valpuesta 2008;Yurin & Springel 2015), which becomes the dominant non-axisymmetric feature after ∼ 4 Gyr of evolution.
Note that the efficiency of the stellar radial migration by recurrent spiral activity in the multi-armed galaxy (disk 1 model) has been recently studied by Vera-Ciro et al. (2014). As first shown in Sellwood & Binney (2002), the authors confirmed that spiral activity in a low-mass disc causes stars to diffuse radially, with stars near corotation of a spiral pattern exchanging angular momentum and moving to a new radius without adding random motion (see also Daniel & Wyse 2015). Indeed, the mech- -Radial profiles of total circular velocity, V cir , Toomre parameter, Q Toomre , surface mass density, Σ, and vertical velocity dispersion, σz, for the two studied galaxy models: disk 1 (left) and disk 2 (right). tn refers to the time of integration in our simulations in Gyr, with t 0 being the the initial time and t 6 corresponding to 6 Gyr.
anism of scattering of stars at corotation predicts a lack of significant disk heating and leaves the surface density profile unchanged. These two features are reproduced in our disk 1 model as shown in Fig. 2 (left panels). In the MW-like galaxy model (disk 2) the bar formation causes greater angular momentum changes than those of a multi-armed spiral structure (Minchev & Famaey 2010). Furthermore, because the bar persists after its formation, the associated angular momentum changes might differ from those caused by recursive spiral arms. Therefore the surface density profile of the disk is expected to change. Some disk heating is also expected to occur with time.

Measuring the vertical action
The action-angle variables can be used to describe the evolution of orbits in static potentials. The actions, j, are constants of motion in the unperturbed field, or when the potential changes slowly so that they are adiabatic invariants (Arnold 1978;Goldstein et al. 2002). The vertical action is defined as: where z and v z describe the position and velocity along the orbit, respectively. There are a number of ways to estimate of the vertical action of a star particle in a The disk 1 simulation develops a large number of low-amplitude spiral structure that perpetuate for at least 6 Gyr (left panels); disk 2 develops approximately a four-fold rotational symmetry and becomes bar-unstable after ∼ 3 Gyr. nearly axisymmetric potential. Solway et al. (2012), for instance, calculated the quadrature in Eq. (3) by measuring the area enclosed by the orbit in the z − v z plane. We adopted here the prescription as described in Appendix A. This method allows us to identify and reject stars trapped in resonances, a case where the adiabatic invariance breaks (Pfenniger 1984;Pfenniger & Friedli 1991).
In Fig. 3 we show an example for a random particle in our disk 1 simulation at two different times. The green points show the surface of section of the integrated orbit and the red line is the connected loop used to calculate Eq. (3). The blue line represents the zero-velocity curve, and gray points are other orbits, shown here to depict the structure of the phase-space in our runs. The area defined by the red-loop in Fig. 3 is the vertical action j z , for this particular example, the vertical action measured at t 1 is a factor of ∼ 80 larger than it is at time t 3 , this is clear evidence that the value of j z is not conserved for individual orbits for slow-varying, low-amplitude perturbations, such as the ones developed in our disk 1 run.

A CONSERVED QUANTITY?
Our previous calculations demonstrated that the vertical action j z is not a conserved quantity for individ- Fig. 3.-Surface of section for a star particle measured at two different times. In the left panel the simulation has already developed spiral arms, the right panel is the surface of section after 2 Gyr. In this specific example, R gui (t 1 ) ≈ R gui (t 3 ) indicates that the particle has not migrated significantly. However, its vertical action has changed by a factor ∼ 80 (δjz = 4.4). For each case the solid red line displays the closed path used to define the action and the blue line indicates the zero-velocity curve. ual stars. This result agrees with findings of previous studies (Solway et al. 2012). We investigated this property in our disk 1 and disk 2 simulations. After 1 Gyr of evolution an ensemble of stars with guiding radius: R gui = 5.00 ± 0.25 kpc was selected and the relative changes in the vertical action were computed at a later time (See Appendix A). This radius is close to two disk scale-lengths in both galaxy models. The vertical action changed significantly over a period of time as shown in Fig. 3. Thus, its change from time t to a later time t > t was quantified in the following way: In addition, we noted that over the same period of time a significant radial migration occurs, with a consequent change of the guiding radii. Thus, a similar definition was also applied to δR gui . Fig. 4 (top panel) displays the outcome of our study for the multi-armed galaxy disk 1 over a period of 2 Gyr. In agreement with previous studies, our findings suggest the vertical action is a conserved quantity on average for an ensemble of stars. Note that the modest fractional changes of δj z in our simulation match the values reported by Solway et al. (2012) 3 Furthermore, Fig. 4 indicates that any variation in the vertical action, δj z , is independent of the radial migration, measured by changes in guiding radii, δR gui . Next, we computed the median of the fractional change in the vertical action µ(δj z ) and the 1σ equivalent dispersion around the median σ(δj z ). We obtained µ(δj z ) = 0, with standard deviation σ(δj z ) = 0.28 (displayed as the solid black line in Fig. 4), also in agreement with the values reported in Solway et al. (2012).
The same analysis performed on the more massive disk 2 produces different results as shown in Fig. 4 (bottom panel). The median of fractional changes in the vertical action is µ(δj z ) = 0.09 with a dispersion σ(δj z ) = 0.64, significantly larger than the case of the multi-armed galaxy. Thus, changes of angular momentum occurred with the formation of the bar, with significant disk heating, combined with the strong variation of the potential of the disk, do not conserve the vertical action at the same level of accuracy.
To demonstrate this unexpected outcome, the evolution by time of the median µ(δj z ) and its dispersion σ(δj z ) of fractional changes in the vertical action are displayed in Fig. 5 for the two galaxy models. Ensembles of stars are selected at after 1 Gyr of evolution (t = t 1 ) with guiding radii R gui = 3, 5, 8 kpc and followed over the evolution time of the disk. While in the multi-armed disk, even after 40 orbital periods the average variation in the vertical action is never larger than 1%, in the MW-like disk galaxy the median reaches values of 15-20% after 20 orbital periods.
Note that in Fig. 5 the dispersion around the median steadily increases as a function of time for both models (left panels), suggesting that the vertical action of an ensemble of stars tends to be less conserved as time progresses. This results holds for both the multi-armed disk and the MW-like disk.
We also checked whether the vertical heating of the disk, measured by the vertical velocity dispersion of the ensemble of stars, could produce changes in the average vertical action. Fig. 5 plots the median µ(δj z ) and dispersion σ(δj z ) of fractional changes in the vertical action against the corresponding variation in the velocity dispersion of the vertical component.
To test for the origin of the changes in the average action we note that the vertical energy of samples of stars proportionally depends on the vertical velocity dispersion j z ∼ σ 2 z (cf. Eq. (7)), and so if the average vertical action increases as a consequence of the vertical heating in the disk, then δj z = 2δσ z . The multiplying factor in this approximation is the vertical epicycle frequency, which remains mostly constant in time (See Fig. 1). Fig. 5 shows that this is not the case, and that the increase in dispersion of the changes in the vertical action is not merely determined by the vertical heating of the disk.

IMPLICATIONS FOR THE VERTICAL STRUCTURE
What does the conservation of the vertical action imply for the structure and evolution of the stellar disk? It has been argued that under the assumption of the conservation of the vertical action j z stars migrating outwards in the disk tend to move far from the mid-plane and populate the thick disk (Roškar et al. 2008;Schönrich & Binney 2009). We address this question here.
The epicycle approximation regulates the motion of stars on near-circular orbits. In the approximation the vertical motion of a star is also decoupled from the inplane motion. The motion perpendicular to the plane is then described by the Hamiltonian (Binney & Tremaine 2008): For an ensemble of star particles at radius R, the average vertical energy is where σ z is the vertical velocity dispersion of the ensemble. Note that σ z in Eq. (6) does not necessarily refer to the value of velocity dispersion of all stars with home radii at R. Indeed, Eq. (6) holds also for any population of stars whose guiding radii have changed over time (e.g. stars migrating from inner or outer locations within the disk) and that are used to calculate it. Thus, in the epicycle approximation, the vertical action is the following j z = H z /ν (Binney & Tremaine 2008): where, α 2 h 2 = z 2 , z 2 is the second moment of the vertical component, and α = π/ √ 12. With this normalization, the parameter h is equivalent to the disk scale height z disk , h ≡ z disk , when the second moment is calculated for all stars sampled from the sech 2 z/z disk model used for the vertical density of our simulations (cf. Eq. (2)).
Once the average vertical action and vertical velocity dispersion 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 different times for both galaxy models (top panels). The solid lines show the average vertical action calculated for a set of stars using the algorithm described in Section 2. The dashed lines illustrate the predicted vertical action using the epicycle approximation for that set of stars. Clearly, the epicycle approximation predicts values of the vertical action 15% higher than the estimates obtained with the numerical procedure. A similar trend is observed in the numerical experiments described in Solway et al. (2012). The scale height profile of the whole disk at t 3 is shown by the solid blue line. In green we show the thickness profile of the migrators at t 1 and in red the thickness at t 3 . The orange dashed line is the result of using Eq. (7) and assuming that the vertical action is conserved j z,epi (t 1 ) = j z,epi (t 3 ) .
Then, Eq. (7) is used to compute the thickness of the disk h at each radius, and the values obtained for the both galaxy models are shown in Fig. 6 (bottom panels). A similar discrepancy is reflected in the estimate of the scale height at different radii. The values derived in the epicycle approximation are 5% larger than the values obtained by evaluating Eq. (3) as previously described.
In both simulations, after the disks evolved for 1 Gyr (t 1 ), we selected a set of stars with guiding radii R gui = 5.00 ± 0.25 kpc (left panels in Fig. 7) and R gui = 8.00 ± 0.40 kpc (right panels), we follow these stars for two billion years longer (t 3 ) and calculated their scale height h profile as a function of their final position in the disk. This would show if a subsample of migrating stars can change their thickness h as they migrate.
We only select stars that are in nearly circular orbits L z /L circ > 0.97, where L z is the vertical angular momentum and L circ is the angular momentum that a star has, should it move in a circular orbit with the same energy. By adding this constraint at both times t 1 and t 3 we ensure that the sample contains stars that migrate while preserving their circular orbits (Sellwood & Binney 2002).
For this set of particles the average vertical action in the epicycle approximation is also j z is computed at both times t 1 and t 3 . The thickness at both times is also included, in general stars that move toward the center tend to reduce their scale height, whereas stars that migrate outwards increase it. For comparison we also show in blue the scale height of the disk as a function of radius, calculated with the normalized RMS value of the z-coordinates of all particles in the disk. This value matches the one reported as z disk in table 1 (blue solid line).
The bias reported in (Vera-Ciro et al. 2014) is observed for both disks. This explains the fact that for stars that migrate outwards, the slight increase in their vertical scale height is not enough to create a thicker component at their final location. We also included in Fig. 7 the predicted thickness h at final t 3 under the assumption that the action is conserved. It can clearly be seen that for our disk 1 no thick distribution of stars is either observed or predicted even when the average vertical action is clearly conserved.
In the MW-like galaxy, disk 2, a bar grows in the inner parts of the disk, in qualitative agreement with simulations of Brunetti et al. (2011). As a consequence of the bar formation and the vigorous spiral activity there are not only significant changes of angular momentum, but also a considerable in-plane heating that forces stars in this simulation to deviate from a near-circular orbits. There are good reasons to believe that in this case the energy of vertical motion is clearly not decoupled from that of the radial part of the motion. Indeed, the epicycle approximation is a poor description of the motion of most particles, for which the radial and vertical oscillations are neither harmonic nor are the energies of the two oscillations decoupled. Furthermore, the actions are constants of motion under the conditions of slow changes of the potential. In the multi-armed galaxy the low-amplitude arms represent tiny perturbations of the potential, hence the vertical action might still be conserved, although after several orbital periods of revolution the vertical action is not conserved at the same level of accuracy. However, in the MW-like galaxy the high-amplitude spiral structure and the bar formation seem to compromise the validity of these conditions.

SUMMARY AND CONCLUSIONS
We have presented quantitative estimates of the vertical action of ensembles of stars in the presence of spiral activity. The vertical action of star particles has been evaluated in isolated disks of galaxies with different spiral morphologies. The two galaxy models consist of a low-mass multi-armed disk (disk 1) and a MW-like disk that develops a bar after ∼ 3.5 Gyr of evolution (disk 2).
For the multi-armed disk galaxy, disk 1, we find that the vertical action for an ensemble of stars is approximately a conserved quantity, with a dispersion of fractional changes around its median of 20%, in agreement with the results reported in the literature. However, as time progresses the action is conserved with less accuracy, as indicated by the increase of the standard deviation up to the value of 50% after 6 Gyr of evolution. This also holds for a subset of stars that radially migrates by resonant scattering at corotation. These results have implications for the vertical structure of the disk when stars radially migrate to the outer part of the disk. In this model stars migrating outwards are a heavily biased subset of stars with preferentially low vertical velocity dispersion. Thus, extreme migrators outwards are characterized by having a have small amplitude of their vertical excursion, independent of the conservation of the vertical action.
In the MW-like galaxy model, disk 2, the bar formation causes large angular momentum changes, with significant in-plane and vertical heating. Thus, in the presence of the formation of the bar and combined with significant disk heating, our results show that the vertical action cannot be assumed a constant of motion. Indeed, in this case the fractional changes in the estimates of the action for a set of stars progressively increase by time, with a standard deviation that approaches 80% after 6 Gyr of disk evolution. If confirmed, this result casts doubts on the possibility of predicting the vertical structure 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 and ATP NASA Grant No NNX144AP53G. ED gratefully acknowledges the support of the Alfred P. Sloan Foundation. We express our appreciation toward the Aspen Center for Physics for their hospitality, funded by the NSF under Grant No. PHYS-1066293. Simulations have been run on the High Performance Computing cluster provided by the Advanced Computing Infrastructure (ACI) and Center for High Throughput Computing (CHTC) at the University of Wisconsin. We thank Daniel Pfenniger, Jerry Sellwood, Martin Weinberg, Victor Debattista, Rob Grand and Jonathan Bird for the very helpful comments and constructive discussions. -Surface of section for a star particle measured at two different times. In the left panel the simulation has already developed spiral arms, the right panel is the surface of section after 2 Gyr. In this specific example, R gui (t 1 ) ≈ R gui (t 3 ) indicates that the particle has not 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 are computed on a polar grid (R, θ, z) centered at the minimum of the (fixed) halo potential. The grid dimension is N R × N θ × N z = 128 × 128 × 128, with the cylindrical coordinates ranging as follows: Φ < R < 30 kpc (logarithmically spaced), −π < θ < π and −6 kpc < z < 6 kpc, where Φ = 50 pc is the softening length. Note that these ranges span six times the scale-length and twenty times the scale height. The potential is numerically calculated at each node of the grid using a tree algorithm that includes up to the quadrupole contribution. We do not impose a symmetry about the galactic mid-plane, in fact, modest asymmetries that can develop over time are naturally captured by the grid.
2. At each point (R, z) the potential is calculated by averaging over the N θ = 128 azimuthal points: Equivalent expressions are derived for ∂Φ/∂R and ∂Φ/∂z. This approximation relies on the fact that the potential is 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 vertical angular momentum L z , the guiding radius R gui as the solution to the equation: 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 in the fixed potential defined by the grid for a time t max = 10 3 T z , where T z = 2π/ν is the period of vertical oscillations, and ν is the vertical epicycle frequency (Binney & Tremaine 2008) 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 −6 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 During the integration of each star, every time the condition R i < R gui < R i+1 , v R,i > 0 is satisfied, the dynamical states are changed from With this change of coordinates, the symplectic structure clearly breaks, therefore we have to use 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 4 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 set obtained with the procedure described above. The blue line represents the zero-velocity curve, and gray points are other orbits, shown here to depict the structure of the phase-space in our runs. At a first glance, the set in the left panel resembles a simple closed structure in the (z, v z )-plane, however, a careful inspection reveals a more complex structure. 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 to devise 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 set of points in the (z, v z )-plane found above as the star moves in the galaxy. To achieve this, the points in the surface of section were connected with the following procedure: 1. In the (z, v z )-plane a random point p 0 ≡ (z 0 , v z,0 ) is selected and its closest neighbor p 1 is found. For simplicity, both coordinates (z, v z ) are rescaled in the interval (−1, 1).
2. The closest N points to p 1 , {q k } N k=1 that have not been already included in a previous set are then selected. For each point q k the angles between the segments p 0 p 1 and p 1 q k are defined.
3. The point with minimum angle weighted by the distance is found and defined as the the new p 1 . We noted that in our implementation the weights ∼ |p 1 q k | 3 perform reasonably well.
4. Steps 2 and 3 are repeated until p 0 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. In the cases where only one loop is present the action is just the area of such cycle and can be calculated using the popular formula based on the Green theorem. Otherwise, the particle is trapped in a resonance and it is no longer considered in our results. This implementation uses N = 16 and rejects any cycle that has fewer points than this threshold. In addition, orbits for which the distance between the last and first point is larger than five times the average distance between the other points in the cycle are rejected. This choice avoids cases in which the algorithm ends before closing a loop, a circumstance that occurs for stars whose orbit takes them close to the zero-velocity curve. If the orbit is composed 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 illustrated by the red line of the inset of the left panel. In our numerical experiments the procedure successfully computes the vertical action for ∼ 85% of the particles located at a distance of R 2 kpc from the disk center, and ∼ 60% of particles orbiting at smaller radii.