Fundamental relations for the velocity dispersion of stars in the Milky Way

We explore the fundamental relations governing the radial and vertical velocity dispersions of stars in the Milky Way, from combined studies of complementary surveys including GALAH, LAMOST, APOGEE, the NASA $Kepler$ and K2 missions, and $Gaia$ DR2. We find that different stellar samples, even though they target different tracer populations and employ a variety of age estimation techniques, follow the same set of fundamental relations. We provide the clearest evidence to date that, in addition to the well-known dependence on stellar age, the velocity dispersions of stars depend on orbital angular momentum $L_z$, metallicity and height above the plane $|z|$, and are well described by a multiplicatively separable functional form. The dispersions have a power-law dependence on age with exponents of 0.441$\pm 0.007$ and 0.251$\pm 0.006$ for $\sigma_z$ and $\sigma_R$ respectively, and the power law is valid even for the oldest stars. For the solar neighborhood stars, the apparent break in the power law for older stars, as seen in previous studies, is due to the anti-correlation of $L_z$ with age. The dispersions decrease with increasing $L_z$ until we reach the Sun's orbital angular momentum, after which $\sigma_z$ increases (implying flaring in the outer disc) while $\sigma_R$ flattens. The dispersions increase with decreasing metallicity, suggesting that the dispersions increase with birth radius. The dispersions also increase linearly with $|z|$. The same set of relations that work in the solar neighborhood also work for stars between $3


INTRODUCTION
Most stars of the Milky Way's disc population are thought to be born on roughly circular orbits with very low velocity dispersion. However, a significant fraction of these stars are observed to have very high velocity dispersion, suggesting that they must have undergone significant dynamical evolution. Studying the velocity distributions for disc stars can therefore not only shed light on the dynamical history of the Galaxy, but also on the dynamical processes that have shaped the present day distribution of stars. Moreover, velocity dispersion relations are an essential ingredient for constructing analytical models of the Galaxy. Multiple studies have sought to characterize the velocity dispersion of stars in the Milky Way disc and to explain them using dynamical models. Despite much progress, however, many open questions remain.

Background
It has been known for a long time that the velocity dispersion of disc stars increases with age in the solar neighbourhood. One of the earliest attempts to explain this observation dates back to Spitzer & Schwarzschild (1951), who examined in-plane motions of disc stars. They showed that the total dispersion of all components, σ tot , increases with age τ due to scattering from massive clouds with a power-law dependence, τ βtot , with exponent β tot = 0.33. Remarkably, they inferred the presence of giant molecular clouds (GMCs) long before they were observed directly. Later, Lacey (1984) generalized this result by including vertical motions. However, the predictions of his model conflicted with observations. First, Lacey concluded that β tot = 0.25, whereas, observations suggest that solar-neighborhood stars have β tot between 0.3 and 0.5 -or more precisely that β z ranges from 0.35 to 0.6 and β R ranges from from 0.19 to 0.35 (Nordström et al. 2004;Aumer & Binney 2009;Sharma et al. 2014;Mackereth et al. 2019). Secondly, the overall heating rate derived (which determines the dispersion of older stars) was too low. Finally, Lacey predicted the ratio of vertical to radial dispersion, σ z /σ R , to be 0.8, which is higher than observed ratio for stars in the solar neighborhood (0.5-0.6). Hänninen & Flynn (2002) used N-body simulations to confirm some of the findings of Lacey (1984). They found that for scattering from GMCs, β tot is indeed less than 0.33, or more precisely that β R = 0.20, β z = 0.26 and β tot = 0.21, which is even less than the value predicted by Lacey. The overall heating rate was also confirmed to be low. Specifically, with the surface density of GMCs set to the presentday value in the solar neighborhood (5 M /pc 2 ), the predicted velocity dispersion for the oldest stars was less than that of the observations by about 15 km/s. Additionally, Hänninen & Flynn (2002) found that the ratio σ z /σ R depends on the number density of GMCs. For the GMC number density of 5 M /pc 2 , σ z /σ R was 0.5 in rough agreement with the observed value, but much less than the value of Lacey (1984). The value of σ z /σ R computed by Lacey (1984) was large because an isotropic distribution of star-cloud impact parameters was assumed. When the anisotropy in the impact parameters is taken into account, σ z /σ R is 0.62 in the steady state (Ida et al. 1993;Shiidsuka & Ida 1999;Sellwood 2008). This value of the ratio agrees well with the observations. For example, for 10 Gyr old stars, Aumer & Binney (2009) report the ratio to be 0.56 using the Geneva Copenhagen survey (GCS - Nordström et al. 2004), while Sharma et al. (2014) report values of 0.59 and 0.65 using the GCS and RAVE surveys respectively.
The inability of GMC scattering models to match the observed data prompted the exploration of other mechanisms to excite random motions. Transient spiral structures are one such mechanism. They lead to potential fluctuations in the disc that can heat up disc stars (Barbanis & Woltjer 1967;Sellwood & Carlberg 1984;Carlberg & Sellwood 1985). De Simone et al. (2004) showed using numerical experiments in two dimensions that spiral arms alone can lead to β R in the range 0.25 to 0.5, and a heating rate such that the value of σ R is consistent with observations. However, spiral scattering is too inefficient to increase the vertical dispersion (Sellwood 2013;Martinez-Medina et al. 2015). This led Jenkins & Binney (1990) to argue that a combination of spiral structure and GMC heating could explain the velocity dispersion of observed stars, though their predicted β z (0.3) was too low and the predicted β R (0.5) was too large.
One way to resolve the discrepancy between the predicted and the observed values of β z is to accommodate a scattering environment that is evolving with time. This consequently means that the velocity dispersion as a function of age (agevelocity relation or AVR) for stars in the solar neighborhood is not same as the evolution of velocity dispersion with time of stars born together (heating history) at a given time in the past. In other words, the AVR is the compilation of the end of the heating history of stellar populations born at different times. Specifically, due to the much higher gas fraction in the early Galactic disc, the contribution of GMC scattering is expected to decrease with time, which has been shown to lead to β z being close to 0.25 for the heating history but greater than 0.4 for the AVR (Aumer et al. 2016a;. There are other physical processes that can heat up the disc, and so shape the AVR. The Milky Way hosts a bar that can heat disc stars, as demonstrated in isolated disc/bar/bulge simulations (e.g., Saha et al. 2010), particularly near the strong Lindblad resonances. The same effect is seen in cosmological N-body simulations, where bars emerge within the evolving disc (Grand et al. 2016). But here there is an added contribution from the disc being bombarded by orbiting satellites (Velazquez & White 1999). The AVR can also be shaped by the fact that the intrinsic velocity dispersion was higher at earlier times, as reported by Hα emission of gas in external galaxies at high redshift (Förster Schreiber et al. 2009;Wisnioski et al. 2015). However, Hα emission tracks ionized gas and it is not yet clear, if stars, which form out of cold gas, also have high velocity dispersion like the ionized gas. There is now strong evidence, both theoretical (Sellwood & Binney 2002;Roškar et al. 2008) and observational, that stars migrate from their place of birth. The observational evidence comes from the presence of low-eccentricity and super-metallicity stars in the solar neighborhood-a realization that dates back to at least Grenon (1972) (see also Kordopatis et al. 2015;Hayden et al. 2020). Because the heating rate is higher in the inner regions than in the outer regions, one has to take migration into account when modelling the AVR at a given location. However, in a comprehensive review of stellar migration, Minchev (2016) argues that migration, on average, generally does not lead to disc heating.

Modelling the Physical Mechanisms of Disc Evolution
Given that various physical processes can play a role in disc evolution, it is imperative to study them both individually (to assess their relative importance) and in combination (to see the full effect of them acting together). The seminal work by Aumer et al. (2016a) highlights the utility of this approach. They analyzed N-body simulations that had spiral arms, GMCs, a bar, and growing discs -and were successful in reproducing the AVR in the solar neighborhood of the Milky Way. However, this model was only compared with observations in the solar neighborhood, and possible dependencies on metallicity and angular momentum were not considered. We show for the first time that these properties taken together can place even stronger constraints on the models.
There are also other good reasons for studying the dependence of velocity dispersion on metallicity and angular momentum. For a given age, the metallicity provides a way to tag the birth radius of a star, which provides leverage on the process of radial migration (Bland-Hawthorn et al. 2010;Frankel et al. 2019). The angular momentum L z provides the mean radius where a star spends most of its history, and so (unlike the present radius R) is a more useful indicator of the amount of scattering the star has undergone. Moreover, there are strong theoretical reasons to prefer L z over R. For an axisymmetric system L z is a constant of motion and, Jeans' theorem tells us that the phase space density of a system in dynamical equilibrium should only depend on constants of motion (Binney & Tremaine 2008). Furthermore, specifying the dispersion as a function of L z paves the way for constructing better analytical models of the Galaxy, e.g., dynamical models based on actions by Binney (2012) and Sanders & Binney (2015).
Given that the velocity dispersion of a population of stars can depend on a number of stellar properties, it is important to come up with a useful way to characterize the velocity dispersion from observations such that it can test theoretical models. The selection function of a survey will tend to leave its imprint on the measured velocity dispersions (Sharma et al. 2014), which makes it difficult to compare and combine results from different surveys, and also to compare observed dispersions with model predictions. If the velocity dispersion σ only depends on a set of observables X, then knowing X is sufficient to characterize the dispersion irrespective of the selection function. This brings us to the question of identifying the fundamental relations governing the dispersion, i.e., what is a suitable choice for the set of observables/variables X, and how does the velocity dispersion depend on them. Intuitively (as discussed above) the dispersion should be governed by age, metallicity and angular momentum. However, the joint dependence of dispersion on these properties has not been studied before, and this is what we address in this paper. Since stars migrate from their place of birth, strictly speaking, the amount of scattering a star experiences will also depend upon the evolutionary history of its angular momentum L z , but we do not consider this as we do not have any observable that tracks this information. Additionally, each observational survey and age estimation technique has its own systematics, and no attempt has been made to characterize such systematics, and we also address this issue.

Disc Evolution in the Age of Massive Galactic Surveys
A number of large observational surveys cataloging the detailed properties of a huge numbers of stars in the Milky Way mean that we are better poised now than ever before to unravel the fundamental velocity dispersion relations. These data sets probe stellar kinematics well beyond the solar neighborhood, and hence provide more coverage of the angular momentum and metallicity dimensions. Additionally, these data sets have large sample of stars that allows the additional dependence on metallicity and angular momentum to be studied robustly. The combination of Gaia DR2 astrometry and accurate radial velocities from ground-based spectroscopic surveys, provides precise six-dimensional phase space information for a large number of stars. Spectroscopic surveys such as GALAH and LAMOST, mean that it is now possible to get reliable age estimates for a large number of main sequence turn-off (MSTO) and subgiant stars in the solar neighborhood, a significant improvement when compared with photometry-based ages (e.g., Bland-Hawthorn et al. 2019). Asteroseismology from missions like Kepler and K2 has opened the door to estimating the ages of intrinsically bright giant stars allowing us to study the velocity dispersions well beyond the solar neighborhood. Ground-based spectroscopic surveys also provide elemental abundances, with which we can tag stellar populations that were born at the same time and same place (Freeman & Bland-Hawthorn 2002).
Some attempts have already been made to characterize the velocity dispersions using the large observational surveys that probe the velocity dispersion beyond the solar neighborhood and below we summarize four such studies. Sanders & Das (2018) used data from multiple observational surveys and multiple stellar types (about 1.2 million stars), but did not consider potential systematics between the different surveys used. They studied the dependence of velocity dispersion on age and radius R, but ignored the dependence on angular momentum, metallicity, and z. They found that velocity dispersion decreases exponentially with R out to the solar Galactic radius, and that beyond this σ z tends to increase, while σ R tends to flatten out. The velocity dispersion was found to grow as a power law with age, with exponent β R ∼ 0.3 and β z ∼ 0.4.  and Mackereth et al. (2019) both studied the relationship between age and kinematics using APOGEE-DR14 data and ages estimated using a neural network model, which was trained on asteroseismic ages from the NASA Kepler mission. Ting    studied the dependence of the expectation value of the vertical actionĴ z as a function of age and average radius R GC , which was defined to be the mean of the birth radius and the current radius. Crucially, they did not consider the dependence on birth radius and angular momentum separately, and also they did not study the in-plane kinematics. They found that the expectation value of birth actionĴ z,0 is constant until about R GC = 10 kpc but rises beyond that. Assuming the following approximate relations, which are valid under epicycle approximation,Ĵ z ∝ σ 2 z /ν, and ν ∝ √ Σ, where ν is vertical oscillation frequency and Σ is mass surface density, we interpret their results as follows. The vertical dispersion falls off exponentially with R GC until 10 kpc, but beyond that it flattens or increases. They also found that the vertical dispersion increases with age as a power law with exponent β z ranging from 0.5 to 0.65. Mackereth et al. (2019) studied the velocity dispersion as a function of Galactocentric cylindrical coordinates R and z for stars binned by age, [Fe/H] and [α/Fe]. A quadratic model was assumed for dependence on z, and an exponential model for dependence on R. We note that the exponential model might be inappropriate, given that Sanders & Das (2018) find the dependence of dispersion on R to be exponential only for R < 8 kpc, but flat and even rising for R > 8 kpc. Mackereth et al. (2019) found that for young stars the dispersions and the ratio σ z /σ R increase with height |z|, which they attribute to stronger heating by spiral structure in the plane and the relatively longer time scale for GMCs to redirect random in-plane motion to vertical motion. For a given age they find that the dispersions are higher for those monometallicity populations that have a larger mean orbital radius. But given that they study their stars by binning them up in mono-metallicity populations and the fact that metallicity is anti-correlated with mean orbital radius, we note that the observed trend with mean orbital radius is indistinguishable from a trend with metallicity. Minchev et al. (2018) studied the vertical dispersion σ z , of about 500 solar-neighborhood stars as a function of birth radius and age. The dependence of σ z on angular momentum was not studied. For a given age, the σ z was found to vary with birth radius such that it has a slope, which is positive for old stars (age greater than 8 Gyr), flat for intermediate age stars, and slightly negative for young stars (age less than 4 Gyr).

DATA
In this paper, we mainly make use of data from the LAM-OST Zhao et al. 2012) and GALAH spectroscopic survey (De Silva et al. 2015). We also use the APOGEE-DR14 spectroscopic survey (Majewski et al. 2017), but only for the purpose of studying systematic effects. We used the LAMOST-DR4 value added catalog from Xiang et al. (2017b), for radial velocity, T eff , log g, [Fe/H], [α/Fe], and distance. For LAMOST stars, we used two types of stars, the MSTO stars and the red-giant (RG) stars. The ages for the LAMOST-MSTO sample were taken from Xiang et al. (2017a) and for the LAMOST-RG-CN sample were taken from Wu et al. (2019). The LAMOST-RG-CN sample consists only of red giant branch stars (red clump stars are not included), with ages derived from spectroscopic C and N features. For the GALAH survey, we also used two types of stars, the MSTO stars and the RG stars. More precisely, we make use of the extended GALAH catalog (GALAH+), which also includes data from TESS-HERMES (Sharma et al. 2018) and K2-HERMES  surveys that use the same spectrograph and observational setup as the GALAH survey. The RG stars that we use have asteroseismic information from the NASA K2 mission and their spectroscopic followup was done by the K2-HERMES survey, hereafter they are referred to as GALAH-RG-K2. We note that the spectroscopic analysis of GALAH-DR2 ) was based on a machine learning model trained on a set of 10,605 stars analysed in detail using the SME code (Valenti & Piskunov 1996;Piskunov & Valenti 2017). In this paper, we exploit parameters from GALAH-iDR3, an internal data release where every star has been analysed using SME and incorporates Gaia-DR2 distance information (Gaia Collaboration et al. 2018;Lindegren et al. 2018). A full discussion will be presented in a forthcoming paper and the results will be available as part of GALAH-DR3.
To select stars with reliable ages, we adopt the following selection function for MSTO stars, (1) For the RG stars, we adopt the following selection criteria, The ages and distances for the GALAH-MSTO and GALAH-RG-K2 stars are computed with the BSTEP code (Sharma et al. 2018). BSTEP provides a Bayesian estimate of intrinsic stellar parameters from observed parameters by making use of stellar isochrones. For results presented in this paper, we use the PARSEC-COLIBRI stellar isochrones (Marigo et al. 2017). For the GALAH-MSTO stars, we use the following observables, T eff , log g, [Fe/H], [α/Fe], J, Ks and parallax. For the GALAH-RG-K2 stars, in addition to the above observables, we use the asteroseismic observables ∆ν and ν max . These stars were observed by the NASA K2 mission as part of the K2GAP program (Stello et al. 2015) and includes stars from campaigns 1 to 15. The asteroseismic analysis is conducted with the method by Kallinger et al. (2010Kallinger et al. ( , 2014, known as the CAN pipeline. ∆ν and ν max for the model stars in the isochrones are determined with the ASFGRID code Sharma et al. (2016) that incorporates corrections to the ∆ν scaling relation suggested by stellar models. A summary of different data sets used in this paper is given in Table 1.
Transformation from heliocentric to Galactocentric coordinates is done assuming R = 8.0 kpc (Reid 1993), z = 0.025 kpc, Ω = 30.24 kms/s/kpc (Reid & Brunthaler 2004), U = 10.96 km/s and W = 7.53 km/s (Sharma et al. 2014). We use a right handed UVW coordinate system, where U points towards the Galactic center, V is in the direction of the Galactic roation and W points towards the North Galactic pole. The transformation is carried out with the following heliocentric quantities: Gaia-DR2 angular position and proper motions, spectroscopic radial velocities, and spectro-photometric distances. Where needed we assume the circular velocity at Sun, Θ , to be 232 km/s (Sharma et al. 2014).

METHOD
The dispersion σ v of velocity v (for either v R or v z ), is assumed to depend on the stellar age τ , angular momentum L z , metallicity [Fe/H], and vertical height from the disc midplane z, via the following multiplicatively separable functional form Here, X = {τ , L z , [Fe/H], z} is a set of observables that are independent variables and and The functional forms were chosen based on a preliminary analysis of the trends with respect to each observable. The σ v has a power law dependence on age, with β v denoting the exponent. The age relation has a finite birth dispersion for stars younger than 0.1 Gyr. The σ v falls off exponentially with L z with scale λ L,v , but at large L z it is allowed to rise as L 2 z (to account for flaring) and this rise is controlled by α L,v . The Dispersion  Dispersion |z| with gradients γ [Fe/H],v and γ z,v respectively. The σ 0,v is a constant that denotes the velocity dispersion for stars lying in the midplane with solar metallicity, solar angular momentum (L z, = Ω R 2 ) and an age of 10 Gyr. The likelihood of the observed velocities v = {v 0 , ...v N } for a sample of N stars can be written as (8) with vi being the uncertainty corresponding to the observed velocity v i of the i−th star. Here, N (v|µ, 2πσ 2 denotes the distribution of a random variable v sampled from a normal distribution with mean µ and variance σ 2 . We find the maximum likelihood estimate (MLE) of θ v by using the Nelder-Mead algorithm as implemented in the python package scipy.optimize.minimize. The MLE values of θ v for the velocity components v z and v R are given in Table 2. Also given alongside are uncertainties, which were estimated using bootstrapping. The data used for estimating θ v contained an equal number of stars from the LAMOST and GALAH surveys.
To see the velocity dispersion profiles f x (x) that result from the observed data corresponding to each independent variable x, we must factor out the dependence on other indepen- Dispersion dent variables in the set X. We accomplish this by binning the stars in x and computing the dispersion of v/σ v, [x] in each bin to get a profile of the dispersion as a function of x. Here, is the complementary velocity dispersion that includes all independent variables in set X except x. Dispersion the bottom panels explore the radial dispersion. The dispersion is assumed to be a function of multiple independent variables, and from left to right each panel shows the effect of one independent variable at a time. The dispersions decrease exponentially with L z up to about the solar angular momentum, beyond that σ z starts to increase ( Figure 1a) whereas σ R flattens ( Figure 1e). The dispersions increase with age and are well described by a power law, with the exponent being higher for σ z as compared to σ R (Figure 1b Figure 2 demonstrates that all of the above discussed trends are independent of the location R and are valid for 3 < R/kpc < 20. Data from different surveys and stellar types are shown separately in Figure 1. Figure 1a and Figure 1b show some systematic differences, but overall the different data sets are all found to be consistent with the same relationship. Agreement between the GALAH and LAMOST results suggests that their spectroscopic parameters do not not have any strong systematics with respect to each other. Agreement between data sets of different stellar types (MSTO and RG), suggests there are no strong systematics related to stellar types. This is very reassuring and useful given the fact that for different stellar types very different age estimation techniques are used.

Variation of basic trends with respect to other independent variables
In Figure 3 to Figure 11, we show the residual dependence of each relation on other independent variables. In each figure, we plot the observed relation corresponding to one independent variable, such as the relation f τ , by binning stars in another independent variable, e.g., L z or [Fe/H]. In general we find that the relations vary very little with respect to other independent variables, suggesting that modelling the dispersion as a product of multiple independent functions, as given by Equation 3, is a good approximation. However, slight variations can be seen. Interestingly, the variations that we see are in most cases systematic and we discuss these systematic trends below.
The f Lz relation varies with age such that the younger stars have higher relative dispersion (Figure 3). Here by relative dispersion we mean dispersion relative to the derived relations. The f Lz relation varies with [Fe/H], such that for high L z , the metal rich stars have systematically higher relative dispersion (Figure 4). The f τ relation does not vary for stars with age less than 8 Gyr, but for older stars the relative dispersion is systematically lower for the high-L z stars ( Figure 5) and for high-metallicity stars ( Figure 6). This could be because both old high-L z stars and old high-metallicity stars are rare and hence an old star bin is more likely to be contaminated by young stars (due to age uncertainties), which will lower the overall dispersion in that bin given that young stars have low dispersion.
The f [Fe/H] relation flattens with increasing L z (Figure 7). The relation does not vary much with age, with the exception that for old stars the relative dispersion is lower at the high metallicity end. (Figure 8). This again could be due to old and high-metallicity stars being rare and hence an old star bin is more likely to be contaminated by young stars. The f z relation for vertical velocity dispersion does not seem to vary Dispersion much with L z (Figure 9a), however, the f z relation for radial velocity dispersion flattens with increase in L z (Figure 9b). The f z relation does not vary much with age ( Figure 10) or [Fe/H] (Figure 11). However, for stars younger than 4 Gyr the relationship is much steeper (Figure 10).

Systematics between surveys
In Figure 1 we presented results from four different data sets. Two were based on the LAMOST spectroscopic survey, with one made up of MSTO stars and other made up of RGB stars. The other two data sets were based on the GALAH spectroscopic survey, with one made up of MSTO stars and other of RG stars having asteroseismic information from K2. Another large spectroscopic survey that we did not use in Figure 1 was APOGEE. In Figure 12a, we plot APOGEE results for the asteroseismic sample from Kepler (Pinsonneault et al. 2018). It shows that the observed AVR is shifted with respect to our empirical relations. Examination of stars in between the APOGEE and GALAH/LAMOST data sets, revealed that the APOGEE iron abundances were systematically higher by 0.1 dex. However, this is still not enough to account for the difference seen in Figure 12a. A further increase in the age of APOGEE stars by 10% is required to bring the the APOGEE-RG-KEPLER data set into agreement with the other datasets. These systematic offsets are the reason APOGEE data was not used in the analysis presented in Figure 1. Figure 12b shows that with these changes the APOGEE-RG-KEPLER data set can also be brought into agreement with the GALAH-RG-K2 data.
In Figure 1, the LAMOST-MSTO and LAMOST-RG-CN data sets show significant flattening of the AVR for age greater than 10 Gyr, while such a flattening is not seen for the GALAH-MSTO stars. The flattening for LAMOST data sets could be due to larger uncertainties on age estimates in them. Figure 12c shows results for the LAMOST-MSTO and LAMOST-RG-CN data sets, using age estimates from Sanders & Das (2018). No flattening is seen here. This could be because the Sanders & Das (2018) ages are more precise than LAMOST ages for the older stars. However, it should be noted that Sanders & Das (2018)  Dispersion Figure 9. Velocity dispersion as a function of distance |z| from the mid-plane of the Galaxy, for stars lying in different angular momentum bins. The shaded region denotes the 16 and 84 percentile confidence interval. The dependence on other independent variables (Lz, age, and [Fe/H]) have been factored out. The description of the panels and the data set used are same as in Figure 3. The relationship for radial dispersion shows a mild variation with the Lz, with the slope becoming flatter with the increase of Lz. on height above the plane, which can increase the age precision for older stars, but is not ideal to study trends with height above the plane, as we do in this paper.

The role of angular momentum in shaping the solar-neighborhood AVR
Angular momentum plays a critical role in shaping the AVR of stars observed in the solar neighborhood. It has been claimed in some previous studies that the AVR deviates from a power law with an abrupt increase for old stars (Freeman 1991;Edvardsson et al. 1993;Quillen & Garnett 2001). Since, most old stars in the solar neighborhood belong to the thick disc, this tentatively suggests that the kinematics of the thick disc stars is different from that of the thin disc stars. We show in Figure 13 that this apparent break is due to a systematic variation of angular momentum with age, since older stars have low angular momentum and low angular momentum stars have higher velocity dispersion. Once Dispersion Figure 10. Same as Figure 9 but for stars lying in different age bins. The slope is higher for younger stars. again, when a variation of angular momentum is allowed for, all stars seem to be consistent with a universal AVR. However, it is still not clear as to why the angular momentum decreases with age. It could be due to inside out formation of the disc and radial migration of stars from the inner disc and needs to be investigated in future.

Dependence of dispersions on age
We find β z = 0.44 and β R = 0.26, which is in good agreement with predictions of simulations by Aumer et al. (2016a), where the effects of spiral arms, GMCs and a bar is taken into account. For the GALAH-MSTO stars, even the old thick disc stars satisfy this relationship. The LAMOST-MSTO, LAMOST-RG-CN and GALAH-RG-K2 all show saturation for age greater than 10 Gyr. This could be due to significant uncertainties in ages for the older stars in samples other than GALAH-MSTO. For example, the uncertainty of asteroseismic ages is known to increase with age and is predicted to be around 30% for RGB stars and even higher for red-clump stars (Silva Aguirre et al. 2018). Since old stars are typically rare, an old star bin is more likely to be contaminated by young stars (due to age uncertainties), which will lower the overall dispersion in that bin given that young stars have low dispersion.
As discussed in Aumer et al. (2016a), the exponents β z and β R of the AVR depend in a complex way upon the whole dynamical history of the Galaxy. This is because there are at least two major scattering agents (spiral structures and GMCs) and the strength of scattering due to them changes with time. Spiral structure is mainly responsible for in-plane scattering, while GMCs contribute to both in-plane and vertical scattering. Spiral structure drives up σ R fairly rapidly, increasing the Toomre stability paramater Q and making the disc stable. This makes the ratio σ z /σ R very small initially. Thereafter, on a longer time scale, scattering from GMCs increases both σ R and σ z . For scattering from stationary fluctuations the exponent β is predicted to be around 0.25, with β z being slightly higher than β R (Hänninen & Flynn 2002). This is seen in the heating history of coeval populations. However, β z for the AVR is much higher because the overall efficiency of heating due to GMCs is higher at earlier times ( Figure-7 of Aumer et al. 2016a). At earlier times, the star formation rate is high and the stellar disc mass is low and this makes the GMC mass fraction higher, which in turn increases the efficiency of GMC scattering. Dispersion

Dependence of dispersions on angular momentum
Without any loss of generality, in what follows, we discuss our results in terms of guiding radius rather than angular momentum, with the guiding radius being defined as R g = L z /Θ (angular momentum divided by circular velocity at Sun). We find that for R g < R both dispersions fall off exponentially with R g . In general, the strength of the secular heating processes, like those due to spiral arms or GMCs, are expected to be proportional to the surface density of stars Σ, so the dispersion is expected to fall off with radius. Using some simple physically motivated arguments we now predict the radial scale length R σ of the exponential fall of dispersion with R. The vertical dispersion is expected to vary with surface density Σ, and scale height h z , as σ z ∝ √ Σh z (van der Kruit 1988). If h z is constant then the scale length of vertical dispersion, R σz , should be related to the scale length of stellar surface density, R d , as R σz = 2R d . Using this we estimate R σz = λ L,vz /Θ = 4.9 ± 0.2 kpc (adopting Θ = 232 km/s from Sharma et al. 2014), which is in good agreement with the theoretical prediction of 5.0 kpc (adopting R d = 2.5 kpc, see e.g., Robin et al. 2003;Jurić et al. 2008). For the radial dispersion we expect, σ R ∝ ΣR and hence R σR = R d . This follows from assuming the Toomre stability parameter to be constant throughout the disc and the rotation curve to be flat, which implies κ ∝ 1/R. However we find R σR to be 9.9 kpc, which is about four times larger than R d (assuming R d = 2.5 kpc). Hence, the observed scale length of σ R cannot be explained by a constant Q. For R g > R kpc, both dispersions break away from being purely exponential functions of R g , with σ z increasing and σ R flattening with R g . This is indicative of flaring in the outer disc for mono-age and mono-metallicity populations. Now, outside the solar radius significant flaring has been reported for all mono age populations, with the flaring being strongest for the youngest population (Mackereth et al. 2017). Flaring implies that σ z should fall off slower than √ Σ, be constant, or rise. Hence, our overall reported rise of σ z with R g in the outer disc, and also the fact that the rise is stronger for younger stars (Figure 3) is consistent with the findings of Mackereth et al. (2017) related to flaring.
The flattening of velocity dispersion with R g is easy to understand. Stars are thought to be born out of the inter-stellar medium with a birth dispersion of about 10 km/s, due to turbulence in the medium driven by the injection of energy from newly forming stars. Due to this non-zero lower bound on the dispersion of newly forming stars, at large R g the dispersion cannot keep on falling exponentially but will hit a floor and flatten.
We see in Figure 1b and Figure 1d that the flattening occurs at a smaller value of R g for σ z as compared to σ R . This is also easy to understand. Both σ z and σ R are exponential functions of R g , but the overall proportionality constant for σ R is larger than σ z . Additionally, the scale length for σ R is larger than that for σ z . Consequently, flattening due to a constant birth dispersion will set in at a lower value of R g for σ z than for σ R . The youngest stars also have the lower overall proportionality constant for dispersion, so they are expected to flatten earlier and this is visible in Figure 3. The fact that non-zero birth dispersion can lead to flattening of dispersions and consequently flaring has been nicely demonstrated by Aumer et al. (2016b) using N-body simulations of discs having spiral arms, GMCs and a bar.
What causes the dispersions to rise for σ z and why doesn't it also rise for σ R ? This needs to be investigated in future. Simulations by Aumer et al. (2016a), incorporating the effects of spiral perturbations, a bar, and GMCs, only predict a monotonic fall or flattening for σ z , but no rise of dispersion with radius (see their Figure 4). This suggests that some additional processes might be at play. For example, the interaction of the disc with orbiting satellites (Kazantzidis et al. 2008;Villalobos & Helmi 2008;Bournaud et al. 2009) is known to cause flaring in the outer disc. The infall of misalinged gas (Roškar et al. 2010;Sharma et al. 2012;) and reorientation of the disc axis  is also known to cause warps and consequently flaring. Interestingly, one of the disc galaxies (Au18) simulated in a cosmological context by Grand et al. (2016) shows a rise of σ z with R. Since, this simulation does not have GMCs, Grand et al. (2016) attribute the vertical heating to the bar and the effect of orbiting satellites.

The shape of the velocity ellipsoid
Our best fit relations (see Equation 3 and Table 2) can be used to estimate the ratio σ z /σ R . These relations suggest that there is a strong dependence of the ratio σ z /σ R on age, metallicity, angular momentum, and height above the plane and this is shown in Figure 14. More precisely, for 10 Gyr-old stars that are in the plane and have solar metallicity, the ratio first decreases with guiding radius to a minimum of 0.53 at R g = 7.25 kpc, and then increases with R g . The ratio is greater than 0.7 for R g > 12.0 kpc. Also, the ratio increases monotonically with decrease in metallicity and increase of height. If the vertical velocity dispersion is governed by scattering from GMCs, an equilibrium ratio of 0.62 is predicted, which can be attained by relatively old stars that had enough time to scatter. Simulations by Aumer et al. (2016a)  The ratio σz/σR as a function of guiding radius for stars of different age and metallicity. The ratio is estimated using the analytical model described by Equation 3 and parameters in Table 2. that spiral perturbations and GMCs can only lead to ratios σ z /σ R in range 0.5 to 0.7. But the fact that we find σ z /σ R to be greater than 0.7 in certain regions of the disc indicates that processes other than spiral structure and GMCs could be affecting the vertical dispersion of stars there.

Dependence of dispersions on metallicity
We see that velocity dispersions increase with decreasing metallicity for any given age and angular momentum (Figure 7 and Figure 8). Our current understanding of disc formation suggests that the ISM probably had a negative metallicity gradient for a significant fraction of its lifetime (Schönrich & Binney 2009;Minchev et al. 2018). This suggests that at any given age the metallicity should decrease with the birth radius of a star. Consequently, velocity dispersions should increase with birth radius. This result is counter intuitive, as naively we expect the dispersion to decrease with any type of radius (Section 5.2). Interestingly and importantly, we observe the dispersions to increase with birth radius for stars of all ages and angular momentum. One reason for this could be the conservation of vertical action, which will happen for stars of any age or angular momentum. Solway et al. (2012) demonstrated that vertical action is conserved for migrating stars. This conservation of vertical action leads to adiabatic heating/cooling of stars moving inwards/outwards (Minchev et al. 2012). To demonstrate this, let E z be the vertical energy, ν the vertical oscillation frequency, Σ the surface density and σ 2 z the vertical velocity dispersion. Using some standard assumptions and approximations it is easy to show that the vertical action J z = E z /ν = σ 2 z / √ 2πGΣ (for details see Section 1.3 and Minchev et al. 2012). When action is conserved, σ z ∝ Σ 1/4 . Hence, a stellar population born at radius R b with vertical velocity dispersion σ zb after migrating to or-bits with guiding radius R g will end up with a dispersion σ zb given by where R d is the scale length of surface density distribution. From this expression it is easy to see that outward migration leads to cooling and inward migration leads to heating. The question we are interested in is, for stars with a given angular momentum and age, is the dispersion of migrated stars higher or lower than that of non-migrated stars? The answer to this will depend on how σ zb varies with R b in Equation 11. For simplicity let us assume that σ zb varies with R b just as σ z varies with R g , i.e., Equation 3. Given that we find that the dispersion is flat or rising for R g > R , the inward migrators with R b > R are expected to be hotter than non-migrators. The outward migrators are also expected to be hotter because σ zb increases with decreasing R b with a scale length (about 2R d see Section 5.2) that is smaller than 4R d . But it was demonstrated by Vera-Ciro et al. (2014), using an idealized simulation of a galaxy with spiral arms, that migrating stars are preferentially of low vertical dispersion, given that they spend more time in the plane (see also Daniel & Wyse 2018). This bias was also shown to be present for discs in cosmological simulations (Grand et al. 2016). The bias makes it possible for the outward migrators to be cooler compared to non-migrators. Other processes can also lead to an increase of dispersions with birth radius. For example, effects like, disc satellite interaction, infall of misalignment gas or reorientation of the disc axis, and warping (that as discussed earlier lead to flaring in the disc), all predict the dispersions to increase with birth radius. Minchev et al. (2018) had also reported an increase of vertical dispersion with birth radius. However, they found the effect to be most prominent for older stars. The slope of variation of dispersion with radius was found to be positive for stars older than 8 Gyr and then flatten to zero at 6 Gyr and eventually turn negative for stars younger than 4 Gyr. In contrast, we find the slope to be positive for all ages, and additionally we also find a positive slope for the radial velocity dispersion, which they did not study. Mackereth et al. (2019) using mono age and mono metallicity populations reported an increase of vertical dispersion with mean orbital radius for low-[α/Fe] stars. For the radial dispersion they reported an increase with mean orbital radius only for stars younger than 4 Gyr. Given that mean orbital radius decreases monotonically with the metallicity for their populations, this means that their results can also be interpreted as an increase of dispersion with birth radius. In that case, for stars younger than 4 Gyr the slopes are opposite of that of Minchev et al. (2018).
Importantly, both Minchev et al. (2018) and Mackereth et al. (2019) had not factored out the dependence on angular momentum, which can have the opposite effect, because for stars with angular momentum less than solar angular momentum, the dispersion increases with decrease of angular momentum. This could be responsible for the differences between the above studies and differences with the results presented here.
If [α/Fe] abundance is assumed to be a good proxy for age, then our metallicity (or birth radius) trends can also be considered to be consistent with the findings of Hayden et al. (2020). They studied the velocity dispersion as a function of [α/Fe] abundance in different [Fe/H] bins using GALAH-DR2 data. They found that the vertical dispersion increases with decrease of metallicity for any given [α/Fe] just as we find the same effect for any given age.

Dependence of dispersions on height
We find that velocity dispersions increase with height for all angular momentum (Figure 9), ages ( Figure 10) and metallicities ( Figure 11). A positive slope is present for all ages but it is much steeper for stars younger than 4 Gyr (Figure 10). Also the slope for σ z is higher than that for σ R by about a factor of 2. This suggests that the ratio of σ z /σ R also increases with height. A non-zero slope implies that the populations defined by a specific age and metallicity are nonisothermal. Mackereth et al. (2019) also report a positive slope for low [α/Fe] populations, which was found to flatten with age. For high-[α/Fe] populations the slope was found to be zero, whereas we find the high-[α/Fe] stars to have a positive slope. As suggested by Mackereth et al. (2019), the non-isothermality could be related to the relatively large time scale for GMC heating as compared to the relatively fast in-plane heating by spiral arms. This is something that can be easily tested in idealized simulations by Aumer et al. (2016a). However, as shown in van der Kruit (1988), isothermality is not necessary for constructing an equilibrium distribution. They show that for a self gravitating disc whose vertical density distribution is exponential, the vertical dispersion is found to increase with |z|.

Dependence of dispersions on [α/Fe] abundance
Contrary to claims by Mackereth et al. (2019) that the velocity dispersion properties of the high-[α/Fe] population is different from that of the low-[α/Fe], we find very little difference between the two populations. We demonstrate this in Figure 15, where we plot the dispersion as a function of various different independent variables. The dashed lines show the best fit relation obeyed by all stars, which can be considered as the relationship obeyed by low-[α/Fe] stars, as the sample of all stars is dominated by them. The high-[α/Fe] stars are found to closely follow the dashed lines. A constant shift of about 10% can be seen between the dashed and colored lines, indicating that the overall normalization may be slightly different, but the profile shapes are very similar. For old stars, the high-[α/Fe] AVR seems to flatten more strongly than for the low-[α/Fe] AVR. We note that this effect is only prominent for samples other than GALAH-MSTO; which also happen to have larger age uncertainties as compared to the GALAH-MSTO sample. Given that high-[α/Fe] stars are expected to lie in a narrow range in age, large uncertainties in age can easily flatten the AVR. Hence, large age uncertainties seem to be the most likely reason for the apparent flattening of the AVR. 5.7. Testing the accuracy of the asteroseismic ages Giants are intrinsically bright and hence for a given apparent magnitude limit they can probe a much larger Galactic volume as compared to MSTO stars. However, it is difficult to estimate the ages of giants from purely spectroscopic parameters. Over the past decades, asteroseismology has attempted to break this barrier, backed by very precise time series photometry from space missions like CoRot, Kepler, and K2. However, measuring ages of a large numbers of stars often requires the use of asteroseismic scaling relations, which are empirical. Testing the accuracy of these scaling relations is complicated by the fact that it is difficult to get independent and precise measurements of mass or age of giants. A few techniques that have been used to verify the asteroseismic ages are, to assume that metal poor stars ([Fe/H]< −1) are older than 10 Gyr (Epstein et al. 2014), to use eclipsing binaries for estimating the mass (Gaulme et al. 2013;Brogaard et al. 2018), or to use cluster members to estimate the age Brogaard et al. (2012). While these studies suggest that the asteroseismic scaling relations overestimate masses by about 10%, they are severely hampered by small number statistics.
An indirect means to verify asteroseismic ages is to rely on ensemble statistics -for example, by comparing the mass distribution of stars against predictions of populationsynthesis-based models of the Galaxy. Sharma et al. (2016) and Sharma et al. (2017) using Kepler data suggested that the asteroseismic scaling-based masses were overestimated by about 10% compared to model predictions. However, a followup study by Sharma et al. (2019) using data from both Kepler and K2 showed that much of the tension between observations and predictions is reduced after updating the metallicity of the thick disc to recent iron abundance measurements, and additionally taking α-element abundance into account. However, we still do not have a Galactic model with all its free parameters tightly constrained. In fact certain parameters are degenerate. Hence, it is useful to look for alternative methods to verify the asteroseismic ages.
Here we provide another method to verify asteroseismic ages based on ensemble statistics -comparing the velocity dispersion of stars conditional on age, metallicity, angular momentum and distance from the plane. The underlying principle is that the velocity dispersion is a global Galactic property. Hence, groups of stars having same age, metallic-ity, angular momentum and distance from the plane, should have same velocity dispersion, irrespective of their stellar type, target selection and age estimation technique. Using the above method, we find that the conditional velocity dispersion of the GALAH-RG-K2 asteroseismic sample is in agreement with that of the GALAH-MSTO sample. The APOGEE-RG-KEPLER sample was found to have an offset with respect to the GALAH-MSTO sample. However, part of this offest is due to APOGEE iron abundance being higher by about 0.1 dex than GALAH. To account for the rest of the observed offset, the APOGEE-RG-KEPLER asteroseismic ages had to be increased by about 10%. This bias is consistent, both in direction and amount, with earlier analysis that compared the mass distribution of the Kepler sample with predictions from stellar-population-synthesis based models .
We now demonstrate that a 10% systematic in age can easily stem from inaccuracies in measurement of average seismic parameters ∆ν or ν max . Pinsonneault et al. (2018) had shown that different methods for measuring ν max and ∆ν can have systematics of up to a few percent. Given that the K2 light curve is much shorter (3 months as compared to 4 years) and is more noisy, we can also expect biases of upto a few percent in the seismic parameters estimated from them. We show below that even a 1% change in either ∆ν or ν max can lead to a change of 10% in age. The age of a red giant star is primarily determined by the time it spends on the main sequence and is roughly τ MS ∝ M/L(M) ∝ M −3.8 for stars with M < 2M (Binney & Merrifield 1998). According to asteroseismic scaling relations M ∝ ν 3 max /∆ν 4 , implying that the age depends on ν max and ∆ν with a power greater than 10.

SUMMARY AND CONCLUSIONS
We have explored the fundamental relations governing the radial and the vertical velocity dispersions of stars in the Milky Way and discussed the dynamical processes that might be responsible for them. For the first time, we present the joint dependence of the vertical and radial velocity dispersions on age, angular momentum, metallicity and distance from the plane. We compare and contrast results from three different spectroscopic surveys, (GALAH, LAMOST, and APOGEE) and three different stellar types (MSTO, asteroseismic giant, and RGB stars).
Vertical and radial velocity dispersions depend upon at least 4 independent variables, and these are age, angular momentum, metallicity and distance from the plane. The joint dependence is well approximated by a separable functional form that is a product of univariate functions, with each function corresponding to one independent variable. In other words, the dependence of the dispersions on each independent variable is almost independent of the other variables. Dispersion  The velocity dispersions increase with age following a power law, with exponent β z = 0.441 ± 0.007 for σ z and β R = 0.251 ± 0.006 for σ R . These exponents are in good agreement with idealized simulations of Aumer et al. (2016b) where the disc heating is due to scattering by bar, spiral arms and GMCs.
The velocity dispersions show a non-monotonic behaviour with L z . They decrease with L z until about solar angular momentum, thereafter, σ R flattens, while σ z increases. The flat-tening at large L z could be due to a non-zero floor on the intrinsic birth dispersion of stars. However, the cause for the rise of σ z at large L z is not fully understood. Idealized simulations having a bar, spiral structure and GMCs do not show such an effect. However, cosmological simulations by Grand et al. (2016) do show such an effect, where the heating is primarily attributed to a bar and orbiting satellites. However, other factors, like warps, the infall of misaligned gas and re-orientation of the disc, can also be responsible for the above effect.
The velocity dispersions decrease almost linearly with metallicity, or in other words the velocity dispersion increases with birth radius. We show that this can be explained by the conservation of vertical action in stars undergoing radial migration. However, for this to work, stars migrating outwards from the inner regions should be preferentially of low velocity dispersion -the so called "provenance bias" as discussed by Vera-Ciro et al. (2014) and also reported by others (Grand et al. 2016;Daniel & Wyse 2018).
The velocity dispersions increase almost linearly with distance from the plane. This effect is more prominent for younger stars. Additionally, the effect is stronger for σ z than for σ R . This agrees with findings of Mackereth et al. (2019) using APOGEE giants. Spiral arms are responsible for in-plane scattering, while GMCs are though to redirect the planar motion into the vertical direction (Jenkins & Binney 1990). As suggested by Mackereth et al. (2019), the longer time scale associated with GMC heating as compared to spiral heating could be responsible for the observed nonisothermality.
A particularly useful aspect of identifying the set of independent variables that govern the velocity dispersion, is that if the dispersion is characterized in terms of these variables then it is almost independent of the target selection function. This provides a means to not only compare results from different observational data sets (e.g to test systematics in spectroscopic stellar parameters between different surveys, or systematics between different age estimation techniques), but also to compare observational results with theoretical predictions. We take advantage of this fact to show that GALAH and LAMOST results are in agreement with each other and that results from different stellar types (MSTO and giant stars) are also in agreement with each other. The ages of giant stars have been estimated using the asteroseismic scaling relations either directly (for K2 stars) or indirectly (for LAMOST RGBs). It is difficult to verify the asteroseismic scaling relations due to a shortage of independent estimates of stellar mass or age. In this sense, we provide a new technique based on ensemble statistics to verify the accuracy of the asteroseismic ages.
The velocity dispersion of the APOGEE data set of asteroseismic giants from Kepler was found to be systematically different from our derived relations. We identify two possible reasons for this. First, the metallicity of APOGEE giants is systematically lower by about 0.1 dex with respect GALAH and LAMOST. Second, it is possible that the average asteroseismic parameters derived from Kepler data have some systematics with respect to those derived from the K2 data, given that the light curves from K2 are much shorter (3 months as compared to 4 years) and noisier.
Finally, we find that all stars, irrespective of them being old or having high-[α/Fe], follow the same relations for velocity dispersion. In other words, no special provision is needed to accommodate the thick disc stars. The AVR of stars in the solar neighborhood does show a break from a pure power law for stars older than 8 Gyr. However, when the angular momentum and metallicity of these stars is taken into account no such break is seen. The apparent break is due to older stars having systematically lower angular momentum.
SS is funded by a Senior Fellowship (University of Sydney), an ASTRO-3D Research Fellowship and JBH's Laureate Fellowship from the Australian Research Council (ARC). JBH's research team is supported by an ARC Laureate Fellowship (FL140100278) and funds from ASTRO-3D. MJH is supported by an ASTRO-3D 4-year Research Fellowship. DS is the recipient of an ARC Future Fellowship (project number FT1400147). SB and KL acknowledge funds from the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award endowed by the Federal Ministry of Education and Research. KL acknowledges funds from the Swedish Research Council (Grant 2015-00415_3) and Marie Sklodowska Curie Actions (Cofund Project INCA 600398). JK, KC and TZ acknowledges financial support from the Slovenian Research Agency (research core funding No. P1-0188). DMN was supported by the Allan C. and Dorothy H. Davis Fellowship. JZ acknowledges support from NASA grants 80NSSC18K0391 and NNX17AJ40G. DH acknowledges support from the Alfred P. Sloan Foundation and the National Aeronautics and Space Administration (80NSSC19K0108). The GALAH Survey is supported by the ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project CE170100013. This work has made use of data acquired through the Australian Astronomical Observatory, under programs: GALAH, TESS-HERMES and K2-HERMES. We acknowledge the traditional owners of the land on which the AAT stands, the Gamilaraay people, and pay our respects to elders past and present. 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 made use of data from SDSS-III. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/. This work has made use of Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope LAMOST) which is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.