The Age-Dependent Vertical Actions of Young Stars in the Galaxy

Disk galaxies, including the Milky Way, are thought to build their stellar mass mostly via star formation, i.e. either from their own gas or from gas brought in by mergers (e.g. Tinsley & Larson, 1978; Fall & Efstathiou, 1980; Steinmetz & Muller, 1995). The dynamical evolution and spatially varying star formation histories shape these galaxies after initial mass accumulation. At the late stages of galaxy evolution, stars are thought to form on nearly circular orbits from thin disks of gas. The time-varying and non-axisymmetric structures present within disks can subsequently heat stellar orbits by increasing their radial actions (or eccentricities) and their vertical actions (or random motions in the vertical direction) (Spitzer & Schwarzschild, 1953). These processes have been observed in simulations (e.g. Sellwood, 2014; Jenkins & Binney, 1990) and have evidence in external galaxies (e.g. Courtois et al., 2015; Vasiliev et al., 2022; Vera-Ciro & D’Onghia, 2016).

Mechanisms like minor mergers and external perturbations can abruptly increase stars’ vertical motions, whereas others, such as interactions with giant molecular clouds (GMCs) or other massive structures, contribute to a more gradual heating effect (Qu et al., 2011; Spitzer & Schwarzschild, 1951; Lacey, 1984; Carlberg, 1987; Jenkins & Binney, 1990). In addition, vertical breathing waves induced by the bar and spiral arms (Barbanis & Woltjer, 1967; Masset & Tagger, 1997; Monari et al., 2015; Grand et al., 2016) may cause heating.

In the Milky Way, star-by-star astrometric and spectroscopic data open a window to a detailed picture of disk galaxy evolution. This data reveals that both stars and the molecular gas from which they form exhibit low velocity dispersion initially (e.g. Bovy et al., 2012; Li et al., 2013). In the solar neighborhood, older stars present higher vertical velocity dispersion ($\sigma_{z}$) than younger stars (e.g. Hayden et al., 2020; Yu & Liu, 2018). As mentioned before, different dynamical processes can induce heating, i.e., acquisition of vertical motion by stars over time. The secular heating of stars has been studied through age-velocity dispersion relations (Wielen, 1977; Quillen & Garnett, 2001; Aumer & Binney, 2009; Aumer et al., 2016; Sanders & Das, 2018), suggesting that the present-day vertical-velocity dispersion continuously rises and scales as $\sigma_{z}\sim\tau^{0.5}$, where $\tau$ represents stellar ages (Nordström et al., 2004; Holmberg et al., 2009). The present-day distribution of vertical motions of stars has also been modeled as a combination of their birth velocity dispersion, followed by a heating period that includes radial dependence (Ting & Rix, 2019).

In recent years, the Gaia collaboration has released the most extensive census of positions, velocities, and other stellar parameters of billions of stars in the Milky Way, including luminous, massive, and hot stars.

These stars are massive, young, and can be used to study the spiral structure and the dynamics of the young Milky Way disk (Xu et al., 2018; Zari et al., 2021; Poggio et al., 2021).

Stars’ orbits at birth can yield insights into the dynamic state of the gas from which they formed. Subsequent vertical evolution may provide details regarding the heating mechanisms in the disk. However, the characterization of stellar heating on short timescales remains an open question. Addressing these questions involves studying young stars in the Milky Way’s disk. The relationship between the vertical velocity dispersion ($\sigma_{z}$) and stellar age in the galaxy has been a subject of debate primarily due to the challenges in determining stellar ages and the varying selection criteria across different surveys. However, precise parallax data from the Gaia mission (Gaia Collaboration et al., 2016, 2021), combined with robust spectroscopic stellar parameters delivered from large spectroscopic surveys, yields uniform age estimates precise to the $\sim 10$% level for a large number of stars (e.g. Xiang & Rix, 2022), enabling new insights into the vertical evolution history of the Milky Way. Understanding the age-vertical action relation is a crucial step in understanding the evolution of the galaxy.

In this work, we use the positions and velocities of young OBA stars from Gaia EDR3 as presented by Zari et al. (2021) to map the present-day $J_{z}$ as a function of stellar age $\tau$ and the current Galactocentric radius. This approach offers new insights into the dynamical state of the disk. We use $J_{z}$ because $\sigma_{z}$ can change due to both gradual variations in the mid-plane baryon density and radial migration. Therefore, to better characterize the global changes in stellar orbits, it is preferable to describe the vertical motion of stars by their vertical actions and to understand heating as an increase in $J_{z}$ (Ting & Rix, 2019).

The paper is organized as follows. Section 2 introduces the OBA star sample and its age distribution. Section 3 describes the theoretical model for the vertical action with respect to galactocentric radii and the subsequent Bayesian inference. In Section 4, we present the posteriors of the model and we expand the results that show a linear behavior of young stars’ vertical action with age. In section 5 we briefly discuss the robustness of the results, data selection, and dust extinction and possible physical explanation.

Our objective is to investigate the vertical actions of stars at their birth within the Milky Way disk. Therefore, we require a sample with estimable ages of young stars, specifically, those younger than a few dynamical periods ($\lesssim 0.5$ Gyr). Additionally, this sample should cover a substantial portion of the Galactic disk and possess complete 6D phase space information, which is essential for accurately estimating the orbital actions.

We constructed our sample based on the catalog of OBA star candidates (Z21) presented by Zari et al. (2021). They used a combination of Gaia EDR3 photometry and astrometry, along with 2MASS photometry, to select candidate OBA stars and exclude contamination from red giant branch and asymptotic giant branch stars. For their filtered sample, they further selected stars in proximity to the mid-plane ($|z|<300$ pc) and belonging to a low-velocity dispersion component in the vertical direction (see sections 2 and 4 of Zari et al., 2021).

We then cross-matched our data set with the LAMOST catalog presented in Xiang et al. (2022) to obtain stellar parameters, including the radial velocity $v_{\mathrm{los}}$ (see Figures 1 and 2). The values for $v_{\mathrm{los}}$ in Xiang et al. (2022) are derived from the redshift estimates of the LAMOST spectra. However, the realistic $v_{\mathrm{los}}$ uncertainty is primarily due to wavelength calibration error, given the low spectral resolution of LAMOST. As a result, the uncertainty of radial velocity reported in Xiang et al. (2022) is likely underestimated. In their comparison of velocity measurements across multiple visits for common stars, Xiang et al. (2022) found a typical scatter in $v_{\mathrm{los}}$ of about 10 km/s. We adopt this value for the noise model of $v_{\mathrm{los}}$ for all stars in our study. The observed scatter from repeated measurements also includes contributions from binary stars. Since binary stars are not directly modeled in this work, we account for their effects on the velocity by marginalizing over the observed broadening of the velocity dispersion attributable to binary star populations.

We estimate ages by employing a Bayesian isochrone fitting procedure as presented in Xiang & Rix (2022). The inputs for our method include parallax, the LAMOST spectroscopic stellar parameters ($T_{\mathrm{eff}}$, $\log g$, [Fe/H]), and the multi-band photometry from Gaia ($G$, $BP$, $RP$) and 2MASS ($J$, $H$, $K$). For age inference, we adopt MIST stellar isochrones (Dotter, 2016; Choi et al., 2016). To ensure the high precision of the spectroscopic stellar parameters, and thus the accuracy of the age estimates, stars with a LAMOST spectral S/N lower than 30 are excluded from our sample.

Using the astrometric parameters ($l,b,\mu_{l},\mu_{b}$), distances, and radial velocities, we compute the six-dimensional phase-space coordinates in the cylindrical Galactocentric reference frame ($R,\phi,Z,V_{R},V_{\phi},V_{Z}$) and the subsequent vertical actions. We adopt a left-handed coordinate system for the frame transformations, in which the x-axis is directed towards the Galactic center, and the y-axis increases to the right when viewed from the North Galactic Pole, opposite the direction of Galactic rotation. For the frame transformations, we adopt a vertical distance of the Sun above the plane of $Z_{\odot}=20.8$ pc (Bennett & Bovy, 2019), a distance of the Sun to the Galactic center of $R_{\odot}=8.15$ kpc (Reid et al., 2019), and a circular velocity at the Sun’s radius of $V_{c}(R_{\odot})=236$$\,\mathrm{km\,s}^{-1}$ (Reid et al., 2019). We assume a peculiar velocity of the Sun with respect to the local standard of rest of $(U_{\odot},V_{\odot},W_{\odot})=(11.1,11,7.25)$ km s^-1 (Schönrich et al., 2010; Reid et al., 2019).

For distances smaller than $\sim 4$ kpc, the dominant contributors to the fractional error distribution are parallax uncertainties. The distance fractional errors are below $10\%$ for $d<5$ kpc, and below $20\%$ for $d_{kin}<10$ kpc. The distances were estimated using a model designed to reproduce the properties of the data set in terms of spatial and luminosity distribution. $d_{kin}$ refers to the ‘astro-kinematic’ distances calculated by combining parallaxes and proper motions with a model for the expected velocity and density distribution of young stars as explained in Zari et al. (2021).

As previously mentioned, our focus is on young stars. Consequently, we have restricted the sample to stars younger than 500 Myr. Beyond 13 kpc, the reliability of our data significantly diminishes. Therefore, we have focused on stars within a Galactocentric radius range of 8 kpc to 13 kpc, which offers a sufficiently large sample size. To ensure data quality, we performed several filtering steps. We removed outliers with exceptionally large radial velocities ($>300,\mathrm{km/s}$), large errors in radial velocity ($>60,\mathrm{km/s}$), and high vertical actions ($J_{z}>20,\mathrm{kpc,km/s}$), leaving about $\sim 7900$ OBA stars. The resulting filtered sample will be referred to as Z21-500 in subsequent sections (i.e., a clean version of the sample of Zari et al. (2021) with ages less than 500 Myr).

The computed phase-space coordinates enable us to visualize the spatial distribution of our star sample. In Figure 1, we show in the left panel the spatial distribution of stars in our sample Z21-500 in the $R-Z$ plane and in the right panel the spatial distribution of stars in our sample in the $X-Y$ plane. The latter was calculated using a kernel density estimator (KDE) with a Gaussian kernel and a bandwidth of $0.1$ kpc. The observed low-density gap at the center of the map (around $x,y\approx 8,0$ kpc) may stem from the magnitude limit of LAMOST observations at the bright end, the intrinsically low population of OBA stars in the solar neighborhood, or the Sun’s location in the Local Bubble (see Zari et al., 2021). Without a comprehensive understanding of the joint selection function of the surveys, it is challenging to ascertain whether the observed overdensities and gaps represent artificial features or actual physical structures. The pencil-beam overdensities emanating from the Sun’s position are indicative of the LAMOST survey’s footprint.

To calculate the vertical actions $J_{z}$:

and radial actions $J_{R}$ of the stars in the sample, we use the axisymmetric Milky Way potential MWPotential2014 (Bovy, 2015), scaled to have a circular velocity of 236 km$,\mathrm{s}^{-1}$, and the Galpy package’s implementation of the Staeckel Fudge (Binney, 2012). We binned the data by Galactocentric radius, using a bin size of 1 kpc and spanning a range from $R_{\mathrm{GC}}=8$ to $13$ kpc. The distribution of $R$, ages, $J_{R}$, and $J_{z}$ is shown in Figure 2. Most stars are centered close to $10$ kpc, with vertical and radial actions close to $0$. Figure 3 displays the vertical action as a function of Galactocentric radii, colored by different ages. $J_{z}$ increases with age, suggesting that young stars are born on cold orbits, i.e., with low vertical actions, and undergo subsequent heating after birth. Motivated by these findings, we aim to analyze the initial actions of stars at birth and their heating over time.

In this section, we construct a physically motivated model to describe the vertical action at birth as a function of the Galactocentric radius and the subsequent temporal evolution of $J_{z}$ as a function of stellar age and radii¹¹1We could have equally well chosen to work with angular momentum rather than Galactocentric radius, but they are almost equivalent for young stars on near-circular orbits..

The model is then quantitatively compared to the data, taking into account observational uncertainties.

In this section, we summarize the results of our Bayesian parameter inference. Figure 5 presents a detailed overview of the posterior distributions and the correlations between parameters. Crosshairs in the figure mark the best-fit values, which are defined as the means of the marginalized distributions. We summarize these best-fit parameters in Table 1. The model parameters are highly non-degenerate and the MCMC chain converges to the best-fit model.

Figure 6 illustrates the comparison between the observed data and the best-fit model predictions (with parameters listed in Table 1, derived from Eq. 4). We bin the data in Galactocentric radius (9, 11, and 12 kpc) and age intervals (0.1 to 0.5 Gyr). The shaded areas represent the standard deviation of the observed data. This was calculated using a bootstrap analysis with 500 samples, where each sample was generated by resampling the original data with replacement. For each bootstrap sample, we calculated the mean and median $J_{z}$ values within defined radial and age bins. We use the MCMC chains to get the uncertainties on the model, depicted as vertical lines around the dashed lines.

The vertical action $p\left(J_{z}\mid{R}_{\mathrm{GC}}\right)$ increases linearly with age. Even though we assumed stars are born cold, they show non-zero vertical action at birth. The vertical heating rate increases at large $R_{GC}$ (Figure 6).

The left panel of figure 7 shows the values of $J_{z}$ calculated from the data. The central panel shows the median $\widehat{J}_{z}$ calculated by the model with the best-fit parameters shown in Table 1. Both panels show a gradient in the diagonal direction as $J_{z}$ increases with respect to the ages and Galactocentric radii. The right panel shows the difference between the model and the data, normalized with respect to the value of the model $\widehat{J}_{z}$.

Using a substantial data set of young massive stars drawn from Zari et al. (2021) and cross-matched with LAMOST to obtain radial velocities and ages, we measured the vertical motions of young stars ($<500$Myr) in the Galactic disk over an unprecedentedly wide range in Galactocentric radii (8 - 13 kpc). Previous studies over such a wide range in $R_{\mathrm{GC}}$ were based on red clump stars (Ting & Rix, 2019), which represent older populations (1.5-10 Gyr). Our comparisons with Ting & Rix (2019) (see Figure 6) are not directly equivalent due to methodological differences. Specifically, our model does not account for selection effects, and neither our model nor the Ting & Rix (2019) observations account for dust extinction. This is the primary reason we initially did not consider the potential flaring in our analysis; instead, our focus was on the evolutionary trends in $J_{z}$ rather than absolute values. Note that we fit only one model with functional radial dependence rather than radial bins, so any rigidity of the model might result in a small offset between the mean of the data and the best-fit model, as can be seen for $R=12$ kpc. Overall, the agreement between our model predictions and our observed data is satisfactory; see Figure 7. Notice, the model of Ting & Rix (2019) is also shown as star points and differs by a factor $\sim 1.4$ ; this difference could result from using different tracers or a different survey selection.

Our analysis centered on characterizing vertical motions through vertical action ($J_{z}$), fitting the data with a model that predicts the mean vertical action $\widehat{J}_{z}(\tau~{}|~{}R_{GC})$. We found that (1) the stars’ mean vertical action $\widehat{J}_{z}(\tau~{}|~{}R_{GC})$ increases nearly linearly with age at all $R_{GC}$, and (2) that $\widehat{J}_{z}(\tau~{}\mid~{}R_{GC})$ increases towards larger Galactocentric radii for any given age, including at birth. Here, we interpret only the former: the age-dependence of the mean vertical action, but not the latter, which can arguably be tied to extinction decreasing our probability to observe young stars in the midplane or to the Galactic warp.

To effectively model the observed relationship between mean vertical action with respect to age and Galactocentric radius, we explored several models, starting with a simplified approach that included only three terms: a constant term, a linear term with respect to radius, and a quadratic term with respect to radius. Subsequently, we extended our investigation to more intricate models incorporating quadratic dependencies on age. Among the range of models explored, the model presented in this paper emerged as the optimal choice, giving the most robust and consistent fit to the data.

As stated in earlier sections (see Section 2), we have chosen a fixed uncertainty value of $10\,\text{km/s}$ for the line-of-sight velocities ($v_{\mathrm{los}}$) of all the stars in our analysis. Nevertheless, as part of our investigation, we conducted a test by incorporating the uncertainties provided by LAMOST into the MCMC analysis and found the same best fit parameters within 1 sigma.

It is tempting and traditional to interpret the observed increase in vertical action with age in terms of heating because the near-linear age dependence of the mean vertical action ($\gamma=0.73\approx 1$) is consistent with heating by giant molecular clouds and with observations of older stellar populations (Ting & Rix, 2019). Such vertical heating could be driven by a combination of giant molecular clouds, external perturbations such as satellites and flybys, or by internal asymmetries. However, internal heating mechanisms do not appear to be likely explanations of the data, or at least not uniquely, because the data require that such heating would need a heating rate constant with $R_{GC}$, even though the density of potentially scattering molecular clouds should drop off with $R_{GC}$ (e.g., Rice et al., 2016).

There may be alternative explanations for our observed $\widehat{J_{z}}(\tau~{}|~{}R_{GC})$. We discuss them below; however, exploring these would necessitate a careful, in-depth analysis and model selection that extends beyond the scope of this work. For instance, stars at increasingly larger radii may originate from gas that is not co-planar with the inner Galactic disk, such as in a warped gas disk (refer to Figure 8 in Soler et al., 2022). Alternatively, the mean increase in vertical action with age could be attributed to the azimuthal mixing of stars over time. Azimuthal mixing refers to the redistribution of stars along their orbital paths around the Galactic center due to differential rotation. If a warp is present in the disk, the mid-plane is not centered on $z=0$ everywhere, but it is a function of radius and azimuth. Stars born in this warped region will have initial vertical positions ($z$) that are not zero, and therefore they will have computed actions $J_{z}$ that are not zero. Their angular velocities decrease with increasing radius ($\Omega_{\phi}\propto 1/R$). Consequently, stars in the outer disk take longer to reach the solar neighborhood. So if there is a warp, then the vertical action at birth ($J_{z}$) is also a function of azimuth ($\phi$). Even in the absence of heating, if we select stars in the solar neighborhood, there will be a relation between their age and their birth azimuth ($\phi=\phi_{0}+\Omega_{\phi}t$), which translates into a relation between $J_{z}$ and age ($t$). Pursuing further explanation would require modeling or simulation analysis, which lies beyond the scope of this paper. Should there be a trend of $J_{z}$ with azimuth at birth, subsequent azimuthal mixing could manifest as an apparent age trend within a localized sample. Pursuing this explanation would require modeling or simulation analysis, which lies beyond the scope of this paper.

Massive stars are often found in binary systems (Sana et al., 2014). Consequently, some stars in our sample could be part of such systems, potentially influencing the kinematics. However, binarity does not seem to be a plausible source of the observed trend of $J_{z}$. There is no evidence to suggest that the impact of binarity escalates with advancing age or at greater Galactocentric radii. Similarly, after massive binaries form, one of the companions may explode in a supernova before the other, thereby unbinding from its companion and leaving it walking (or running) at the speed it had when they were bound. Typical ejection velocities range from 5-20 km$,\mathrm{s}^{-1}$ and exhibit only a weak dependence on the mass (and thus, the age) of the ejected star (Renzo et al., 2019). Thus, this probably not be the origin of our observed $J_{z}-\tau$ trend.

Our results are likely caused by a combination of factors, including vertical heating by giant molecular clouds and the structural evolution of the Galactic disk. There is a complex interplay of multiple mechanisms and this complexity warrants a cautious interpretation. In this study, we have not considered the potential under-representation of mid-plane stars in our sample due to interstellar extinction or effects of LAMOST selection function. Figure 1 shows some lack of data around $z\approx 0$ and the decrease in the number of stars beyond 11 kpc. Qualitatively, interstellar extinction, leading to an increasing scarcity of mid-plane stars with distance from us, would result in an increase of the mean vertical action with distance from the Sun (and thus, in this case, towards the outer disk). Consequently, without modeling of dust and selection functions, the radial dependence of mean vertical motion should be interpreted with caution. However, using a toy model, we confirmed that our main finding—that the mean vertical action increases nearly linearly with age across all Galactocentric radii—remains robust. We produced a toy disk model, generating stellar 3D positions and velocities, and where the mean vertical action increases with age. Then, we applied a toy dust model where stars with galactic latitudes larger than some threshold do not make it into the mock dataset. From this experiment, we found (1) a non-physical increase in $\langle J_{z}\rangle(R)$ with radius, as expected and (2) no effect on the $\langle J_{z}\rangle(\tau)$ relation, which gives us confidence that our main result is robust.

Finally, we explored only a small volume near the Sun due to the reduced sky coverage of the sample, as determined by the footprint of the LAMOST survey. The all-sky coverage of future spectroscopic surveys, such as SDSS-V (Almeida et al., 2023; Kollmeier et al., 2017), will allow to extend this to a larger volume of the Milky Way disc. Enhanced coverage of the disk in Galactic azimuth will facilitate distinguishing whether the observed trends arise from a combination of heating processes or intrinsic birth properties, such as disk asymmetries or warping.

This work was partially carried out during the MPIA Summer Internship, funded by the Max Planck Society.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

This work has used data from the Guoshoujing Telescope, where the LAMOST survey was executed. LAMOST is a National Major Scientific Project built by the Chinese Academy of Sciences(Zhao et al., 2012; Deng et al., 2012).

NF acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), [funding reference number 568580] through a CITA postdoctoral fellowship and acknowledges partial support from an Arts & Sciences Postdoctoral Fellowship at the University of Toronto.

This research made use of Astropy (Robitaille et al., 2013; Collaboration, 2022, 2018, 2013), matplotlib (Hunter, 2007), numpy (Harris et al., 2020), and SciPy (Virtanen et al., 2020). We extend our gratitude to the reviewer for their insights and constructive feedback.