Dynamical heating across the Milky Way disc using APOGEE and $\it{Gaia}$

TL;DR

This study maps the kinematic evolution of the Milky Way's disc across various radii and ages, revealing distinct heating histories for different stellar populations using APOGEE and Gaia data.

Contribution

It provides the first detailed measurements of age-velocity-dispersion relations across the Milky Way disc for mono-age, mono-metallicity populations, highlighting differences between high and low alpha sequences.

Findings

High alpha disc has flat AVRs and constant velocity dispersion ratio.

Outer disc populations show flatter radial AVRs than inner regions.

Outer disc stars likely influenced by satellite perturbations and disc warping.

Abstract

The kinematics of the Milky Way disc as a function of age are well measured at the solar radius, but have not been studied over a wider range of Galactocentric radii. Here, we measure the kinematics of mono-age, mono-$\mathrm{[Fe/H]}$ populations in the low and high $\mathrm{[\alpha/Fe]}$ discs between $4 \lesssim R \lesssim 13$ kpc and $|z| \lesssim 2$ kpc using 65,719 stars in common between APOGEE DR14 and $\it{Gaia}$ DR2 for which we estimate ages using a Bayesian neural network model trained on asteroseismic ages. We determine the vertical and radial velocity dispersions, finding that the low and high $\mathrm{[\alpha/Fe]}$ discs display markedly different age--velocity-dispersion relations (AVRs) and shapes $\sigma_z/\sigma_R$. The high $\mathrm{[\alpha/Fe]}$ disc has roughly flat AVRs and constant $\sigma_z/\sigma_R = 0.64\pm 0.04$, whereas the low $\mathrm{[\alpha/Fe]}$ disc has large variations in this ratio which positively correlate with the mean orbital radius of the population at fixed age. The high $\mathrm{[\alpha/Fe]}$ disc component's flat AVRs and constant $\sigma_z/\sigma_R$ clearly indicates an entirely different heating history. Outer disc populations also have flatter radial AVRs than those in the inner disc, likely due to the waning effect of spiral arms. Our detailed measurements of AVRs and $\sigma_z/\sigma_R$ across the disc indicate that low $\mathrm{[\alpha/Fe]}$, inner disc ($R \lesssim 10\,\mathrm{kpc}$) stellar populations are likely dynamically heated by both giant molecular clouds and spiral arms, while the observed trends for outer disc populations require a significant contribution from another heating mechanism such as satellite perturbations. We also find that outer disc populations have slightly positive mean vertical and radial velocities, likely because they are part of the warped disc.

Figures

Equations

$$L_{\mathrm{fakemag}}$$

$$=\varpi 10^{\frac{1}{5}{K}_{s,0}},$$

$${\sigma }_{[R,z]}(R,z)$$

$$\times \exp \left(\frac{-(R-R_0)}{h_{{\sigma }_{[\mathrm{R},\mathrm{z}]}}}\right)+{\delta }_{v_{[R,z]},i}^2$$

$${\sigma }_{[R,z]}(z,R_0)=a_{[\mathrm{R},\mathrm{z}]}(z-z_{1/2}{)}^2+b_{[\mathrm{R},\mathrm{z}]}(z-z_{1/2})+{\sigma }_{[\mathrm{R},\mathrm{z}]}(R_0,z_{1/2})$$

$$\tan (\alpha )={\alpha }_0+{\alpha }_1\frac{z}{R},$$

$${\sigma }_{Rz}^2(R,z)=\left({\sigma }_R^2(R,z)-{\sigma }_z^2(R,z)\right)\frac{\tan \alpha }{(1-{\tan }^2\alpha )}+{\delta }_{v_{Rz},i}^2$$

$$\Delta =\left[\begin{array}{cc}{\sigma }_R^2(R,z)& {\sigma }_{Rz}^2(R,z)\\ {\sigma }_{Rz}^2(R,z)& {\sigma }_z^2(R,z)\end{array}\right].$$

$$\ln \mathcal{L}(O|R,z,v_R,v_z)=\sum _i\ln \left[\frac{(1-ϵ)}{2\pi |\Delta {|}^{1/2}}\exp \left(-\frac{1}{2}{\vec{v}}_i{\Delta }^{-1}{\vec{v}}_i\right)+ϵp_{backgr.}(z,R)\right]$$

$$R_{\mathrm{mean}}=\frac{(r_{\mathrm{peri}}+r_{\mathrm{ap}})}{2},$$

$$\Sigma (R)\propto \frac{{\sigma }_z^2}{h_z}$$

Full text

††pubyear: 2018††pagerange: Dynamical heating across the Milky Way disc using APOGEE and Gaia–B

1 Introduction

The present day kinematics of stars in the Milky Way are immutably tied to the dynamical history of the Galaxy. As a result, the measurement of the kinematics of the stellar disc as a function of age or element abundances offers tight constraints on models for the formation and evolution of the Galaxy (Rix & Bovy, 2013). In the new era of Gaia (Gaia Collaboration et al., 2016), our access to the kinematics of stars has increased by orders of magnitude. In addition, the recent advent of large-scale spectroscopic surveys (e.g. Gaia-ESO, Gilmore et al. 2012; APOGEE, Majewski et al. 2017; GALAH Martell et al. 2017) has delivered stellar spectra, and therefore element abundance measurements, on an unprecedented scale. However, the measurement of accurate stellar ages, especially those of red-giant stars, has been problematic (see Soderblom, 2010, and references therein). On this front, recent advances, combining spectroscopic surveys with high-quality asteroseismic data (e.g. Pinsonneault et al., 2014, 2018), have lead to great improvements in age measurements (e.g. Martig et al., 2016; Silva Aguirre et al., 2018), leading Galactic astrophysics into a truly multi-dimensional mode of operation.

Pioneering work which studied the variation of stellar velocities with indicators of their age demonstrated that there exists a positive correlation between age and velocity dispersion in the solar vicinity (Strömberg, 1946; Roman, 1950a, b). This property of the Galactic disc has since been well measured and characterised in the solar neighbourhood (e.g. Wielen, 1977; Seabroke & Gilmore, 2007; Soubiran et al., 2008; Casagrande et al., 2011). Observational limitations have meant that this relationship is not yet well constrained throughout the disc which would allow a detailed assessment of its dynamical history and the heating processes which have shaped it.

Barbanis & Woltjer (1967) were among the first to show that transient spiral arms could contribute to disc heating, while Spitzer & Schwarzschild (1951, 1953) demonstrated similar effects from small scale irregularities in the disc potential, such as those caused by Giant Molecular Clouds (GMCs). Jenkins & Binney (1990) considered the combined effect of both these heating agents, finding that increasingly prominent spiral perturbations (over those from GMCs) tend to flatten the axis ratio of velocity dispersions (e.g. decreasing $\sigma_{z}/\sigma_{R}$) over time. However, these models all struggled to reconcile the predicted slope of the vertical AVR (usually modelled as a power law $\sigma_{z}\propto\tau^{\beta_{z}}$, where models predict $\beta_{z}\sim 0.25$, e.g. Lacey, 1984; Hänninen & Flynn, 2002) with that observed at the solar vicinity (where $\beta_{z}$ is closer to $\sim 0.5$, e.g. Wielen, 1977; Seabroke & Gilmore, 2007; Soubiran et al., 2008). Increased heating from GMCs early in the history of the disc can be invoked to overcome these discrpeancies (e.g. Aumer et al., 2016b, a), and a recent study of the vertical actions of APOGEE-Gaia Red Clump giants as a function of age in the low $\mathrm{[\alpha/Fe]}$ disc seemed to confirm this notion (Ting & Rix, 2018).

High resolution cosmological zoom simulations have also revealed that other non-axisymmetric features such as bars may be a dominant heating agent, alongside (rare) perturbations from relatively massive ($M\gtrsim 10^{10}\ \mathrm{M_{\odot}}$) satellite mergers (Grand et al., 2016). This is particularly prescient given that it has become apparent that the Milky Way likely underwent such a merger $\sim 10$ ago (e.g. Belokurov et al., 2018; Deason et al., 2018; Helmi et al., 2018; Kruijssen et al., 2018; Mackereth et al., 2018a). The effect of this accretion event on the dynamics of the disc itself is yet to be fully considered.

Understanding the effect of different heating agents is of course instructive, but any constraints therein on galaxy formation are limited without consideration of the velocity dispersion that the stars formed with. It is clear, from cosmological simulations at least, that discs do not tend to form thin and cool, but instead follow an "upside-down" formation (e.g. Brook et al., 2004; Brook et al., 2012; Stinson et al., 2013a; Bird et al., 2013; Stinson et al., 2013b; Grand et al., 2017; Ma et al., 2017; Navarro et al., 2017), where the stars form with the same velocity dispersion as the settling gas discs. Disentangling the combined outcome of the equilibrating gas disc and the subsequent heating of its stellar populations using observational data (e.g. Leaman et al., 2017) will lead to the greatest insights into the evolution of the Milky Way.

Whilst it has been informative to look at how the disc kinematics change with either age or abundances in isolation, it has become clear that stellar ages and element abundances in the solar neighbourhood (and beyond) do not simply correlate with one another, with perhaps the clearest example of this being $\mathrm{[\alpha/Fe]}$-$\mathrm{[Fe/H]}$ as a function of age (e.g. Martig et al., 2016; Ness et al., 2016). The discovery of the separated high and low $\mathrm{[\alpha/Fe]}$ disc components in the solar vicinity (e.g. Fuhrmann, 1998; Prochaska et al., 2000), the subsequent mapping of them through the galaxy (Adibekyan et al., 2012; Anders et al., 2014; Nidever et al., 2014; Hayden et al., 2015), and the characterisation of their structure (Bovy et al., 2012a; Bovy et al., 2012b; Haywood et al., 2013; Bovy et al., 2016), has lead to a new view of our Galaxy. In particular, it seems that the commonly assumed linkage between the high and low $\mathrm{[\alpha/Fe]}$ disc and the thick and thin disc is not as clear cut (Bovy et al., 2012a). In fact, any structural dichotomy between the high and low $\mathrm{[\alpha/Fe]}$ discs is more apparent radially, with the high $\mathrm{[\alpha/Fe]}$ disc being centrally concentrated, and the low $\mathrm{[\alpha/Fe]}$ populations occupying donut shaped annuli (Bovy et al., 2016), which change shape as a function of age and metallicity (Mackereth et al., 2017). Haywood et al. (2016) present similar arguments, proposing an inner/outer disc divide to be more constraining to the evolution of the disc.

Models for the formation of the bimodality in $\mathrm{[\alpha/Fe]}$ generally predict that high $\mathrm{[\alpha/Fe]}$ stars form in rapid and intense star formation, precluding enrichment by Type Ia SNe. These high star formation efficiency environments are either brought about by rapid infall of gas (e.g. Chiappini et al., 1997, 2001; Chiappini, 2009), or are present in the innermost, high density regions of modelled discs (e.g. Schönrich & Binney, 2009a, b). Alongside such analytic models, $\mathrm{[\alpha/Fe]}$ bimodality has also been recently realised in fully self-consistent cosmological simulations. Grand et al. (2018) found that $\mathrm{[\alpha/Fe]}$ bimodality in the discs of the AURIGA zoom-in simulations (Grand et al., 2017) arises due to a double peaked star formation history, interjected by a shrinking of the gas disc. They also showed that $\mathrm{[\alpha/Fe]}$ bimodality can appear without this process, via an initial rapid burst of star formation in the inner galaxy, but noted that the bimodality was not ubiquitous among the simulations. Mackereth et al. (2018b) further showed that $\mathrm{[\alpha/Fe]}$ bimodality in the EAGLE simulations (Schaye et al., 2015; Crain et al., 2015) is very rare (occuring in $\sim 6\%$ of Milky Way mass galaxies), and driven largely by the atypical assembly history of the haloes hosting the bimodal galaxies, which accreted mass faster at earlier times (and slower at late times) than their non-bimodal counterparts. Importantly, they showed that low and high $\mathrm{[\alpha/Fe]}$ populations form chemically separated from each other and usually overlap in age, becoming co-spatial later in the history of the galaxies.

These models for the formation of the high and low $\mathrm{[\alpha/Fe]}$ components make qualitative predictions for the kinematic structure of these populations, which are becoming more easily testable as the available data becomes more sophisticated. For example, the ‘radial migration’ model of Schönrich & Binney (2009a) predicts that outwardly migrating high $\mathrm{[\alpha/Fe]}$ stars, formed in the inner disc at higher $\sigma_{z}$ act to increase the slope of the age-velocity dispersion relationship in the solar neighbourhood (relative to models where stars do not migrate). Therefore, isolating the low $\mathrm{[\alpha/Fe]}$ disc AVR should yield a reconciliation with the relation predicted by models of time-dependent scatterers ($\beta_{z}\sim 0.25$, e.g. Lacey, 1984; Hänninen & Flynn, 2002). ‘Two-infall’ type models require a very rapid infall of gas onto the proto-Milky Way to form the high $\mathrm{[\alpha/Fe]}$ stars, which may correspond to the clumpy galaxies observed at high redshift (e.g. Noguchi 1998, and see also the recent preprint by Clarke et al. 2019), and rapid merging at similar times seen in simulations (e.g. Brook et al., 2004; Bournaud et al., 2009), consistent with the Milky Way-like galaxies from Mackereth et al. (2018b). Rapid, early accretion and subsequent intense star formation would likely produce a centrally concentrated and very kinematically hot high $\mathrm{[\alpha/Fe]}$ population. Such models require that the high and low $\mathrm{[\alpha/Fe]}$ populations become co-spatial, and so require that the high $\mathrm{[\alpha/Fe]}$ populations are heated radially somehow, bringing stars to the solar radius. The distinction in kinematics between radial migration and heating may be slight, but would place a very strong constraint on these models. At present, the main distinction between radial migration and heating can be made by measuring the flaring of the disc (e.g. Bovy et al., 2012b), as it is expected that the flaring resultant from each process is different (e.g. Minchev et al., 2012). Understanding whether the flaring, the disc surface density, and its kinematics are commensurate with one another would help to disentangle these processes.

Since the second Gaia data release (DR2 Gaia Collaboration et al., 2018a), there has been a rapid development in our understanding of the kinematics and dynamics of the disc. Gaia Collaboration et al. (2018b) mapped the kinematics of the largest portion of the Milky Way to date, revealing a rich structure and unveiling the clear non-axisymmetries of the disc, providing the first insight into our new view of the Galactic disc. Antoja et al. (2018) showed that Gaia DR2 data reveals a ‘phase-spiral’ feature in the solar neighbourhood, which is now widely considered to be the lasting harmonic relic of a recent perturbation by a satellite flyby (e.g. Bland-Hawthorn et al., 2018; Binney & Schönrich, 2018; Laporte et al., 2019). Such a wave-like oscillation is also apparent as an asymmetry in number counts as a function of height above the Galaxy midplane (Bennett & Bovy, 2018). It is worth noting also that the effect of phase-mixing and warping on phase space in the disc was studied prior to Gaia by a number of works (e.g. Minchev et al., 2009; Quillen et al., 2009; Gómez et al., 2012; de la Vega et al., 2015). The general picture which emerges is that the disc which is likely somewhat out of equilibrium due to recent events in its history.

As mentioned briefly above, a study of the vertical kinematics of Red Clump (RC) giant stars in common with APOGEE and Gaia DR2 recently suggested that gradual orbit scattering by small perturbations like those from GMCs may be enough to explain the observed trends between vertical actions and age - showing that a model of birth temperature scaling with star formation rate did not fit the data well (Ting & Rix, 2018). That study focused mainly on the younger, low $\mathrm{[\alpha/Fe]}$ dominated populations, whereas we extend the modelling of the vertical kinematics to the high $\mathrm{[\alpha/Fe]}$, old stellar populations, as well as examining the radial kinematics. Ting & Rix (2018) found that the slope of vertical action-age relationship (which can be understood as the adiabatically invariant counterpart to the AVR) increased with Galactocentric radius, but found this increase to also be in-line with expectations from GMC scattering.

In this paper, we present a dissection of the disc kinematics in age, $\mathrm{[Fe/H]}$ and $\mathrm{[\alpha/Fe]}$, exploiting an unprecedented, multi-dimensional data set consisting of stars with element abundances and ages based on APOGEE spectra, whose distances are estimated to high precision using neural network modelling of these spectra trained on Gaia DR2 data, which also offers high-quality proper motion information. In Section 2, we present and describe this data set, and describe the procedure by which we estimate ages for a large portion of the APOGEE catalogue, using a Bayesian Convolutional Neural Network trained on APOGEE spectra and asteroseismic data from the APOKASC catalogue. We complete that section by describing the features seen in the data in the age-$\mathrm{[\alpha/Fe]}$-$\mathrm{[Fe/H]}$ plane. In Section 3, we present the procedure used to model the velocity dispersions in the vertical and radial direction in the Galactic disc for mono-age, mono-$\mathrm{[Fe/H]}$ populations in the low and high $\mathrm{[\alpha/Fe]}$ discs. Section 4 presents the main results from modelling the velocity dispersion of the mono-age, mono-$\mathrm{[Fe/H]}$ populations. In Section 5, we assess the results in the context of previous work on constraining the effects of disc heating, and discuss the implications of our findings on models for the origin of the high and low $\mathrm{[\alpha/Fe]}$ populations. We summarise and conclude the paper in Section 6.

2 Data

We use a catalogue of stellar positions, velocities, element abundances and estimated ages, comprising a cross match of the catalogues from the fourteenth data release (DR14, Abolfathi et al., 2018) of the SDSS-IV APOGEE-2 survey, and the second data release (DR2, Gaia Collaboration et al., 2018a) of the ESA-Gaia mission. Ages are estimated from a neural network based model (described below) trained on data from the APOKASC catalogue, which contains stars observed both spectroscopically by APOGEE and asteroseismically by the Kepler mission.

2.1 The APOGEE-2 DR14 catalogue

APOGEE (Majewski et al., 2017) is a spectroscopic survey of the Milky Way in the near infra-red H-band (1.5-1.7$\mathrm{\mu}$m), which has observed over 200,000 stars at high signal-to-noise ratio (SNR > 100 pixel*-1*) and high resolution (R $\sim$ 22 500), measuring over 15 element abundances. We use the DR14 data (Abolfathi et al., 2018), which consists of a combination of stars observed between the first and second iterations of the survey, APOGEE-1 (forming part of SDSS-III; Eisenstein et al., 2011) and 2 (in SDSS-IV; Blanton et al., 2017). In this paper we refer to these surveys collectively as APOGEE. All APOGEE data products used in this paper are those output by the standard data analysis pipeline, the APOGEE Stellar Parameters and Chemical Abundances Pipeline (ASPCAP García Pérez et al., 2016), which uses a precomputed spectral library (Zamora et al., 2015), synthesised using a customised $H$-band linelist (Shetrone et al., 2015), to measure stellar parameters and element abundances. A full description and examination of the analysis pipeline is given in Holtzman et al. (2018). The individual element abundances are well tested against samples from the literature (Jönsson et al., 2018), and found to agree very well.

We refer here only to the abundances of the $\alpha$-elements included in DR14: Oxygen, Magnesium, Silicon, Sulphur, and Calcium, and the abundance of Iron. All abundances are those calculated by ASPCAP, and included in the APOGEE catalogue. We combine the $\alpha$ elements to attain the ratio of the mean $\alpha$ element abundance to that of Iron, $\mathrm{[\alpha/Fe]}$. The element abundances are determined by ASPCAP via a two-step process. The best fit stellar parameters $T_{\mathrm{eff}}$, $\log(g)$, $\nu_{\mu}$ (microturbulent velocity), [M/H], $\mathrm{[\alpha/M]}$, [C/M] and [N/M] are determined by a global fit to the grid of synthetic stellar spectra. Individual element abundances are then estimated by fitting windows in the spectrum to synthetic spectra with varying [M/H] (or $\mathrm{[\alpha/M]}$, [C/M] or [N/M], for $\alpha$-elements, Carbon and Nitrogen, respectively). For DR14, a small external calibration is applied to all the abundances, which forces the abundance ratios of solar metallicity stars in the solar vicinity to be equal to solar (Holtzman et al., 2018).

We do not use Gaia parallaxes directly as distance indicators, but instead use spectro-photometric distances obtained by training a neural network to predict the luminosity of a star from its infrared spectrum using luminosities for a training set obtained from Gaia DR2. This procedure is described in detail in Leung & Bovy (2019a), but we re-iterate it here for clarity. The neural network has a similar architecture to the ApogeeBCNN() network used in Leung & Bovy (2019b): it is composed of 2 convolutional layers and 2 dense layers with rectifier activation (which introduce non-linearity in the hidden layers and prevent non-physical negative luminosity in the output); the network is optimized with the ADAM optimizer (Kingma & Ba, 2014) and is implemented in astroNN (Leung & Bovy, 2019b), which itself relies on the TensorFlow framework (Abadi et al., 2016). To train the network in a way that allows us to take the Gaia parallax uncertainty into account, we follow Anderson et al. (2018) and define the output to be the following transformation of the absolute magnitude $M_{K_{s}}$

[TABLE]

where $\varpi$ is the parallax in units of mas and $K_{s,0}$ is the extinction-corrected apparent magnitude in the $K_{s}$ band (for all stars, we obtain extinction corrections directly from the APOGEE catalogue). We then train the neural network to predict the ’pseudo-luminosity’ $L_{\mathrm{fakemag}}$ from the continuum-normalized APOGEE spectrum using a $\chi^{2}$ objective function that takes the uncertainty in the $L_{\mathrm{fakemag}}$ of the training set (due to the uncertainty in $\varpi$), into account. When training the neural network, we adopt a $0.0562$ mas Gaia DR2 parallax offset, because this value leads to approximately unbiased luminosities at small $\varpi$ (see Leung & Bovy, 2019a, for more details on this procedure). We train the network using APOGEE DR14 (Holtzman et al., 2018) spectra that have SNR$>200$, no APOGEE_STARFLAG flags set, radial velocity scatter smaller than $1\,\mathrm{km\,s}^{-1}$, and $\varpi_{Gaia}/\sigma_{\varpi,Gaia}>1$.

Once trained, we use the network to predict luminosities $L_{\mathrm{fakemag}}$ for the entire APOGEE DR14 sample and, combined with extinction-corrected $K_{s,0}$, we obtain spectro-photometric parallaxes using Equation (2). We obtain uncertainties on the predicted $L_{\mathrm{fakemag}}$ using dropout variational inference as an approximation of a Bayesian neural network (Leung & Bovy, 2019b) and these uncertainties are propagated to the spectro-photometric parallax. The distances are well tested in Leung & Bovy (2019a) against the APOGEE red clump star (RC) catalogue (Bovy et al., 2014), alternative spectrophotometric distances from Santiago et al. (2016) and also those from (Hogg et al., 2018), and found to be in very good agreement, if not with improved precision. We also test the spectrophotometric distances against the recently published catalogue of distances of Schönrich et al. (2019). We find that the distances agree to within $\sim 10\%$, comparable with the uncertainties quoted in the Leung & Bovy (2019a) catalogue.

2.2 Gaia DR2

Gaia DR2 contains astrometric parameters for over 1.3 billion sources, collected over 22 months from July 2014. Many major improvements were made over the initial data processing, released in the DR1 and Tycho-Gaia Astrometric Solution (TGAS) catalogues, which are summarised in Gaia Collaboration et al. (2018a). In particular, Gaia DR2 is the first data release which is not tied to any external catalogue to provide proper motion measurements (DR1 was tied to the Hipparcos and Tycho-2 catalogues), and uses its own refernce frame based on quasars (Gaia-CRF-2, Gaia Collaboration et al., 2018c). Other improvements included better modelling of the spacecraft attitude, and a significant improvement to the source detection algorithm which limits the number of spurious sources, and removes many duplicate sources which were included in the DR1 source catalogue (for this reason the catalogues between data releases are treated separately). As well as removing spurious sources, these improvements have led to significant increases in the astrometric precision, affording typical residuals on the astrometric solution of $\sim 0.2$ to 0.3 mas in the middle of the magnitude range, decreasing to $\sim 2$ mas for the faintest sources.

We cross match the full APOGEE DR14 catalogue with the Gaia DR2 source catalogue, searching in a cone of radius $0.5"$ around each source for its corresponding object in the Gaia catalogue using the CDS X-match service111http://cdsxmatch.u-strasbg.fr/xmatch. We find that there are 254,789 objects in common between APOGEE and Gaia DR2. After removing stars from APOGEE which have warning or bad flags set, selecting only stars which have proper-motion measurements from Gaia, and removing duplicate entries from the APOGEE catalogue, 83,189 stars with 6D phase space information remain. In the following analysis, we only use stars with $1.8<\log(g)<3.0$ (removing dwarf and subgiant stars, and upper RGB stars) and those with $\mathrm{[Fe/H]}$ $>-0.5$, below which the APOKASC catalogue (used in the following section for training our age estimation model) has poor coverage. In addition to this, it was recently shown that extra mixing becomes significant at approximately this metallicity, which would substantially change the relationship between mass (and therefore the inferred ages) and $\mathrm{[C/N]}$ (Shetrone et al., in prep.), which is likely to be the parameter which drives the relationship between the spectroscopic data and age. This final sample contains 65,719 red giant stars, whose spectra are measured with a minimum SNR of $\sim 80$, and a median SNR of $\sim 170$. The median uncertainty on the proper motion measurements from Gaia is $\sim 0.5$ mas, corresponding to a median uncertainty (after transformation) on the Galactocentric radial and vertical velocities of $\sim 2$ and $\sim 1\ \mathrm{km\ s^{-1}}$, respectively. The median uncertainty on the distances from astroNN (described in Section 2.1) is also remarkably low, at $\sim 0.24$ kpc, better than 10% in most cases.

We transform the positions, proper motions, and radial velocities into the Galactocentric cylindrical frame, adopting the radial and vertical solar motion relative to the local standard of rest of Schönrich et al. (2010) and the tangential motion of the sun relative to the Galactic center of $245.6\ \mathrm{km\ s^{-1}}$ computed using the proper motion and distance to Sgr A* (Reid & Brunthaler, 2004; Gravity Collaboration et al., 2018). We propagate the observational uncertainties and their covariance matrix into this frame also. The data extends in Galactocentric radius roughly between $4\lesssim R\lesssim 13$ kpc, and between $-1\lesssim z\lesssim 2$ kpc above and below the plane (here, and throughout the paper, we assume $R_{0}=8$ kpc and $z_{0}=0.025$ kpc). We show the distribution of the data in $R$ and $z$ in Figure 1. This figure demonstrates the strong spatial selection biases which are present when using APOGEE data. For example, there is a large overdensity of stars above the plane, just inside the solar radius, that corresponds to the heavily-observed Kepler field. Correcting for these biases is possible, but not necessary for the kinematic modelling performed here, because neither the APOGEE nor Gaia data are kinematically biased. As we are interested only in modelling the velocity distributions as a function of position, these spatial biases do not affect our procedure.

2.3 Asteroseismic data and age estimates for APOGEE red giants

To make estimates of the ages of stars in the APOGEE-Gaia data described above, we make use of data from the APOKASC-2 catalogue (Pinsonneault et al., 2018) of asteroseismic data for stars in common between APOGEE and Kepler (Borucki et al., 2010). We use the APOKASC-2 data set to train a Bayesian Convolutional Neural Network (BCNN) model, as implemented in the astroNN python package (Leung & Bovy, 2019b), to predict stellar ages from the APOGEE spectra. Our technique for modelling ages using the BCNN is described fully in Appendix A. Briefly, we expect that the ability of the BCNN model to predict ages from the spectra to a good degree of accuracy stems from the presence of molecular bands of Carbon and Nitrogen in the APOGEE spectra. It is well documented in the literature that there is likely to be a relationship between these element abundances and the stellar mass (and therefore age, e.g. Salaris et al., 2015), and this has been exploited in previous work to estimate ages from APOGEE abundances and spectra (Martig et al., 2016; Ness et al., 2016; Das & Sanders, 2019). The BCNN uses the full information content of the APOGEE spectra to estimate the ages, and allows for the proper propagation of uncertainties into the analysis, to predict ages with good accuracy and with reasonable error estimation. Using the method, we generate a catalogue of ages with a median uncertainty between $\sim 30$ and $35\%$ accross the full range of ages. We fully discuss the limitations of the predicted ages in Appendix A, but note here for clarity that we expect that ages above $\sim 10$ Gyr are likely to be underestimated, and subject to very large errors, such that it becomes difficult to distinguish, for example, an 11 Gyr old star from one at 13 Gyr old.

2.3.1 The age-$\mathrm{[\alpha/Fe]}$-$\mathrm{[Fe/H]}$ relation of disc stars

For completeness, in Figures 2 and 3 we demonstrate the complexity of the sample in age, $\mathrm{[\alpha/Fe]}$ and $\mathrm{[Fe/H]}$ space. Figure 2 shows the $\mathrm{[\alpha/Fe]}$-$\mathrm{[Fe/H]}$ plane as a 2D histogram, coloured by the mean age of the stars in each bin. The opacity of the bins reflects the density of stars in that region of $\mathrm{[\alpha/Fe]}$-$\mathrm{[Fe/H]}$ space. The solid and dashed lines indicate the ‘by-eye’ cut used to define high and low $\mathrm{[\alpha/Fe]}$ stars later in the paper. Stars inside the dashed lines, the gap between which is equal to the median uncertainty in $\mathrm{[\alpha/Fe]}$ ($\sim 0.05$ dex), are not used in the analysis, to avoid contamination of either sample by the other. The high $\mathrm{[\alpha/Fe]}$ population is dominated by stars $\gtrsim 7$ Gyr old, whereas the low $\mathrm{[\alpha/Fe]}$ population has a wide range of ages from the youngest at $\sim 1-2$ Gyr up to $\sim 6-7$ Gyr.

The age-$\mathrm{[\alpha/Fe]}$ relation, coloured by $\mathrm{[Fe/H]}$, is shown in the top panel of Figure 3. The age-$\mathrm{[Fe/H]}$ relation, coloured by $\mathrm{[\alpha/Fe]}$, is shown in the lower panel. There are a number of features evident in these planes which warrant discussion, and provide a useful context for the further results which we present in the paper. The age-$\mathrm{[\alpha/Fe]}$ relationship at the solar radius has been shown to have very little scatter at old ages, extending to high $\mathrm{[\alpha/Fe]}$ in a tight sequence, which becomes slightly more scattered at low $\mathrm{[\alpha/Fe]}$ and young age (e.g. Haywood et al., 2013). However, this appears to be dependent on the sample and age measurement method used, as asteroseismic analysis of stars in the Kepler field (within 2 kpc of the Sun) suggests the age-$\mathrm{[\alpha/Fe]}$ relationship may actually be broader at high $\mathrm{[\alpha/Fe]}$ (Silva Aguirre et al., 2018). Our extension of the data to a much wider range of Galactocentric radii somewhat complicates this picture and provides some insight into this (as yet, unsettled) issue. The scatter in $\mathrm{[\alpha/Fe]}$ is increased at all ages, and a number of interesting features become apparent.

First, it is clear from the top panel that the high $\mathrm{[\alpha/Fe]}$ ‘sequence’, as is well characterised in the literature (e.g. Bensby et al., 2005; Hayden et al., 2015; Nidever et al., 2014), which is apparent here extending between $6\lesssim\mathrm{age}\lesssim 10$ Gyr, does overlap slightly in age with the low $\mathrm{[\alpha/Fe]}$ stars with ages $\sim 6$ Gyr old. The spread in $\mathrm{[\alpha/Fe]}$ of old stars ($\gtrsim 6$ Gyr) is significant, and bridges the gap between the low and high $\mathrm{[\alpha/Fe]}$ stars, with intermediate stars assuming roughly solar values of $\mathrm{[Fe/H]}$. This high $\mathrm{[\alpha/Fe]}$ sequence feature manifests itself in the lower panel as the old ‘branch’ that extends from $\mathrm{(age,[Fe/H])}\sim(8,-0.4)$ to $\mathrm{(age,[Fe/H])}\sim(7,0.0)$, joining the large ‘cloud’ of low $\mathrm{[\alpha/Fe]}$ stars at the most metal-rich end (the low $\mathrm{[\alpha/Fe]}$, $\mathrm{[Fe/H]}$ rich, old blob at the lower-right of the low $\mathrm{[\alpha/Fe]}$ population in the top panel). We note here that the spread in age of the oldest high $\mathrm{[\alpha/Fe]}$ stars is likely to be inflated by the age uncertainties and therefore, any trends discussed among these populations is subject to this significant caveat. Our results do appear, however, to be qualitatively similar to those of Silva Aguirre et al. (2018), and show a spread in age at fixed $\mathrm{[\alpha/Fe]}$ for higher $\mathrm{[\alpha/Fe]}$ stars.

The low $\mathrm{[\alpha/Fe]}$ stars which are identifiable separately from the high $\mathrm{[\alpha/Fe]}$ ‘sequence’ in the top panel demonstrate a number of interesting features in both projections. In age-$\mathrm{[\alpha/Fe]}$, it is clear that there is some scatter in $\mathrm{[\alpha/Fe]}$ at fixed age, most evident at $\mathrm{age}\sim 4$ Gyr, which is significant over the characteristic uncertainties on the $\mathrm{[\alpha/Fe]}$ abundances ($\sim 0.05$ dex). This spread correlates with $\mathrm{[Fe/H]}$, such that the stars with slightly enhanced $\mathrm{[\alpha/Fe]}$ have the lowest values of $\mathrm{[Fe/H]}$. These $\mathrm{[\alpha/Fe]}$ enhanced low $\mathrm{[\alpha/Fe]}$ population stars are evident in the lower panel at $\mathrm{(age,[Fe/H])}\sim(4,-0.4)$, with a similar $\mathrm{[Fe/H]}$ as the oldest stars in the high $\mathrm{[\alpha/Fe]}$ population. The oldest stars in the low $\mathrm{[\alpha/Fe]}$ population exhibit the highest $\mathrm{[Fe/H]}$ values. These old, low $\mathrm{[\alpha/Fe]}$ population stars are clear in the lower panel at $\mathrm{(age,[Fe/H])}\sim(5,0.4)$. Understanding the nature of these stellar populations is likely to shed some light on the history of the disc.

In the following analysis, we aim to put these features into the context of the evolution of the Milky Way disc, using their kinematics to inform their origins and subsequent evolution. We will show that the different populations present in these planes are the result of a complex history of disc formation, which we can begin to disentangle and re-assess using the multi-dimensional data set at hand. A future paper will complement the work here by studying the detailed orbital properties of these stars.

3 Modelling the Galactic kinematics of mono-age, mono-$\mathrm{[Fe/H]}$ populations

We model the radial and vertical velocity dispersion profiles of mono-age, mono-$\mathrm{[Fe/H]}$ populations as multivariate Gaussian distributions, explicitly fitting for the covariance terms, allowing us not only to assess the variation of the velocity dispersion as a function of age and $\mathrm{[Fe/H]}$, but also to study variations in the tilt of the velocity ellipsoid with these parameters. This methodology also allows us to account for outliers and readily convolve the model with the observational uncertainties which, for Gaia, can be strongly correlated with one another.

We model the radial and vertical velocity dispersion with smoothly varying functions of $R$ and $z$, allowing them to vary exponentially in the $R$ direction, and following a quadratic relation in the $z$ direction (as in Bovy et al. 2012c) such that

[TABLE]

where $a_{[R,z]}$ and $b_{[R,z]}$ are, respectively, linear and quadratic coefficients in the quadratic function

[TABLE]

which is centered on the median $z$ value, $z_{1/2}$. Therefore, $\sigma_{[R,z]}(R_{0},z_{1/2})$ gives the velocity dispersion at $z_{1/2}$ and $R_{0}$ (which here we define as the solar radius, assumed to be 8 kpc). The $R$ dependence is presumed to be exponential, with a scale length $h_{\sigma_{[R,z]}}$. The addition of $\delta^{2}_{v_{[R,z]},i}$ accounts for the convolution of the model with the error on the velocity measurements. These assumptions on the $R$ and $z$ dependence of the velocity dispersion are not intended to be prescriptive of the actual underlying dependence, for which we have little prior expectation, but are intended to approximate it with smooth functions.

We assume that the tilt angle $\alpha$ is a function of $R$ and $z$ such that

[TABLE]

which includes the cases of a velocity ellipsoid that is always aligned with the Galactocentric cylindrical coordinate system ($\alpha_{0}=\alpha_{1}=0$) or the Galactocentric spherical coordinate system ($\alpha_{0}=0$ and $\alpha_{1}=1$). The covariance between the radial and velocity dispersion is then

[TABLE]

where we again convolve the model with the covariance between the radial and vertical velocity $\delta^{2}_{v_{R,z},i}$. This then allows the construction of the covariance matrix which describes the velocity ellipsoid in the $R$ and $z$ directions

[TABLE]

The likelihood of the parameters $O=[\sigma_{[R,z]}(R_{0},z_{1/2}),a_{[R,z]},b_{[R,z]},h_{\sigma_{[R,z]}},\alpha_{0},\alpha_{1},v_{R,0},v_{z,0},\epsilon]$ given the data for a given age-$\mathrm{[Fe/H]}$ bin can then be expressed as

[TABLE]

where $p_{\mathrm{backgr.}}(R,z)$ describes the velocity distribution of an interloper model that is contributed by a fraction of stars $\epsilon$. We use an interloper model which is also a normalised multivariate Gaussian, but with no covariance terms and with $\sigma_{R}=\sigma_{z}=100\,\mathrm{km\ s^{-1}}$. This interloper model is also convolved with the observational uncertainty for each star. We find that for all mono-age, mono-$\mathrm{[Fe/H]}$ bins, $\epsilon$ is less than a few percent. The velocity distribution is a multivariate Gaussian centered on $v_{[R,z],0}$, where $v_{[R,z],0}$ are allowed to vary freely. We find that allowing this variation in $v_{[R,z],0}$ provides a better fit to the data, and also reveals information regarding the non-axisymmetric velocity structure of the disc.

We minimise the negative log-likelihood using a downhill simplex algorithm (Nelder & Mead, 1965), and use this optimal solution to initiate an Markov Chain Monte Carlo (MCMC) sampling of the posterior PDF of the parameters O using an affine-invariant ensemble MCMC sampler (Goodman & Weare, 2010) as implemented in the python package emcee (Foreman-Mackey et al., 2013). In the following, all reported parameter values and uncertainties are the median and standard deviation of the resulting MCMC chain. For all MCMC fits, we use 100 ‘walkers’ over 500 iterations, generating 50,000 samples, of which we remove 300 from each chain, leaving 20,000 clean samples of the posterior PDF in each age-$\mathrm{[Fe/H]}$ bin. Cutting many samples allows the MCMC chains to be very well ‘burnt-in’, ensuring an accurate sampling of the underlying PDF.

We perform this fitting procedure across 8 bins in age with width $1.5$ Gyr between $1<\mathrm{age}<13$ Gyr, and 10 bins in $\mathrm{[Fe/H]}$ of width 0.1 dex, between $-0.5<\mathrm{[Fe/H]}<0.5$ dex (the adopted age-$\mathrm{[Fe/H]}$ bins are shown in Figure 3), for stars in the low and high $\mathrm{[\alpha/Fe]}$ populations, as defined in Figure 2. Modelling is only performed for bins with $N>200$ stars to ensure good sampling of the underlying velocity distribution. This results in 35 populations modelled in the low $\mathrm{[\alpha/Fe]}$ population and 13 in the high $\mathrm{[\alpha/Fe]}$. In both the low and high $\mathrm{[\alpha/Fe]}$ populations, the oldest bin ($11.5<\mathrm{age}<13$ Gyr) is never sufficiently filled at any $\mathrm{[Fe/H]}$ and so we disregard it in the further discussion, reminding the reader that it is likely that the ages of the oldest stars in our sample ($\sim 10-11$ Gyr) are likely underestimated. The median number of stars in each bin for the low $\mathrm{[\alpha/Fe]}$ populations is 1078, with a minimum of 204, and a maximum of 3916. The median number of stars in high $\mathrm{[\alpha/Fe]}$ bins is lower, at 503, with a maximum of 771 and a minimum of 227. We test the robustness of the fits inside the low N bins by fitting stars in a test age-$\mathrm{[Fe/H]}$ bin (in the low $\mathrm{[\alpha/Fe]}$ population, at $2.5<\mathrm{age}<4.0\ \mathrm{Gyr}$ and $-0.2<\mathrm{[Fe/H]}<-0.1$ dex), using 1500,1000,500 and 200 random samples of stars from the bin to perform the fit. We find that the best fit values of all of the parameters do not vary, within their uncertainties, due to the lower number statistics. Thus, we can confidently use the bins with $N>200$ stars.

We demonstrate the resulting posterior PDFs for parameters $O$ for an example bin in age-$\mathrm{[Fe/H]}$ in Appendix B.

4 Results

4.1 Velocity dispersion profiles

We first examine the velocity dispersion profiles for mono-age mono-$\mathrm{[Fe/H]}$ populations. The radial and vertical velocity dispersion profiles as a function of height above the midplane $|z|$ at $R_{0}$ for each mono-age mono-$\mathrm{[Fe/H]}$ population are shown in Figure 4. We show the variation of $\sigma_{z}$ (left) and $\sigma_{R}$ (right) for each mono-age mono-$\mathrm{[Fe/H]}$ population between the 5th and 95th percentile of the $|z|$ distances for that bin. The median height $z_{1/2}$ is indicated for each population by the scatter points. Although these profiles were not constrained to be flat, and $\sigma_{R}$ and $\sigma_{z}$ could vary smoothly with $|z|$, we find that a generally flat profile is fit in almost all cases, within the uncertainties, meaning that for most bins, the parameters $a_{[R,z]}$ and $b_{[R,z]}$ are fit to be very small, and with low uncertainty. This means that these populations are well approximated as isothermal discs. While this was also shown to be true in $\sigma_{z}$ for mono-abundance populations (MAPs) in SDSS/SEGUE by Bovy et al. (2012c), here, we extend that finding to show that mono-age, mono-$\mathrm{[Fe/H]}$ populations (for which MAPs are near analogues, aside from a few caveats, see Minchev et al., 2017) show very little variation in $\sigma_{z}$ and $\sigma_{R}$ as a function of $|z|$. It is worth noting at this point that Bovy et al. (2012c) used the isothermality of MAPs to demonstrate the precision of the SEGUE element abundances. While we do not rigourously perform a similar test here, our finding that most mono-age bins demonstrate isothermality over a large range in $|z|$ is likewise indicative that the precision in our age measurements is likely similar to the size of bins that we use (given the assumptions made in Bovy et al., 2012c).

An interesting exception to the isothermality is apparent in the younger populations shown in Figure 4. Both $\sigma_{z}$ and $\sigma_{R}$ in these populations does appear to show some variation as a function of $|z|$, such that the velocity dispersion in both directions is slightly increased at higher $|z|$. This is unlikely to be due to uncertainties in the age measurement, as the ages are well constrained at the youngest ages. A slightly increasing $\sigma_{z}$ as a function of $|z|$ was also seen in the lowest $\mathrm{[\alpha/Fe]}$ MAPs in Bovy et al. (2012c), which are likely roughly equivalent to these young populations. We return to this in Section 5.4.

Further to the generally isothermal nature of the mono-age, mono-$\mathrm{[Fe/H]}$ populations, we also find that the dependence of $\sigma_{R}$ and $\sigma_{z}$ on $R$ is generally very weak. The scale lengths of the assumed exponential dependencies are long, with median values across all age and $\mathrm{[Fe/H]}$ bins of $h_{\sigma_{R}}=15_{-4}^{+11}$ kpc and $h_{\sigma_{z}}=16_{-5}^{+19}$ kpc. However, the constraints on the scale lengths are relatively poor in comparison to the stronger constraints placed on the other parameters of the model. This is likely due to the fact that in any given age-$\mathrm{[Fe/H]}$ bin, the radial extent of the data is relatively small. Regardless of the poor constraints on the scale length of the $R$ dependence, we note that although we allowed the data to be fit with either increasing or decreasing $\sigma$ as a function of $R$, we find that it is best fit by a slowly decreasing profile in all cases.

Finally, we find that the data are fit in most of the age-$\mathrm{[Fe/H]}$ bins by models where the tilt angle of the velocity ellipsoid $\alpha$ is aligned between the Galactocentric spherical and cylindrical coordinate system (i.e. we find that $\alpha_{1}$ lies between 0 and 1, with $\alpha_{1}=0$ corresponding to cylindrical and $\alpha_{1}=1$ corresponding to spherical), with a large uncertainty, such that most bins are consistent with $\alpha_{1}=0$. We compared our fits to the work of Büdenbender et al. (2015), who found that stars between $0.5\lesssim|z|\lesssim 2.0$ kpc were well fit by a relationship close to a spherical alignment, finding that $\alpha=0.9\arctan(|z|/R_{\odot})-0.01$. On directly comparing our fits with their results (plotting $\alpha(z_{1/2})$ against $z_{1/2}$), we find that the mono-age, mono-$\mathrm{[Fe/H]}$ populations at $z_{1/2}\gtrsim 0.5$ kpc agree well with Büdenbender et al. (2015), albeit with very large uncertainties. Younger bins, at lower $z_{1/2}$, also appear to lie on the Büdenbender et al. (2015) relation. Binney et al. (2014) found similarly consistent results for the tilt angle in the RAVE data.

4.2 Age - velocity dispersion relations in the high and low $\mathrm{[\alpha/Fe]}$ discs

We now turn to examining the trends between $\sigma_{R}$ and $\sigma_{z}$ with age. Quoted and displayed values of $\sigma_{[R,z]}$ in the following are those evaluated at $R_{0}$ and the median height $z_{1/2}(\mathrm{age,[Fe/H]})$ which, due to the weak dependence of $\sigma_{[R,z]}$ with $z$, essentially corresponds to $\sigma_{[R,z]}(R_{0})$. Figure 5 shows $\sigma_{z}$ (left) and $\sigma_{R}$ (right) against age for the mono-age, mono-$\mathrm{[Fe/H]}$ populations in the low and high $\mathrm{[\alpha/Fe]}$ populations. At each age, we show all mono-$\mathrm{[Fe/H]}$ bins in which there were more than 200 stars. We colour the points by the median of the mean orbital radii of the stars in each bin $\langle R_{\mathrm{mean}}\rangle$, calculated as

[TABLE]

where $r_{\mathrm{peri}}$ and $r_{\mathrm{ap}}$ are the peri- and apocentre radius of the orbits, which are determined using the orbital parameter estimation method of Mackereth & Bovy (2018), as implemented in the galpy python package for galactic dynamics (Bovy, 2015). The mono-age, mono-$\mathrm{[Fe/H]}$ populations in the high $\mathrm{[\alpha/Fe]}$ disc are displayed as open points, whereas the low $\mathrm{[\alpha/Fe]}$ populations are solid points.

The velocity dispersions $\sigma_{z}$ and $\sigma_{R}$ both increase roughly monotonically as a function of age, with $\sigma_{z}$ going from $\sim 10$ to $\sim 40\,\mathrm{km\ s^{-1}}$ and $\sigma_{R}$ increasing from $\sim 30$ to a maximum of $\sim 70\,\mathrm{km\ s^{-1}}$ in the oldest, high $\mathrm{[\alpha/Fe]}$ bins. The spread in $\sigma_{z}$ at fixed age is greater than that in $\sigma_{R}$ for most age bins, but these become similar at the oldest ages. For the younger, low $\mathrm{[\alpha/Fe]}$ populations, with $\mathrm{age}<6$ Gyr, the bins with higher $\sigma_{z}$ correspond to those with the largest $\langle R_{\mathrm{mean}}\rangle$, indicating that these stars are typically orbiting outside the solar circle. This indicates that $\sigma_{z}$ increases with $R$ at all ages. Given that this is a somewhat unexpected finding, we briefly discuss its implications in Section 5.

The increasing $\sigma_{z}$ vs. $\langle R_{\mathrm{mean}}\rangle$ trend does not appear in the high $\mathrm{[\alpha/Fe]}$ populations, although the spread in $\sigma_{z}$ is similar. For low $\mathrm{[\alpha/Fe]}$ bins, the spread in $\sigma_{R}$ at fixed age is very small, such that the bins with high $\langle R_{\mathrm{mean}}\rangle$ have an excess of $\sigma_{z}$, as discussed further below. The high $\mathrm{[\alpha/Fe]}$ populations have lower $\langle R_{\mathrm{mean}}\rangle$ than the low $\mathrm{[\alpha/Fe]}$, but have the higher $\sigma_{R}$ and $\sigma_{z}$. Importantly, the considerable spread in $\sigma_{R}$ at fixed age in the high $\mathrm{[\alpha/Fe]}$ populations (and lack thereof for low $\mathrm{[\alpha/Fe]}$ bins) is very clear.

We approximate the shape of the radial and vertical age-velocity dispersion relations (AVRs) by fitting simple power-law relationships $\sigma_{[R,z]}(\mathrm{age})\propto\mathrm{age}^{\beta_{[R,z]}}$ (similarly to those fit by, e.g. Aumer et al., 2016b; Ting & Rix, 2018) to the run of $\sigma_{z}$ and $\sigma_{R}$ as a function of age in each $\mathrm{[Fe/H]}$ bin. We only fit the relationship in $\mathrm{[Fe/H]}$ bins where there are 3 or more age bins with more than 200 stars, and use a simple maximum likelihood procedure to estimate $\beta_{R}$ and $\beta_{z}$. The resulting $\beta$ for the vertical and radial AVRs as a function of $\langle R_{\mathrm{mean}}\rangle$ is shown for the low $\mathrm{[\alpha/Fe]}$ mono-$\mathrm{[Fe/H]}$ bins in Figure 6. The left panel shows that low $\mathrm{[\alpha/Fe]}$ populations have a similarly shaped vertical AVR as a function of $\mathrm{[Fe/H]}$ (and therefore roughly constant with $R_{\mathrm{mean}}$ also), such that the mean $\beta_{z}=0.5\pm 0.1$ when these populations are combined. We also fit the high $\mathrm{[\alpha/Fe]}$ populations (but do not show these in Figure 6). For those that were fit, we find that although less well constrained, these populations have seemingly flat vertical AVRs, with $\beta_{z}=0.02^{+0.12}_{-0.03}$. We strongly emphasize, however, that age uncertainties are much larger for high $\mathrm{[\alpha/Fe]}$ stars. Therefore, contamination between these age bins is likely and this contamination would flatten any intrinsic AVR. Furthermore, it is likely that the age spread in the high $\mathrm{[\alpha/Fe]}$ stars is artificially enhanced by these uncertainties, and so the AVR for these populations has little significance.

While $\beta_{z}$ is roughly constant as a function of $\mathrm{[Fe/H]}$ (and therefore $R_{\mathrm{mean}}$) in the low and high $\mathrm{[\alpha/Fe]}$ discs, the shapes of the radial AVRs for the low $\mathrm{[\alpha/Fe]}$ populations change significantly. We find that the high $\mathrm{[\alpha/Fe]}$ population radial AVRs are consistent with one another, and again with being nearly flat, having $\beta_{R}=0.2\pm 0.2$. However, for the low $\mathrm{[\alpha/Fe]}$ disc populations we find that the radial AVRs do appear to change shape with $\mathrm{[Fe/H]}$. The values of $\beta_{R}$ against $\langle R_{\mathrm{mean}}\rangle$ (coloured by $\mathrm{[Fe/H]}$) fit for the low $\mathrm{[\alpha/Fe]}$ disc are shown in the right panel of Figure 6. It is clear from this figure that the outer disc populations (those with the lowest $\mathrm{[Fe/H]}$) have significantly flatter radial AVRs than those populations that reside closer to the solar radius. The index $\beta_{R}$ appears to increase to a maximum $\sim 0.35$ at $\langle R_{\mathrm{mean}}\rangle\sim 8$ kpc, and then decline slightly again at $\langle R_{\mathrm{mean}}\rangle$ inside the solar radius. However, within the uncertainties, all the mono-$\mathrm{[Fe/H]}$ populations with $\langle R_{\mathrm{mean}}\rangle$ $\lesssim 9$ kpc are roughly consistent with the same $\beta_{R}$. We discuss the implication of this finding in Section 5.

4.3 The shape of the velocity ellipsoid as a function of age

To examine the kinematic structure of the disc as a function of age, in Figure 7 we display the axis ratio $\sigma_{z}/\sigma_{R}$ as a function of age for the mono-age, mono-$\mathrm{[Fe/H]}$ populations. As in Figure 5, the points are coloured by $\langle R_{\mathrm{mean}}\rangle$, and high $\mathrm{[\alpha/Fe]}$ populations are plotted as open points. In the two oldest age bins ($\gtrsim 6$ Gyr), where the high $\mathrm{[\alpha/Fe]}$ stars reside, the relation between $\sigma_{z}/\sigma_{R}$ and age is roughly flat, and has very little scatter. Combining all the MCMC samples for these bins, we find $\sigma_{z}/\sigma_{R}=0.64\pm{-0.04}$. This value is in rough agreement with the measurements of $\sigma_{R}/\sigma_{z}$ in LAMOST-TGAS stars by Yu & Liu (2018).

For low $\mathrm{[\alpha/Fe]}$ populations, the relation between age and $\sigma_{z}/\sigma_{R}$ is more complex. The apparent excess of $\sigma_{z}$ seen in Figure 5 is very clear: bins at the same age with different $\langle R_{\mathrm{mean}}\rangle$ have different $\sigma_{z}/\sigma_{R}$, such that the mono-age, mono-$\mathrm{[Fe/H]}$ populations residing outside the solar radius have $\sigma_{z}/\sigma_{R}\gtrsim 0.7$, with many bins having higher $\sigma_{z}/\sigma_{R}$ than the oldest, high $\mathrm{[\alpha/Fe]}$ bins. At fixed age in the low $\mathrm{[\alpha/Fe]}$ populations, greater values of $\sigma_{z}/\sigma_{R}$ correspond to larger $\langle R_{\mathrm{mean}}\rangle$. This trend is not apparent in the high $\mathrm{[\alpha/Fe]}$ disc, where all populations have a similar $\langle R_{\mathrm{mean}}\rangle$, and correspondingly similar $\sigma_{z}/\sigma_{R}$ (as discussed above). That outer disc populations are those with the highest $\sigma_{z}/\sigma_{R}$ due to an excess of $\sigma_{z}$ may seem unsurprising at first, given that it is well known that at young ages (and low $\mathrm{[\alpha/Fe]}$), the disc is strongly flared (Bovy et al., 2016; Mackereth et al., 2017) and that the disc has a significant warp in its outer regions (e.g Poggio et al., 2018). However, we discuss the possible effects of this warping and other disc heating agents further in Section 5. The more novel result here is that there is an age (between $\sim 6$ and 8 Gyr ago) at which mono-age, mono-$\mathrm{[Fe/H]}$ populations exhibit a significant change in their velocity structure, which coincides with the age that roughly divides the high and low $\mathrm{[\alpha/Fe]}$ disc populations. This may have important implications for the history of formation and assembly of the Galaxy, which we will also discuss in Section 5.

4.4 The mean velocity – $\langle R_{\mathrm{mean}}\rangle$ relationship of populations as a function of age

By allowing the mean velocities of the model velocity distributions, $v_{R,0}$ and $v_{z,0}$, to be free parameters in the fits, we can assess the degree by which the velocity distribution is subject to non-axisymmetries in the disc. We find that the general results pertaining to $\sigma_{z}$ and $\sigma_{R}$ are robust to holding the mean velocities fixed at $0\ \mathrm{km\ s^{-1}}$, and that the majority of the mono-age, mono-$\mathrm{[Fe/H]}$ populations are well fit by models with $v_{z,0}=v_{R,0}=0$. However, we find that some mono-age, mono-$\mathrm{[Fe/H]}$ populations are better fit by models that have a significant departures from $v_{[R,z],0}=0$.

In Figure 8, we show that there are clear trends between the $v_{z,0}$ and $v_{R,0}$ of mono-age, mono-$\mathrm{[Fe/H]}$ populations and the typical orbital radius $\langle R_{\mathrm{mean}}\rangle$ of the populations: young ($\lesssim 4$ Gyr), low $\mathrm{[\alpha/Fe]}$ populations have a small but significant increasing trend of their mean vertical velocity with $\langle R_{\mathrm{mean}}\rangle$, increasing to $\sim 5\ \mathrm{km\ s^{-1}}$ for the outermost populations. These populations retain the same high $\sigma_{z}/\sigma_{R}$ values (within the uncertainties) as when fit with the fixed mean velocity model. This trend of $v_{z,0}$ with $\langle R_{\mathrm{mean}}\rangle$ in young low $\mathrm{[\alpha/Fe]}$ stars is roughly consistent with the position and magnitude of the warp seen in TGAS by Schönrich & Dehnen (2018) and in Gaia DR2 by Poggio et al. (2018) (their Figure 3). The scale of the warping found in our best fit models is also roughly consistent with that found in the 2MASS data (Momany et al., 2006).

We find a different trend of the mean radial velocity with $\langle R_{\mathrm{mean}}\rangle$ for the young, low $\mathrm{[\alpha/Fe]}$ populations, such that $v_{R,0}$ has a wave-like pattern with increasing $\langle R_{\mathrm{mean}}\rangle$, dipping down to $\sim-3\ \mathrm{km\ s^{-1}}$ at $\langle R_{\mathrm{mean}}\rangle\sim 8.5$, before rising again to $\sim 5\ \mathrm{km\ s^{-1}}$ at the greatest $\langle R_{\mathrm{mean}}\rangle$. Interestingly, the older populations ($\sim 6$ Gyr) in the low $\mathrm{[\alpha/Fe]}$ disc do not show any significant trends in their mean velocities with $\langle R_{\mathrm{mean}}\rangle$. Similarly, we find that all the high $\mathrm{[\alpha/Fe]}$ populations are consistent with $v_{z,0}=v_{R,0}=0$, although these populations are far more centrally concentrated than the older, low $\mathrm{[\alpha/Fe]}$ populations. We further discuss this feature of the velocity distribution of the disc and it’s possible bearing on its history in Section 5

5 Discussion

We have shown that when regions beyond the solar vicinity are considered, the relationship between the age and the different components of the velocity dispersion of stars in the disc is complex, and co-dependent on the stars’ mean orbital positions and their element abundances. In the following section, we contextualise our results within the existing work on the velocity structure of the disc and discuss the constraints that these results place on models for the formation and evolution of the Galactic disc.

5.1 The increasing trend of $\sigma_{z}$ with $R_{\mathrm{mean}}$ in the low $\mathrm{[\alpha/Fe]}$ disc

Given that our finding that $\sigma_{z}$ shows a clear increase with $R_{\mathrm{mean}}$ at all ages is somewhat surprising, as the expectation is that discs are most likely to show a decreasing trend, we briefly discuss here some possible explanations for this finding, and present some comparison to other observational results.

For a disc in equilibrium:

[TABLE]

where $\Sigma(R)$ is the disc surface density, and $h_{z}$ is the disc scale height. It follows that if $h_{z}$ is not a function of $R$, then $\sigma_{z}$ should trace $\Sigma(R)$ and therefore be a declining function of $R$. This relationship between $\sigma_{z}$ and $\Sigma(R)$ is commonly used to measure the surface density of external galaxy discs, making an assumption of a radially constant $h_{z}$, which is difficult to measure in face-on galaxies where $\sigma_{z}$ is more readily measurable. It is commonly seen in these studies that $\sigma_{z}$ declines with $R$ (see, e.g. Martinsson et al., 2013; Aniyan et al., 2018). The assumption of constant $h_{z}(R)$ is usually based on observations of edge-on spiral galaxies with small bulges (e.g. Yoachim & Dalcanton, 2006), likely smaller than that of the Milky Way, whose bulge and bar extends out to  3-4 kpc (e.g. Wegg et al., 2015). Lastly, when comparing these results with those from external galaxies, where stellar populations are not divided by age or abundances, it is important to note that we do make such divisions here. Minchev et al. (2015) showed that measured trends in $h_{z}$ can be quite different between mono-abundance populations and the total population, and it is possible that these effects also play a role in the $\sigma_{z}$-$R$ relationship of mono-age, mono-$\mathrm{[Fe/H]}$ populations.

Aside from these points, the Milky Way does appear from these results, however, to diverge from external galaxies in this respect. Flattening and a slight increase in $\sigma_{z}$ with $R$ outside the solar radius was also found by Sanders & Das (2018). It is known that mono-age populations in the disc flare considerably (e.g. Bovy et al., 2016; Mackereth et al., 2017). Since $h_{z}$ then depends on $R$, $\sigma_{z}$ can stay constant or increase with $R$. The existence of the extended bar in the Milky Way may also affect this relationship, and the specific details of its structure and dynamics, and interaction with disc stars, is as yet not well determined. Our finding is perhaps more logical still, given that the disc is likely not in equilibrium, but being perturbed by external forces such as interactions with satellites.

5.2 Heating agents and the velocity dispersion history throughout the disc

Our understanding of heating in the Milky Way disc has thus far largely been informed by observational constraints that were confined relatively near to the Sun. Our data set extends well beyond the solar vicinity, out to $R\gtrsim 12$ kpc in the outer disc and into the inner disc down to $R\lesssim 5$ kpc. This allows us to build a picture of heating and its connection to Galactic evolution over the extent of the disc. While a more robust understanding of these results will likely only be afforded by detailed numerical simulations, we discuss the implications of our findings using the extensive body of existing work on this topic. We attempt to understand trends in our results from the context of the more classical results on this topic, which, although outdated in some regards by numerical simulations, provide some key insights into the processes at hand. Providing reference to the most up-to-date theoretical work on the topic (using numerical simulations, e.g. those examined in Aumer et al., 2016b, a, 2017; Aumer & Binney, 2017) allows us to consolidate these findings into a coherent picture of the heating processes which have shaped the Milky Way disc.

The majority of studies that aim to understand the history of the disc from its velocity structure concentrate on two main observables, which are those that we have measured and discussed in Section 4: a) the shape of the age-velocity dispersion relations, usually characterised by a power law with index $\beta$, and b) the ratio of the vertical to radial velocity dispersion $\sigma_{z}/\sigma_{R}$. Along with these observables, two main heating agents are generally considered: scattering by small, localized perturbations in the potential such as GMCs (e.g. Spitzer & Schwarzschild, 1951, 1953; Lacey, 1984) and scattering by spiral arms (e.g. Sellwood & Carlberg, 1984; Carlberg & Sellwood, 1985). The general picture that has emerged is that GMCs are responsible for much of the heating in the vertical direction, whereas spiral arms act to heat the disc radially. Spirals contribute little to the vertical heating because the scale lengths of the spiral irregularities ought to be much longer than the vertical scale of the stellar disc (e.g. Jenkins & Binney, 1990).

5.2.1 Velocity ellipsoid shape

Analytical formulations of the (vertical) heating due to GMCs by Lacey (1984) showed that the velocity dispersion should increase as $\sigma_{z}\propto\tau^{0.25}$, and that the axis ratio $\sigma_{z}/\sigma_{R}$ should approach a value of $\sim 0.8$. Later work, using more sophisticated calculations, found that the axis ratio resulting from GMC heating is more likely somewhere between $\sim 0.5$ to $0.6$, but found a similar power-law index of $\sim 0.3$ (Ida et al., 1993; Hänninen & Flynn, 2002). We find that all high $\mathrm{[\alpha/Fe]}$ populations are consistent with having $\sigma_{z}/\sigma_{R}=0.64\pm 0.04$, in rough agreement with this revised value for GMC heating (we return to the AVR later in Section 5.2.2). In addition to this, it is interesting to note that all high $\mathrm{[\alpha/Fe]}$ populations have similar, low $\langle R_{\mathrm{mean}}\rangle$, but large variations in velocity dispersion with $\mathrm{[Fe/H]}$ in a fixed age bin. That these populations were presumably formed at similar times and have different radial and vertical velocity dispersions but a similar ratio of these may indicate that these stars in fact formed in a turbulent ISM (external galaxy observations between redshift 0.7 and 2.3 indicate that the ISM is indeed turbulent at these early times, e.g. Wisnioski et al., 2015), and have simply retained their kinematics. Understanding this will, however, require further modelling of the formation of these stellar populations in self-consistent cosmological simulations.

In the low $\mathrm{[\alpha/Fe]}$ disc, the ratio $\sigma_{z}/\sigma_{R}$ varies strongly as a function of $\mathrm{[Fe/H]}$ (and therefore $\langle R_{\mathrm{mean}}\rangle$) and age, with many of the low $\mathrm{[Fe/H]}$ (outer disc) populations having $\sigma_{z}/\sigma_{R}$ exceeding that predicted by GMC models, and the higher $\mathrm{[Fe/H]}$ (inner disc) populations falling well below the predicted value. The fact that each $\mathrm{[Fe/H]}$ bin in the low $\mathrm{[\alpha/Fe]}$ disc also shows an increasing trend between $\sigma_{z}/\sigma_{R}$ and age also indicates that the kinematics of this population cannot be explained by GMC heating alone, because models predict a flat $\sigma_{z}/\sigma_{R}$ with age if a disc is heated by a single stationary process (this is seen, for example, in the work of Jenkins & Binney, 1990). Predictions of $\sigma_{z}/\sigma_{R}$-age relations from more recent numerical simulations which include a changing mass fraction of GMCs over time show clearly similarly increasing trends to those seen in our results (Aumer et al., 2016b). In those simulations, GMC scattering is important for the early disc, and slowly gives way to spiral arm heating, as the disc grows in mass. Our result that $\sigma_{z}/\sigma_{R}$ increases with age in the low $\mathrm{[\alpha/Fe]}$ disc is clearly in-line with this picture of heating by non-stationary processes. A similar scenario was proposed also by Yu & Liu (2018), whose results on the velocity ellipsoid shape are consistent with those presented here.

Jenkins & Binney (1990) showed that increasing the importance of spirals over GMCs reduced $\sigma_{z}/\sigma_{R}$, such that models that have the most dominant spiral arm perturbations have the lowest $\sigma_{z}/\sigma_{R}$. This is also borne out in numerical simulations (e.g. Aumer et al., 2016b), where spirals heat generally only in the plane (increasing $\sigma_{R}$), whereas GMCs heat both in-plane and vertically. By this logic $\sigma_{z}/\sigma_{R}$ becomes a good measure of the relative importance of these processes in shaping the disc, but it is important to bear in mind also the external effect of satellite interactions which should also perturb and heat the disc (such interactions are not included, for example, in Aumer et al., 2016b). We find that in the low $\mathrm{[\alpha/Fe]}$ disc, the outermost disc populations (at low $\mathrm{[Fe/H]}$) have the highest $\sigma_{z}/\sigma_{R}$ at fixed age. This is consistent with the picture of declining importance of spiral arm heating with Galactocentric radius. However, an important finding here is that the highest values of $\sigma_{z}/\sigma_{R}$ in the outer disc far exceed that expected from GMC or spiral heating. The resulting $\sigma_{z}/\sigma_{R}$ for the low $\mathrm{[\alpha/Fe]}$ populations in the outer disc also exceeds that found in many of the idealised simulations of Aumer et al. (2016b). It is logical to assume that the number of GMCs in the outer disc is less than that in the inner disc (indeed, this is readily observed in the Milky Way, e.g. Miville-Deschênes et al., 2017), and likely even less at earlier times, while the disc was still growing. Therefore, it seems unlikely that GMCs or spiral arms could have heated these outer disc populations to their very high velocity dispersion. It is notable that this increasing vertical velocity dispersion with R is also realised in the results of Sanders & Das (2018) (their Figure 9). Fully understanding this excess of vertical velocity dispersion will almost certainly require additional modelling of the spatial density structure of the disc, which can then be compared with the kinematic modelling performed here.

5.2.2 The age-velocity dispersion relation

It is also notable that we find that $\sigma_{z}\propto\tau^{0.5}$ for the low $\mathrm{[\alpha/Fe]}$ disc populations, in general agreement with previous observational work (e.g. Wielen, 1977; Seabroke & Gilmore, 2007; Soubiran et al., 2008) and in disagreement the classical theoretical value for GMC heating. This discrepancy has been widely noted in the literature, and was cast by Aumer et al. (2016b) as arising from the discrepancy between AVRs and true heating histories, e.g. the evolution of velocity dispersion in co-eval populations which generates the AVR. They find that heating histories (described by a power law index $\tilde{\beta}_{R}$), which are mainly set in the few Gyr after stellar birth, tend to depend on the birth time of stars itself, such that stars born at different times in the Galaxy’s history attain different kinematics due to the changing nature of perturbers with time. This leads to a discrepancy in the AVR from simple analytical values, which assume heating to be a stationary process.

Considering radial heating, the simulations of Aumer et al. (2016b) predict that $\beta_{R}$ should vary with Galactocentric radius, and also predict that the heating histories of younger populations should have greater $\tilde{\beta_{R}}$, likely due to the fact that in their simulations the GMC heating gives over to spiral arm perturbations (this was also predicted by Jenkins & Binney, 1990, although they failed to predict correct indices, for a variety of reasons). We find that $\beta_{R}$ shows a clear trend with $\langle R_{\mathrm{mean}}\rangle$ in the low $\mathrm{[\alpha/Fe]}$ disc, such that for the outer disc stars $\beta_{R}\sim 0.15$. As $\langle R_{\mathrm{mean}}\rangle$ decreases to $\lesssim 10$ kpc, $\beta_{R}$ rapidly increases to values between $\sim 0.3$ and $0.4$. The decreased $\beta_{R}$ in the outer disc stellar populations may be a manifestation of the waning influence of spiral arm heating on the kinematics of low $\mathrm{[\alpha/Fe]}$ stellar populations that reside in the outermost part of the disc. If spiral arm heating is acting strongly inside $\sim 10$ kpc, peaking at $8<\langle R_{\mathrm{mean}}\rangle<9$ kpc, then this may also suggest that radial migration of stars is most efficient in this region of the disc. We also find, at relatively low significance, that $\beta_{R}$ increases slightly between the innermost populations and those at $\sim 8$ kpc.

5.3 Warping of the disc as a heating agent

In the absence of extensive GMC or spiral arm heating, there are few remaining possible causes for the large velocity dispersions that we observe in the outer disc. One interesting possibility is heating from warping of the Galactic disc, caused by satellite interactions (e.g. Kim et al., 2014), bending instabilities (e.g Khoperskov & Bertin, 2017), or misalignment of the angular momenta of the disc and dark halo (e.g. Jiang & Binney, 1999; Debattista & Sellwood, 1999), to name a few proposed processes. Such a warping of the Milky Way disc is well documented (e.g. Momany et al., 2006; Schönrich & Dehnen, 2018) and has been recently confirmed in the Gaia DR2 data (Poggio et al., 2018). Modelling of such bending modes has suggested that this process may act to increase $\sigma_{z}$ (as shown in the early work of Hunter & Toomre, 1969), therefore increasing $\sigma_{z}/\sigma_{R}$ (e.g. Martig et al., 2014). Khoperskov & Bertin (2017) suggested heating by these processes can cause $\sigma_{z}/\sigma_{R}$ to approach values as high as unity. Therefore, it is at least possible that these stars have been heated by warping or bending of the disc.

That the amplitude of warping in the low $\mathrm{[\alpha/Fe]}$ disc is stronger at low ages may suggest that the warp is likely stronger in the star-forming gas disc, causing more recently formed stars to retain the warp kinematics more readily. This is somewhat consistent with the finding of Mackereth et al. (2017) that young stars in the low $\mathrm{[\alpha/Fe]}$ disc are those which have a stronger flare in the outer parts of the disc. The warping present in the kinematics combined with the flattened radial AVR and excess $\sigma_{z}/\sigma_{R}$ of the outer disc populations, strongly suggests that the kinematics of these populations is affected by the warping of the Galactic disc.

5.4 Non-isothermality in young populations

We saw in Section 4.1 that the vertical dependence of the radial and vertical kinematics departs from isothermality for the youngest populations in the disc. This manifests itself in the $\sigma_{[R,z]}$ profiles as a function of $|z|$ as an increasing trend with increased height above the midplane. This feature was first hinted at at low significance in the dependence of $\sigma_{z}$ on $z$ for low $\mathrm{[\alpha/Fe]}$ populations by Bovy et al. (2012c), but is now well established from our analysis in this paper. We now investigate this further by also studying the shape of the velocity ellipsoid $\sigma_{z}/\sigma_{R}$ (as shown in Figure 7) as a function of $|z|$. These profiles for mono-age, mono-$\mathrm{[Fe/H]}$ populations in the low and high $\mathrm{[\alpha/Fe]}$ disc are shown in Figure 9.

It is immediately clear from Figure 9 that the younger, low $\mathrm{[\alpha/Fe]}$ $\sigma_{z}/\sigma_{R}$ profiles increase with $|z|$, whereas the high $\mathrm{[\alpha/Fe]}$, old population profiles are flat over a large range of $|z|$. Following from the discussion of the work of, e.g., Aumer et al. (2016b) in Section 5.2, a decreased $\sigma_{z}/\sigma_{R}$ is likely a signature of stronger heating by spiral structure. This presents a compelling explanation for the departure from isothermality of these young populations, and may also indicate the differing timescale of these heating agents. For example, if spirals can radially heat the stellar population quickly, this would reduce the $\sigma_{z}/\sigma_{R}$ in the plane for young populations. As GMCs likely redirect this radial heating into vertical heat more slowly than the spirals heat the disc, then one might expect also a trend in $\sigma_{z}$ with $|z|$ for the younger populations, as the stars which have been vertically heated are at this point more likely to be above the plane. For older ages, this trend becomes erased as the populations approach isothermality due to prolonged exposure to this combination of heating agents. It is, however, apparent in Figure 9 that some of the older low $\mathrm{[\alpha/Fe]}$ populations still show an increasing $\sigma_{z}/\sigma_{R}$ profile. As these populations are those residing in the outermost regions of the disc, we propose that this may be due again to the waning effect of the spiral structure, reducing the flattening at old ages of these profiles.

That such trends are not seen at all in the high $\mathrm{[\alpha/Fe]}$ populations further suggests a distinct evolutionary history for these populations, which we will discuss further in the following section.

5.5 Kinematic distinction of the high and low $\mathrm{[\alpha/Fe]}$ discs

Our results demonstrate that there are many kinematic differences between the high and low $\mathrm{[\alpha/Fe]}$ discs. Not only do these populations differ in their velocity dispersions $\sigma_{R}$ and $\sigma_{z}$, but they appear to have an entirely different kinematic structure, displaying different trends in their kinematics with age, $\mathrm{[Fe/H]}$, and position in the Galaxy.

It is well known from previous work that the high $\mathrm{[\alpha/Fe]}$ disc populations are old, have the lowest $\langle R_{\mathrm{mean}}\rangle$, and the highest velocity dispersions. In this paper, we have determined that the high $\mathrm{[\alpha/Fe]}$ disc populations are furthermore characterised by near-flat radial and vertical AVRs and that they are consistent with a single $\sigma_{z}/\sigma_{R}=0.64\pm 0.04$, which does not vary with age, $\mathrm{[Fe/H]}$, or $|z|$. These trends are markedly different from those displayed by the younger, low $\mathrm{[\alpha/Fe]}$ populations that show complex trends in $\sigma_{z}/\sigma_{R}$ with age and $\mathrm{[Fe/H]}$, and have increasing values of this ratio with $|z|$. High $\mathrm{[\alpha/Fe]}$ stellar populations appear to be more centrally concentrated, with smaller scale lengths and larger scale height than their low $\mathrm{[\alpha/Fe]}$ counterparts (Cheng et al., 2012; Bovy et al., 2012b; Bovy et al., 2016; Mackereth et al., 2017), and have the smallest $\langle R_{\mathrm{mean}}\rangle$ of all the mono-age mono-$\mathrm{[Fe/H]}$ populations. A centrally concentrated and kinematically distinct high $\mathrm{[\alpha/Fe]}$ component is a key prediction in models which invoke an early, rapid formation of the high $\mathrm{[\alpha/Fe]}$ population, induced by rapid gas inflow or high frequency gas-rich mergers (Chiappini et al., 1997, 2001; Brook et al., 2004; Mackereth et al., 2018b). The slight discontinuity in age vs. velocity dispersion between the high and low $\mathrm{[\alpha/Fe]}$ discs (seen in Figure 5) is also suggestive of a different, more abrupt origin of the high $\mathrm{[\alpha/Fe]}$ populations than that predicted by, e.g. Schönrich & Binney (2009a, b). It was also predicted by Aumer et al. (2016a) that thick disc components do not arise during quiescent disc growth.

The low $\mathrm{[\alpha/Fe]}$ populations are extended in $R$, and have a rich kinematic structure in age and $\mathrm{[Fe/H]}$, consistent with a gradual formation over a longer timescale, and a subsequent evolution that was likely affected by spiral structure, GMCs, and satellite interactions in ways which strongly depend on position in the Galaxy. As discussed in Section 5.2, we find that the low $\mathrm{[\alpha/Fe]}$ populations have vertical AVRs that agree well with those measured previously, where $\sigma_{z}\propto\tau^{0.5}$, and disagree with theoretical predictions that $\sigma_{z}\propto\tau^{0.25}$ for GMC and spiral heating. The Schönrich & Binney (2009a, b) model suggests that the increased slope of the AVR may be due to biases from high $\sigma_{z}$, old, high $\mathrm{[\alpha/Fe]}$ stars which have migrated to the solar radius. This bias should therefore be removed when fitting the AVR of only low $\mathrm{[\alpha/Fe]}$ stars. We do not find that the AVR slope is decreased when fitting these populations separately. As pointed out in Ting & Rix (2018), the AVR may be steepened by invoking a larger number of scatterers (GMCs) in the past. This, contrasted with the relatively flat vertical AVR of the high $\mathrm{[\alpha/Fe]}$ populations, further suggests that these populations had very different origins and subsequent evolutionary histories. The emergent picture then from the kinematics of the low and high $\mathrm{[\alpha/Fe]}$ disc components is that they are quite different, and therefore likely formed and evolved differently and certainly over very different timescales.

6 Summary and Conclusions

In this paper, we have performed a detailed dissection of the Milky Way’s low and high $\mathrm{[\alpha/Fe]}$ disc kinematics as a function of $\mathrm{[Fe/H]}$ and age. We derive a new set of ages for APOGEE DR14 red giant branch stars by applying a Bayesian Convolutional Neural Network (BCNN) model, trained on stars with APOGEE spectra and asteroseismic ages from the APOKASC catalogue. The improved age precision, alongside the exquisite abundances afforded by APOGEE allows for a new view of the Galactic disc in age-$\mathrm{[Fe/H]}$-$\mathrm{[\alpha/Fe]}$ space. The large crossover with the Gaia DR2 catalogue affords us an unparalleled catalogue of stars with measured ages, abundances, positions and kinematics with which to place constraints on models for the formation and evolution of the Galaxy. Here, we have focused on determining the kinematic structure of the disc in age-$\mathrm{[Fe/H]}$-$\mathrm{[\alpha/Fe]}$ space. Our main findings can be summarised as follows:

• •

The radial and vertical velocity dispersions of mono-age, mono-$\mathrm{[Fe/H]}$ populations in the high and low $\mathrm{[\alpha/Fe]}$ discs generally show little variation with $|z|$—i.e., most are close to isothermal—and have slowly decreasing exponential dependencies on $R$. The youngest mono-age populations do however display a departure from isothermality as they slightly increase in $\sigma_{R}$ and $\sigma_{z}$ with $|z|$ (see Figures 4 and 9).

• •

The tilt of the velocity ellipsoid is consistent across all mono-age, mono-$\mathrm{[Fe/H]}$ populations with alignment with the Galactocentric spherical coordinate system, but with very large uncertainties, i.e. the velocity ellipsoid points toward the Galactic centre. We find our fits of the tilt angle as a function of $z$ to be in good agreement with the findings of, e.g. Büdenbender et al. (2015) and Binney et al. (2014), confirming that radial and vertical motions of stars in mono-age, mono-$\mathrm{[Fe/H]}$ groups are coupled.

• •

In the age bins considered, the radial and vertical velocity dispersion increases roughly monotonically with age (Figure 5). High $\mathrm{[\alpha/Fe]}$ mono-age, mono-$\mathrm{[Fe/H]}$ populations show a spread at fixed age in both $\sigma_{R}$ and $\sigma_{z}$, whereas the low $\mathrm{[\alpha/Fe]}$ populations only have a significant spread in $\sigma_{z}$, which is correlated with the mean orbital radius of the populations.

• •

The shape of the vertical age-velocity dispersion relation does not vary between mono-age, mono-$\mathrm{[Fe/H]}$ populations in the low and high $\mathrm{[\alpha/Fe]}$ disc. We find that $\beta_{z}=0.5\pm 0.1$ for the low $\mathrm{[\alpha/Fe]}$ disc, in general agreement with previous studies (e.g Wielen, 1977; Seabroke & Gilmore, 2007; Soubiran et al., 2008) and the recent Gaia DR2 work by Ting & Rix (2018). We find $\beta_{z}=0.02^{+0.12}_{-0.03}$ (consistent with a flat AVR) for the high $\mathrm{[\alpha/Fe]}$ populations albeit at low significance, due to the large uncertainties on age in this regime.

• •

The radial age-velocity dispersion relation changes shape as a function of $\mathrm{[Fe/H]}$ for the low $\mathrm{[\alpha/Fe]}$ disc populations, with the lowest $\mathrm{[Fe/H]}$ populations having flatter radial AVRs (Figure 6). These populations have the largest $\langle R_{\mathrm{mean}}\rangle$, suggesting that the outer disc has undergone a different history of heating to the inner disc. We tentatively attribute this to the declining effect of spiral arm heating with Galactocentric radius, and the increasing effect of external heating agents (e.g. satellite mergers) on the Galactic disc.

• •

The shape of the velocity ellipsoid, defined by the axis ratio $\sigma_{z}/\sigma_{R}$, varies in a complex manner with age, $\mathrm{[Fe/H]}$ and $\mathrm{[\alpha/Fe]}$ (See Figure 7). High $\mathrm{[\alpha/Fe]}$ populations have a constant $\sigma_{z}/\sigma_{R}=0.64\pm 0.04$, which does not vary with age or $\mathrm{[Fe/H]}$. The low $\mathrm{[\alpha/Fe]}$ populations’ $\sigma_{z}/\sigma_{R}$ display clear variations with age, and $\mathrm{[Fe/H]}$, which is manifest also as a trend between the kinematics and their mean orbital radius $\langle R_{\mathrm{mean}}\rangle$. At each age, outer disc (low $\mathrm{[Fe/H]}$) populations have the highest $\sigma_{z}/\sigma_{R}$. This increased $\sigma_{z}/\sigma_{R}$ is due to the ‘excess’ $\sigma_{z}$ in the outer disc, in combination with the flattened radial AVR of these populations, reinforcing the notion that the outer disc is likely to have been heated by some process other than GMC scattering or spiral arm heating, which predict smaller $\sigma_{z}/\sigma_{R}$.

• •

Additionally, we fit models where the mean velocity in the vertical and radial direction is allowed to vary freely. In general, mono-age, mono-$\mathrm{[Fe/H]}$ populations are well fit by models with $v_{z,0}=v_{R,0}=0$. But we find that the young ($\lesssim 4$ Gyr) low $\mathrm{[\alpha/Fe]}$ populations with larger $\langle R_{\mathrm{mean}}\rangle$ seem to be better fit by models with a slightly positive $v_{[R,z],0}$, and have a wave-like trend in $v_{R,0}$ with $\langle R_{\mathrm{mean}}\rangle$ (See Figure 8). Old stars in the low $\mathrm{[\alpha/Fe]}$ disc do not appear to be part of this warping. We contend that these shifts are likely due to the Galactic warp, which may be partially responsible for the excess heating in the outer disc.

Our main findings place novel constraints on the history of dynamical heating across the disc. The finding that the velocity ellipsoid shape is strongly dependent on Galactocentric radius for low $\mathrm{[\alpha/Fe]}$ populations, and that the AVR becomes flatter at greater radii is strongly suggestive that the Milky Way disc has been heated by external perturbers. We also propose that the changing AVR shape with $\langle R_{\mathrm{mean}}\rangle$ provides insight into where heating from spiral structure becomes most important, and may indicate where radial migration is most prevalent in the disc (for example, our results suggest that radial migration should be most efficient at $8<R<9$ kpc; see Section 5.2). We also find that the high and low $\mathrm{[\alpha/Fe]}$ disc components appear to be kinematically distinct, reinforcing the notion that these structures formed differently.

These results pose many new constraints which must be considered in any model which attempts to explain the kinematics of the Milky Way disc. Importantly, they strongly suggest that external effects may play an important role in setting the kinematic structure of the outer disc. Furthermore, our results point to the need for the inclusion of non-axisymmetries in any modelling, as the excess heating seen in the data may be caused by the induced warping of the disc, and the kinematics of the low $\mathrm{[\alpha/Fe]}$ population may be affected by perturbations from spiral structure. Indeed, our results also suggest that fitting the observed disc kinematics with axisymmetric models (as we do here) may be limiting the inferences that can be made from the data.

Particularly, these results make clear the need for self-consistent modelling using large volume numerical simulations in a cosmological context, which provide statistical samples of galaxies evolving in different environments and with different assembly histories, but also where the stellar dynamics and kinematics are adequately resolved. In this regard, the data are somewhat leading the theory, as simulations of large volumes tend to poorly resolve (or not model) the cold gas phases required to form stars in disc galaxies with realistic kinematics and subsequently realistic dynamical heating processes. However, the most recent zoom-in simulations are now beginning to resolve colder gas phases, and self consistently model heating from GMCs, for example (e.g. Wetzel et al., 2016; Grand et al., 2017). An advantage of using such simulations is that non-axisymmetries are also modelled self-consistently, allowing the combined effects of bars, spiral arms, warps and waves on the galaxy kinematics to be assessed.

Acknowledgements

The authors thank the anonymous referee for a useful report. We also thank Ivan Minchev, Jerry Sellwood, Chao Liu and Michael Hayden for helpful discussions on the originally submitted manuscript. We also thank Leandro Silva and Martin Smith for highlighting necessary typographic corrections. JTM is funded by an STFC studentship, and is grateful for funding from the Royal Astronomical Society and the Dunlap Institute for Astronomy and Astrophysics at the University of Toronto and thanks them for their hospitality during an extended visit while preparing this work. JTM and AM acknowledge support from the ERC Consolidator Grant funding scheme (project ASTEROCHRONOMETRY, G.A. n. 772293). JB acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), funding reference number RGPIN-2015-05235, and from an Alfred P. Sloan Fellowship. JT acknowledges that support for this work was provided by NASA through the NASA Hubble Fellowship grant #51424 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. This project was developed in part at the 2017 Heidelberg Gaia Sprint, hosted by the Max-Planck-Institut für Astronomie, Heidelberg, and the 2018 NYC Gaia Sprint, hosted by the Center for Computational Astrophysics of the Flatiron Institute in New York City.

This research made use of the cross-match service provided by CDS, Strasbourg. Analyses and plots presented in this article used iPython, and packages in the SciPy ecosystem (Jones et al., 2001; Hunter, 2007; Perez & Granger, 2007; van der Walt et al., 2011). The study made use of high performance computing facilities at Liverpool John Moores University, partly funded by the Royal Society and LJMU’s Faculty of Engineering and Technology.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://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.

Appendix A Estimating ages for Red giant stars using Bayesian COnvolutional Neural Networks

The APOKASC-2 catalogue (Pinsonneault et al., 2018) consists of asteroseismic and spectroscopic data products measured for 6,676 stars in common between the APOGEE and NASA-Kepler mission. This represents a significantly larger catalogue than the previous release, as well as a number of improvements over the original analysis (presented and examined in Pinsonneault et al., 2014). APOKASC DR2 derives stellar properties by combining spectroscopically measured parameters from APOGEE DR14 with light curves from Kepler, the improved reduction of which is described by Elsworth et al. (2018, in prep.). The asteroseismic analysis is performed using a set of five distinct pipelines, where the mean results of these are used to generate the final catalogue. Pipelines estimate the asteroseismic parameters, the large frequency separation $\Delta\nu$ and the frequency of maximum power, $\nu_{\mathrm{max}}$, which are then used in semi-empirical scaling relations that give the mass and radius of the star, given a known $T_{\mathrm{eff}}$ (the spectroscopic data). Scatter about the ensemble mean of the pipelines is used to infer the random uncertainties on the asteroseismic parameters $\Delta\nu$ and $\nu_{\mathrm{max}}$, which are propagated through to the stellar properties. Any differences between the pipelines that depend on $\Delta\nu$ and $\nu_{\mathrm{max}}$ contribute to a systematic error budget. A key difference between APOKASC-1 and 2 lies in the treatment of the $\Delta\nu$ scaling relation, which is no longer considered to be an exact relation, and is corrected on a star-by-star basis, through the approach described by Serenlli et al. (2018, in prep.). This theoretical correction also relies on the APOGEE DR14 $\mathrm{[Fe/H]}$ and $\mathrm{[\alpha/Fe]}$ values. The property of most importance to this work, the stellar ages, are finally estimated using the derived mass, surface gravity, and abundances with their uncertainties estimated through propagation of the uncertainties on those properties. We use the ‘recommended’ age estimates and uncertainties as tabulated and available in Pinsonneault et al. (2018), without applying any external corrections.

The APOKASC-2 ages and their corresponding combined APOGEE spectra form the basis of the training set which we use here to train a model to predict ages for the remainder of the DR14 red-giant sample. This method relies mainly on the relation between the surface abundances of red giants and their main sequence mass, which is set when the stars undergo first dredge up (FDU). Stars on the main sequence undergo hydrogen fusion by the CNO and CN cycles, in which C, N and O atoms catalyse burning, and so undergo evolution in their own relative abundances. As shown by Salaris et al. (2015), the abundance profile of Nitrogen going inwards to the stellar core (before FDU) reflects first that of the CN cycle, followed by the CNO cycle. The FDU mixes the envelope with the CN layer. Stellar models with increased mass tend to have larger convective envelopes, meaning that they reach deeper into the stellar interior, meaning that the surface $\mathrm{[C/N]}$ decreases with increasing RGB mass. We assume here that such differences in surface abundances of C, N and O are reflected in APOGEE spectra, such that there should be some (likely complicated) relation between the spectra and the stellar mass and, therefore, the implied age. In Martig et al. (2016), this relationship was exploited by fitting a polynomial feature regression model to predict age from APOGEE C, N, $T_{\mathrm{eff}}$, $\log{g}$ and $\mathrm{[Fe/H]}$. Here, we use the data in the full spectrum, applying a neural network based model, in an attempt to improve age estimates by capturing the more complicated aspects of the relationship between the age and the surface abundances, and avoiding the reliance of the age estimation on the individual abundance determinations from the spectra.

Our method uses a Bayesian Convolutional Neural Network (BCNN), implemented in the astroNN python package (Leung & Bovy, 2019b), which wraps the Keras and TensorFlow machine learning architectures. BCNNs treat the more commonly used Convolutional Neural Network (CNN) as a Bayesian regression problem, inferring the probability distributions over the model weights. This not only allows for the propagation of uncertainty into any predictions of the model, but also offers increased robustness to overfitting, which is common in CNNs with small training data-sets (as is the case here). The use of BCNNs for the prediction of abundances for APOGEE spectra is explored in Leung & Bovy (2019b). Here we focus only on their use for predicting stellar ages.

The training data is compiled and loaded using functions in astroNN. Individual APOGEE visit spectra are recombined and continuum normalised using a method similar to that used in the Cannon (Ness et al., 2015), which makes a Chebyshev polynomial fit to specifically selected pixels separately between the different CCD chips in APOGEE. This procedure results in an improved normalisation, which is preferable when training the BCNN model, over the normalisation used in the standard APOGEE data reduction. The model consists of 2 fully connected hidden layers, with 128 and 64 nodes (in that order), activated by the ReLU function (Nair & Hinton, 2010). Before this, two convolution layers process the spectra through 2, then 4 filters. Training data is fed to the model in batches of 32 spectra, and the model is trained for 60 epochs (where one epoch represents 166 batches in total, for a training set size of 5311 spectra, with a set of 590 validation spectra). The model is trained using a mean-squared-error loss function, with an initial learning rate of 0.005, which is reduced when the loss reaches a plateau. 500 spectra are retained from the training altogether, and used as a test set.

Figure 10 shows the BCNN predicted age, $\mathrm{age}_{\mathrm{out}}$, against the input age from APOKASC-2, $\mathrm{age}_{\mathrm{in}}$, for the training (black points) and test samples (shown in blue, with their associated uncertainties, as estimated through the Bayesian dropout variational inference). This data is included as supplementary material. The bottom panel gives the residual $\Delta\mathrm{age}$ defined as $\mathrm{age}_{\mathrm{out}}-\mathrm{age}_{\mathrm{in}}$. The standard deviation of residuals across the whole range is 36%, corresponding to a median $\Delta\mathrm{age}=1.26$ Gyr. For input ages less than 5 Gyr, the scatter is reduced to 30%. It is clear, however that there is a signifcant trend such that young ages are over-estimated, and old ages are under-estimated, with a stronger effect at old age. Stars with input ages of $\sim 10$ Gyr are underestimated in the test set by as much as 3.5 Gyr, with a large scatter, such that no stars are correctly predicted or over-predicted. This behaviour has been seen in other methods that exploit asteroseismology, using simpler polynomial regression method based on abundance measurements to measure ages (e.g. Martig et al., 2016). It is noteworthy that it persists even when the full scope of information in the spectra are used. We also tested the effect on the predicted ages of adjustments to the training data, implying progressively larger cuts in the asteroseismic $\nu_{\mathrm{max}}$ in order to remove stars which lie high on the RGB, where the asteroseismic models are less well tested. We find that implying these cuts makes little difference to the final results of the training, and the output ages from the test set match well with those predicted from the full training set in all cases. Dividing the sample between RC and RGB stars (as determined by the asteroseismology) and training the BCNN separately on these populations also provides comparable results, implying that the neural network is able to capture this information from the spectra and account for it when training on the entire APOKASC-2 sample (the presence of such information in stellar spectra has been discussed by Hawkins et al., 2018; Ting et al., 2018).

We further examine the issue with the absolute scaling of ages predicted from the neural network (and other) models by comparing the predicted age-$\mathrm{[C/N]}$ relationship with that obtained from stellar evolution models. In Figure 11, we show the BCNN predicted ages (upper panel) and the training set APOKASC-2 ages (lower panel) against the APOGEE DR14 $\mathrm{[C/N]}$ for stars in APOKASC-2 with $-0.05<\mathrm{[Fe/H]}<0.05$ dex. The age-$\mathrm{[C/N]}$ relationship for solar metallicity stars, as obtained from the BASTI (Salaris et al., 2015) and STAREVOL (Lagarde et al., 2017) models, is overplotted in the dotted and dashed lines. It is apparent that the models roughly agree with the training and predicted data at young ages ($\lesssim 5$ Gyr). The Salaris et al. (2015) model appears to agree slightly better than that of Lagarde et al. (2017), which has a systematic offset to lower $\mathrm{[C/N]}$ at any given age. At older ages, the spectroscopic ages have higher $\mathrm{[C/N]}$ than is predicted by the models, whereas the asteroseismic input ages seem to agree well over the full range of abundance ratios. This test provides a further indication that the oldest ages are under-predicted by the neural netwok model, as the $\mathrm{[C/N]}$ ratios of these stars (as measured by APOGEE) should correspond an older age, given the stellar evolution models. This implies one of two possible explanations, either: a) the ages of old stars provided in the training set are underestimated, or b) that the spectroscopic information is not representative of the true mass (as measured by asteroseismology) at old age. As the training set appear to match the models relatively well, it would seem that the second explanation is more likely. Fully understanding this effect is beyond the scope of this paper, but will be essential for accurately determining the ages of the oldest stars in the Milky Way using spectroscopy.

In a previous paper in this series, Mackereth et al. (2017), we chose to correct for the under-prediction of old ages using a non-parametric LOWESS fit to the $\mathrm{age}_{\mathrm{out}}-\mathrm{age}_{\mathrm{in}}$ relation. Here, we leave under-predicted ages in the full data set, and simply warn readers that results pertaining to any trends in age at values greater than $\sim 8$ Gyr likely extend from that age to the oldest stars in the sample. The above tests demonstrate that the sample likely contains stars older than the $\sim 10$ Gyr limit seen in the final predicted ages, above which it appears the spectroscopic information is no longer useful to predict ages. The final age catalogue for DR14 stars, as predicted by the BCNN, is included in a supplementary data table.

Appendix B Example posterior PDF of velocity distribution model

To demonstrate the overall quality of the fit to the velocity distribution that can be achieved using the APOGEE-Gaia sample, we provide here an example of the posterior probability distribution gained from the MCMC sampling of the parameters for the velocity distribution of a single exemplary bin in age-$\mathrm{[Fe/H]}$. The bin which we choose to show here is at intermediate age ($4.0<\mathrm{age}<4.5\ \mathrm{Gyr}$ old) and sub-solar $\mathrm{[Fe/H]}$ ($-0.1<\mathrm{[Fe/H]}<0.$) in the low $\mathrm{[\alpha/Fe]}$ population, which has 2178 stars. The posterior distribution of each of the parameters and their covariance is shown in a corner plot, in Figure 12. The posteriors are in most cases well behaved and approximately Gaussian. As expressed in Section 3, decreasing the number of stars used for the fit broadens the distributions, but the median values remain the same within the uncertainties. Checking the posterior distributions for all the mono-age, mono-$\mathrm{[Fe/H]}$ bins