The vertical metallicity gradients of mono-age stellar populations in the Milky Way with the RAVE and Gaia data

TL;DR

This study analyzes how the vertical metallicity gradients of mono-age stellar populations in the Milky Way vary with age, revealing a trend of increasing gradient steepness and metallicity dispersion with age, supporting a flaring disc formation scenario.

Contribution

It provides the first detailed analysis of vertical metallicity gradients across multiple mono-age populations using combined RAVE and Gaia data.

Findings

No vertical metallicity gradient in the youngest stars.

Steeper negative gradients in older populations.

Metallicity dispersion increases with age.

Abstract

We investigate the vertical metallicity gradients of five mono-age stellar populations between 0 and 11 Gyr for a sample of 18 435 dwarf stars selected from the cross-matched Tycho-Gaia Astrometric Solution (TGAS) and RAdial Velocity Experiment (RAVE) Data Release 5. We find a correlation between the vertical metallicity gradients and age, with no vertical metallicity gradient in the youngest population and an increasingly steeper negative vertical metallicity gradient for the older stellar populations. The metallicity at disc plane remains almost constant between 2 and 8 Gyr, and it becomes significantly lower for the $8 < \tau \leqslant 11$ Gyr population. The current analysis also reveals that the intrinsic dispersion in metallicity increases steadily with age. We discuss that our results are consistent with a scenario that (thin) disc stars formed from a flaring (thin) star-forming disc.

Tables

Table 1. Table 1: Results from the linear regression model.

	unweighted				weighted
Age	d[M/H]/d$\left|Z\right|$	$\beta$	$\sigma$		d[M/H]/d$\left|Z\right|$	$\beta$	$\sigma$
Gyr	dex kpc-1	dex	dex		dex kpc-1	dex	dex
0 - 2	$-{0.007}_{0.070}^{0.070}$	${0.016}_{0.015}^{0.015}$	${0.131}_{0.006}^{0.006}$		${0.014}_{0.068}^{0.070}$	${0.003}_{0.017}^{0.016}$	${0.144}_{0.006}^{0.006}$
2 - 4	$-{0.076}_{0.023}^{0.024}$	$-{0.078}_{0.005}^{0.005}$	${0.168}_{0.002}^{0.002}$		$-{0.095}_{0.023}^{0.021}$	$-{0.075}_{0.005}^{0.005}$	${0.174}_{0.002}^{0.002}$
4 - 6	$-{0.169}_{0.037}^{0.037}$	$-{0.125}_{0.006}^{0.007}$	${0.227}_{0.002}^{0.003}$		$-{0.192}_{0.035}^{0.036}$	$-{0.111}_{0.007}^{0.007}$	${0.232}_{0.003}^{0.003}$
6 - 8	$-{0.302}_{0.075}^{0.076}$	$-{0.150}_{0.012}^{0.012}$	${0.281}_{0.004}^{0.004}$		$-{0.411}_{0.070}^{0.068}$	$-{0.113}_{0.012}^{0.012}$	${0.287}_{0.004}^{0.004}$
8 - 11	$-{0.429}_{0.127}^{0.120}$	$-{0.298}_{0.018}^{0.017}$	${0.311}_{0.006}^{0.006}$		$-{0.653}_{0.113}^{0.117}$	$-{0.218}_{0.018}^{0.017}$	${0.321}_{0.007}^{0.007}$
All stars	${0.060}_{0.022}^{0.021}$	$-{0.162}_{0.004}^{0.004}$	${0.238}_{0.001}^{0.002}$		$-{0.026}_{0.021}^{0.020}$	$-{0.140}_{0.004}^{0.004}$	${0.239}_{0.002}^{0.002}$

Equations

$$p_1(\tau ,M,\text{mbox}{[Fe/H]}_{\mathrm{init}}|T_{\mathrm{eff}},\text{mbox}[Fe/H],d)\propto p_{0}L,$$

$$L=\prod _i\frac{1}{\sqrt{2\pi }{\sigma }_{\mathrm{i}}}\times exp\left(-\frac{(O_{\mathrm{i}}-S_{\mathrm{i}}{)}^2}{2{\sigma }_{\mathrm{i}}^2}\right),$$

$${\xi }_{\mathrm{i}}$$

$${\eta }_{\mathrm{i}}|{\xi }_{\mathrm{i}}$$

$$y_{\mathrm{i}},x_{\mathrm{i}}|{\eta }_{\mathrm{i}},{\xi }_{\mathrm{i}}$$

$$x_{\mathrm{i}}$$

$$y_{\mathrm{i}}$$

Full text

††pagerange: The vertical metallicity gradients of mono-age stellar populations in the Milky Way with the RAVE and Gaia data–11††pubyear: 2017††pagerange: The vertical metallicity gradients of mono-age stellar populations in the Milky Way with the RAVE and Gaia data–11

1 Introduction

Galactic archaeology is entering a golden era with the $Gaia$ space mission, which is providing positional and proper motion information for over a billion stars (Gaia Collaboration et al., 2016a). $Gaia$’s Radial Velocity Spectrometer (Cropper & Katz, 2011) will provide radial velocity measurements for the 150 million brightest stars in the Milky Way. Large-scale stellar spectroscopic surveys, such as the Radial Velocity Experiment (RAVE; Steinmetz et al., 2006), the Sloan Extension for Galactic Understanding and Exploration (SEGUE; Yanny et al., 2009), the Large Sky Area Multi-Object Fibre Spectroscopic Telescope (LAMOST; Liu et al., 2014), the Apache Point Observatory Galactic Evolution Experiment (SDSS-III/APOGEE; Eisenstein et al., 2011; Majewski et al., 2015), Gaia-ESO (Gilmore et al., 2012), the GALactic Archaeology with HERMES (GALAH; Martell et al., 2017), the 4-metre Multi-Object Spectroscopic Telescope (4MOST; de Jong et al., 2012) and the WHT Enhanced Area Velocity Explorer (WEAVE; Dalton et al., 2012), aim to complement the $Gaia$ data with chemical and radial velocity information for the faint stars, thus allowing us to probe the Milky Way to unprecedented scale and detail. The ultimate goal of Galactic archaeology in the Milky Way is to reconstruct the Galaxy’s formation and evolution history. For this purpose, one of the most cogent approaches that employ $Gaia$ (and spectroscopic survey data) is to map the chemo-dynamical structure of the Galaxy in chronological order. In this way, we can discern the role of various internal processes in shaping the geometry, dynamics and chemistry of the Galaxy at every evolutionary stage by comparing observations with theoretical models, often in which the time evolution is the clearest indicator.

Stellar age is an essential ingredient in mapping the chrono-chemo-dynamical structure of the Milky Way. Determining the age of stars is a challenging problem. There are several ways to obtain the ages of stars, including isochrone fitting (e.g. Holmberg et al., 2009; Casagrande et al., 2011; Haywood et al., 2013), asteroseismology (e.g. Chaplin & Miglio, 2013, for a review) and the C/N abundance ratio (Masseron & Gilmore, 2015; Martig et al., 2016). For example, by using the age information in conjunction with metallicity and distance measurements, we can characterise the radial and vertical metallicity profile of the Galactic disc at every evolutionary stage. Then, we can make informed inferences about the kinematic and dynamical properties of the Galactic disc at the formation and subsequent evolutionary stages.

Using 418 red-giant stars from the CoRoT and APOGEE combined data close to the Galactic plane (6 < $R\lesssim$ 13 kpc, $\left|Z\right|$ < 0.3 kpc, here $R$ and $\left|Z\right|$ are the Galactocentric radius and vertical height from the Galactic plane, respectively), Anders et al. (2017) investigated the radial metallicity distribution of stars in different mono-age stellar populations, defined as a group of stars with similar ages. They found that, at each radius, the dispersion of the metallicity distribution increases with age. Also, their study revealed that the radial metallicity gradient flattens with age.

These observational trends can be explained as follows. In this paper, we assume that there was no metallicity dispersion for stars formed at a fixed radius as is seen for young B stars in the Solar neighbourhood (Nieva & Przybilla, 2012; Feltzing & Chiba, 2013), but there was a negative radial metallicity gradient. Then, at a fixed radius, radial mixing of stars could have brought more metal-rich stars from the inner region and more metal-poor stars from the outer region (e.g. Sellwood & Binney, 2002; Schönrich & Binney, 2009a, b; Minchev et al., 2013; Grand et al., 2015), and this would increase the metallicity dispersion. Since older stars had more time to be affected by radial mixing, the older stellar population at a fixed radius contains both more metal-poor and rich stars than the younger population, and this explains the broader metallicity dispersion of the older stars. The flattening of the radial metallicity gradient can be explained by the radial metallicity gradient at formation time being flatter at earlier epochs or radial mixing flattening the radial metallicity gradient for the older stars (e.g. Roškar et al., 2008; Sánchez-Blázquez et al., 2009; Minchev et al., 2014; Grand et al., 2015; Kawata et al., 2017).

Another study that measured the metallicity distribution in the Galactic disc as a function of age was done by Xiang et al. (2015). They employed the LAMOST Spectroscopic Survey of the Galactic Anti-Center dataset (LSS-GAC; Liu et al., 2014), and used photometric distances to perform isochrone fitting of turn-off stars to obtain stellar ages for 297 042 stars covering a vast region between 7.5 < $R$ < 13.5 kpc and $\left|Z\right|$ < 2.5 kpc. They found that the measured radial metallicity gradients are approximately constant up to the age of 11 Gyr. They considered an age-defined thick disc as the population of stars with age > 11 Gyr and found that this thick disc population has a flat radial metallicity gradient, distinct from the younger thin disc populations.

Xiang et al. (2015) also investigated the vertical metallicity gradients of mono-age stellar populations by splitting their sample into age bins of width 2 to 3 Gyr between 2 and 11 Gyr and a wider single age bin for the 11-16 Gyr population. Their results revealed that at the Solar radius, the older mono-age stellar populations display increasingly steep negative vertical metallicity gradients from d[M/H]/d$\left|Z\right|\simeq-\ 0.13$ dex kpc*-1* for stars with ages between 2 and 4 Gyr to d[M/H]/d$\left|Z\right|\simeq-\ 0.21$ dex kpc*-1* for the 8 - 11 Gyr population. However, the oldest population with age $>11$ Gyr shows a flatter vertical metallicity gradient of d[M/H]/d$\left|Z\right|=-\ 0.11$ dex kpc*-1*. Here, we use [M/H] to indicate the metallicity defined as [M/H] = $\rm log_{10}(M/H)-log_{10}(M/H)_{\sun}$, with M indicating the metal mass fraction and H the Hydrogen mass fraction in the star.

An indirect approach to determine the vertical metallicity of mono-age stellar populations relies on using [$\alpha$/Fe] as a proxy for age (e.g. Smecker-Hane & Wyse, 1992; Haywood et al., 2013; Hayden et al., 2015). Older stars are thought to be more $\alpha$-enhanced than their younger counterparts. For example, using dwarf stars selected from the SEGUE survey, Schlesinger et al. (2014) found a negligible vertical metallicity gradient for different [$\alpha$/Fe] populations above the Milky Way plane (0.3 $<\left|Z\right|<$ 0.6 kpc), where $\alpha$ and Fe indicate $\alpha$-element and iron abundances, respectively. Assuming [$\alpha$/Fe] is a proxy for age, Schlesinger et al. (2014)’s results imply that mono-age stellar populations should have similarly negligible vertical metallicity gradients, in contradiction with the Xiang et al. (2015) result. However, a recent numerical simulation by Minchev et al. (2017) argues that the [$\alpha$/Fe] may not be an irrevocable proxy for the age. What is more, Mackereth et al. (2017) showed that mono-age stellar populations flare even for the highest [$\alpha$/Fe] (see also Martig et al., 2014). Thus, robust age estimates are needed to map the chemical profile of the mono-age stellar populations of the Galactic disc.

The vertical metallicity gradients as a function of stellar age provide strong constraints on the importance of radial mixing and vertical structure of the star-forming disc at the different epochs. It is a common assumption that the (thin) disc stars formed from a thin high-density molecular gas, where vertical chemical diffusion must be efficient, and thus, the mono-age stellar populations that formed from this gas disc should have no initial vertical metallicity gradients. At the same time, the mono-age stellar populations must have had a negative radial metallicity gradient initially, as was summarised in Anders et al. (2017), as discussed above. Then, the existence of a vertical metallicity gradient in a mono-age stellar population of the current Galactic disc stars indicates that a dynamical mechanism, such as radial mixing, built up the vertical metallicity gradient (e.g. Minchev et al., 2013; Kawata et al., 2017; Schönrich & McMillan, 2017). Measuring the change of the vertical metallicity gradient with age gives strong constraints on the mechanism.

Using an N-body simulation of a simplified Milky Way-like disc model and the working assumptions of no initial vertical metallicity gradient and a negative radial metallicity gradient, Kawata et al. (2017) showed that if the mono-age disc population formed with a constant vertical scale height, then radial mixing could drive a positive vertical metallicity gradient of the mono-age stellar population. This process can happen because the metal-rich (poor) stars from the inner (outer) disc tend to end up at a higher (lower) vertical height due to their original higher (lower) vertical energy in the inner (outer) region of the disc. Then, if radial mixing is slow enough, older populations should show a more positive vertical metallicity gradient.

Kawata et al. (2017) also suggested that a mono-age stellar population can build up a negative vertical metallicity gradient if the scale height of the mono-age disc increases with the radius, i.e. flaring. This is because when the metal-poor stars in the outer region, where stars have higher vertical energy due to their greater scale height, migrate inwards, they oscillate vertically more and become dominant at the higher vertical height in the inner region, which leads to a negative vertical metallicity gradient. In this flaring star-forming disc case, if the radial mixing process is slow enough, the older populations should have a more negative vertical metallicity gradient, because more stars from the outer region migrate into the inner region and occupy the higher vertical height region.

Hence, the vertical metallicity gradient for the mono-age stellar population provides critical information for understanding the dynamical evolution in the Galactic disc. Xiang et al. (2015) did pioneering work in studying the vertical metallicity gradient of the mono-age stellar population. Further research with independent data is needed to follow up on their work. Recently, Gaia Collaboration et al. (2016a, b) released their first $Gaia$ data, including the Tycho-Gaia Astrometric Solution dataset (TGAS; Michalik et al., 2015), which provides accurate astrometry for more than 2 million bright stars. At the same time, the RAVE survey team made their fifth data release (DR5 Kunder et al., 2017), consisting of carefully calibrated metallicity and stellar parameters of more than 200 000 stars that are also present in the TGAS which provides accurate parallax, i.e. distance, information. The combined dataset of TGAS and RAVE is a unique dataset that allows us to determine the vertical metallicity gradient. Hence, this paper presents an independent, new measurement of the vertical metallicity gradients as a function of age for Solar neighbourhood stars. Our sample goes up only a few 100 pc from the plane, for which accurate parallax measurements are available. The typical errors in the metallicity measurements with RAVE are about 0.1 dex. This makes it challenging to reliably determine the vertical metallicity gradient of the level of about d[M/H]/d$\left|Z\right|$ $\sim 0.1$ dex $\rm kpc^{-1}$. Therefore, we employ the Bayesian hierarchical model of Kelly (2007) to perform robust linear regression, which is a very flexible method that enables us to measure not only the vertical metallicity gradient and intercept, but also the intrinsic dispersion of the metallicity of mono-age stellar populations. This more sophisticated approach applied to new data allows us to determine for the first time the vertical metallicity gradient of the youngest population of 0-2 Gyr.

The current paper is structured as follows. In Section 2 we describe the observational data, the method employed to obtain the measurements of vertical metallicity gradients for the mono-age stellar population and how the observational bias due to stellar population is accounted for in the analysis. Then, we define mono-age stellar populations as samples of stars belonging to different age bins. We finally describe how we fit the vertical metallicity gradients using the Bayesian hierarchical model of Kelly (2007). Section 3 presents the results of our measurement of the vertical metallicity gradient for the stars with different ages. Section 4 provides the discussion and summary of our study.

2 Methods

We used the cross-matched TGAS and RAVE DR5 catalogue to use RAVE metal abundance and stellar parameters with the TGAS astrometric data. TGAS is the primary astrometric data set in the first $Gaia$ data release (Gaia Collaboration et al., 2016a) and provides position, parallax and proper motion information for 2 057 050 stars which are in the Tycho-2 catalogue (Høg et al., 2000). The parallaxes have median uncertainties of 0.3 mas, with an additional systematic error of around 0.3 mas (Gaia Collaboration et al., 2016a; Lindegren et al., 2016). RAVE is a magnitude-limited (9 < I < 12) spectroscopic survey of stars in the Milky Way using the 1.2 m UK Schmidt Telescope of the Australian Astronomical Observatory (AAO). RAVE DR5 provides effective temperature, surface gravity, and overall metallicity carefully calibrated against other spectroscopic data and the Kepler-2 (K2) asteroseismic data (Kunder et al., 2017). From the cross-matched TGAS-RAVE DR5, we obtain parallax and metallicity information for 215 590 stars.

For this study, we select stars with Galactic latitude $\left|b\right|>25\degree$ to ensure an unbiased sky coverage by RAVE (Wojno et al., 2016). As recommended in Kunder et al. (2017), we select “normal” stars with RAVE flags that satisfy c1=c2=c3=c4=…=c20=“n", and with ALGO_CONV = 0, frac_c > 0.7, CHISQ < 2000 and SNR > 20 to obtain a sample with high-quality metallicity measurements. These high-quality cuts for metallicity and stellar parameters are complemented by a further cut in relative parallax error less than 20% for accurate distance measurements in TGAS which allows us to treat the distance errors as Gaussian. As shown below, we will minimise the selection bias induced by the selection by limiting our colour-magnitude range and taking into account stellar population bias.

Tycho-2 is 90% complete up to magnitude $V$ = 11.5 mag (Høg et al., 2000). Although the TGAS sample is not the same as the Tycho-2 sample (Arenou et al., 2017), we assume that TGAS provides a good enough level of completeness for this study in this magnitude limit. Because the TGAS selection against the Tycho-2 should not introduce any bias for the metallicity of stars at different heights which are the focus of this paper, we assume that this is a safe assumption.

The magnitudes used in this analysis are extinction-corrected using the same method described in Section 10 of Kunder et al. (2017). To control our selection bias, we select samples in $J-K$, $K$ colour-magnitude diagram (CMD). To obtain accurate age measurements, we focus on dwarf stars in this study. We found that dwarf stars roughly follow $V=K+3.28(J-K)+0.35$ mag. This scaling relation provides that $V$ < 11.5 mag corresponds to $K<11.15-3.28(J-K)$, and this magnitude limit is much brighter than $I$ < 12 mag limit where the RAVE sample is chemically unbiased (Wojno et al., 2016). In addition, to sample dwarf stars, we select stars with $0.2<J-K<0.4$ mag. As discussed below, this colour cut provides a compromised solution to minimise selection bias and maximise the sample size.

The stellar age information has been determined using a Bayesian isochrone fitting method called $\tt Elli$ (Lin et al., 2017). The uncertainties of the age measured with $\tt Elli$ are smaller than those reported in the literature due to the extra constraint from TGAS parallaxes (Lin et al., 2017).To assure reliable age estimates, we use only stars whose relative age error is less than 50%.

The stellar isochrones used in this version of $Elli$ are taken from the Dartmouth Stellar evolution database (Dotter et al., 2008), with $\rm-2.48\leq[Fe/H]_{init}\leq+0.56$ dex and $\rm-0.2\leq[\alpha/Fe]\leq+0.8$ dex. The age and mass grids range from 250 Myr to 15 Gyr, and 0.1 to 4.0 $M_{\rm\odot}$, respectively, with solar abundances adopted from Grevesse & Sauval (1998). The equation of state for stellar evolution tracks is assumed to follow the general ideal gas law for stars heavier than 0.8 $M_{\rm\odot}$, with the Debye-Hückel correction, and the FreeEOS equation of state for stars less massive than 0.8 $M_{\rm\odot}$ 111http://freeeos.sourceforge.net/. The conditions at the surface of the star are computed using the Phoenix model atmospheres with $T_{\rm eff}$ ranging from 2000 to 10 000 K and log g from 0.5 to 5.5 dex. For stars hotter than 10 000 K, Castelli & Kurucz (2003) model atmospheres are employed.

The method consists of sampling from a posterior distribution $p_{\rm 1}$ which is the product of a likelihood $L$ and the prior distribution $p_{\rm 0}$ to construct a probability distribution for the parameters of interest, which in this case, are the age ($\tau$), mass ($M$), and bulk metallicity $\rm[Fe/H]_{\rm init}$. The mathematical description for the posterior probability distribution is given as:

[TABLE]

where the likelihood $L$ is defined as:

[TABLE]

where $O_{\rm i}$ is the observed values, $T_{\rm eff}$, [Fe/H] and $M_{\rm K}$, and $S_{\rm i}$ are model values. Note that the absolute magnitude, $M_{\rm K}$, depends on the distance, $d$.

Similar to our main analysis of determining the vertical metallicity gradient described below, the posterior is sampled using Markov-Chain Monte Carlo (MCMC). The observed values, $\rm O_{\rm i}$, are compared with model values, $\rm S_{\rm i}$, obtained by interpolating $\tau,M,{\rm[Fe/H]_{\rm init}}$ over the corresponding grid values of isochrones. We did not use logg due to its large uncertainty in the RAVE data. $M_{\rm K}$ and the distance were obtained from the 2MASS apparent magnitude and the TGAS parallax, respectively. For each isochrone, the likelihood is then computed following equation (2), assuming the observed parameters have Gaussian errors and the age estimate is obtained by maximising the logarithm of the likelihood function for a given metallicity.

To validate the version of $Elli$ used in this paper, we used the $Elli$ age measurements for the sample of Casagrande et al. (2011), and Fig. 1 provides a comparison between the ages derived with and the ages shown in Casagrande et al. (2011). The difference in the ages between two methods with respect to their errors are shown in the right panel and are indicative of a reasonable agreement. Although the $Elli$ age is systematically larger than Casagrande et al. (2011), the differences are within error, which is not surprising given that different isochrones are applied in these two methods. In Appendix A, we also provide a more extensive comparison with another isochrone fitting method (McMillan et al., 2017). The results in Appendix A also demonstrate that the age measurements suffer from systematic uncertainty. However, it proves that the qualitative trends discussed in this paper are robust.

Fig. 2 shows the $J-K$ vs. $K$ colour-magnitude diagram for our selected sample (red) before applying a cut in relative age error and the matched TGAS-RAVE DR5 (black) in the selected sky coverage of $\left|b\right|>25\degree$. The left and right panels of Fig. 3 display the cumulative distribution functions in the $J-K$ and $K$ distribution, respectively, for our selected sample before applying a cut in relative age error and the TGAS and RAVE cross-matched sample (referred as “TGAS+RAVE sample") within the chosen colour ($0.2<J-K<0.4$ mag) and magnitude range. Fig. 3 shows that for $0.2<J-K<0.4$ mag, the cumulative distribution function of the $J-K$ distribution for our sample follows that of the RAVE-TGAS sample where no selection is applied. The cumulative distribution function of the $K$ distribution for our sample shows no significant deviation that of the RAVE+TGAS sample. Although the $0.2<J-K<0.4$ mag seems to introduce bias, this is our compromised solution after we explored various other colour cut ranges. Our final data set after applying a cut in relative age error includes 18 435 sources populating a limited volume in both the vertical and radial direction of $\sim 500$ pc from the Sun as shown in Fig. 4.

We analyse the vertical metallicity gradient for five mono-age stellar populations, defined as populations of stars in 2-3 Gyr age bins. In this paper, we use age intervals of 0-2, 2-4, 4-6, 6-8 and 8-11 Gyr and call them mono-age stellar populations. Strictly speaking, a 2-3 Gyr age bin cannot be called a mono-age stellar population. However, since we need a statistically significant sample of stars for each population, we use this large age bin. The cut in age at 11 Gyr is also due to the same sample size requirement (see Appendix A).

The vertical metallicity gradients of the mono-age stellar populations are analysed using the linear regression method of Kelly (2007), which employs a Bayesian hierarchical model. This approach enables us to measure the vertical metallicity gradient taking into account heteroscedastic measurement errors and intrinsic scatter in the regression relationship (Kelly, 2007). The equations to describe the probabilistic hierarchical model are given in Kelly (2007) as:

[TABLE]

where $(\xi_{\rm i},\eta_{\rm i})$ are the actual values. The observed values $(x_{\rm i},y_{\rm i})$ are measured with errors ($\epsilon_{\rm x,i},\epsilon_{\rm y,i}$) as:

[TABLE]

In this hierarchical model, $\xi_{\rm i}$ is drawn from a probability distribution described by a set of parameters $\psi$. $\eta_{\rm i}$ is drawn from a normal distribution with mean of $\alpha+\beta\xi_{\rm i}$ and variance $\sigma^{2}$, where $\alpha$, $\beta$, $\sigma$ are the slope, intercept and intrinsic dispersion, respectively. The observed values, $x_{\rm i}$ and $y_{\rm i}$ come from the multivariate normal density with mean $\xi_{\rm i}$ and $\eta_{\rm i}$ and covariance matrix $\Sigma_{\rm i}$, where $\Sigma_{\rm 11,i}=\sigma_{\rm y,i}^{2}$, $\Sigma_{\rm 22,i}=\sigma_{\rm x,i}^{2}$, where $\sigma_{\rm x,i},\sigma_{\rm y,i}$ are the measurement errors on $x$ and $y$ for the data point $\rm i$, respectively. The non-diagonal terms in $\Sigma_{\rm i}$ are zero since we assume non-correlated errors between metallicity and distance for our data. The method is based on deriving a likelihood function for the measured data $p(x_{\rm i},y_{\rm i}|\theta,\psi)=p(y_{\rm i}|x_{\rm i},\theta,\psi)p(x_{\rm i}|\psi)$, where $\rm\theta=(\alpha,\beta,\sigma^{2})$ (Kelly, 2007). In the current analysis, we model $\xi$ as a Gaussian distribution, characterised by a mean $\mu_{\rm 1}$ and variance $\tau_{\rm 1}$. These parameters define the parameter set $\psi$ (Kelly, 2007).

For each mono-age stellar population sample, we obtain the vertical height of stars $\left|Z\right|$ from a distance calculated from the TGAS parallaxes and Galactic latitudes, assuming the vertical offset of the Sun to be $Z_{\sun}=25$ pc (Bland-Hawthorn & Gerhard, 2016). The vertical distance errors are computed as $sin(b)\pi_{\rm e}/\pi^{2}$, where $\pi$ is the parallax and $\pi_{\rm e}$ is the error in parallax from the TGAS data, with a 0.3 mas systematic added in parallax (Gaia Collaboration et al., 2016b). The distance errors are not simply Gaussian with dispersion of $\pi_{\rm e}/\pi^{2}$ (e.g. Bailer-Jones, 2015). However, for simplicity we assume Gaussian errors to be a good approximation for our selected data with $\pi_{\rm e}/\pi<0.2$ and $\pi>0.0$. The metallicity and their associated errors are taken from the RAVE DR5 catalogue.

In the probabilistic model, we also consider the effect of the stellar population bias as suggested in Casagrande et al. (2016). These probabilities are determined as a function of stellar age ($\tau$), metallicity ([M/H]) and distance ($d$), assuming the Salpeter Initial Mass Function (Salpeter, 1955), using the approach of Casagrande et al. (2016) shown schematically in their Fig. 11. To avoid an extrapolation which has a risk to cause too extreme values, we use the probability associated with $\rm[M/H]=0.5$ and $\tau=13$ Gyr for those stars with $\rm[M/H]>0.5$ or $\tau>13$ Gyr, respectively. We then define a weight as $w_{\rm i}={p_{\rm i}({\tau_{\rm i},{\rm[M/H]_{i}},d_{\rm i}})}^{-1}$ for each star, which is then used to weight the log of the likelihood functions: $\frac{Nw_{i}{\rm i}}{\sum_{i}w_{\rm i}}\ln L_{\rm i}({\rm[M/H]}||Z|,\theta,\psi)$ and $\frac{Nw_{\rm i}}{\sum_{i}w_{\rm i}}\ln L_{\rm i}(|Z||\psi)$, where $\ln L_{\rm i}({\rm[M/H]})$ and $\ln L_{\rm i}(|Z|)$ correspond to $p(y_{\rm i}|x_{\rm i},\theta,\psi)$ and $p(x_{\rm i}|\psi)$, respectively. $N$ is the number of stars in the sample. This entropic weighting, i.e. weighting in the log for both likelihood functions is allowed since $\rm[M/H]$ and $\left|Z\right|$ are assumed to be independent measurements. For simplicity, we do not take into account the uncertainties associated with the weights, and instead use the weights from the measured age, [M/H] and distance, i.e. $1/\pi$, without considering their errors. To avoid an unreasonably large weight value due to a very low probability, which can happen if the measured values are significantly different from their actual value, we cap the weight at $w_{\rm i}<20.0$, corresponding to a probability $p_{\rm i}>0.05$. In this paper, we will refer to this weight as the stellar population weight and will perform an analysis with and without the weights. It is important to note that the stellar population weight does not correct the stellar population bias perfectly. The computed probabilities are dependent on the isochrones adopted in the stellar population model, and the current analysis does not take into account the observational errors of the assigned probabilities. Therefore, the results of the weighted case are intended to provide the effect of the stellar population bias in a qualitative manner.

Inference of the parameters is explored using a MCMC approach with an Adaptive Metropolis-Hastings algorithm, which allows us to sample posterior distributions and provides the full marginalised probability of the parameters. We implement this method in pymc, which facilitates modelling the complex dependencies of the priors employed in the current analysis. The prior space consists of the model parameters, $\alpha,\beta,\sigma$, as well as those describing the $\xi$ distribution, and a full account is given in Kelly (2007). They are drawn from the distributions considered in Kelly (2007) apart from $\mu$ and $\tau$, which are sampled from uniform distributions with more restrictive bounds.

3 Results

In Fig. 5, we present the vertical distribution in metallicity of the mono-age stellar populations and of the full sample in the whole age range of $0<{\rm age}\leq 11$ Gyr. As discussed in Section 2, in our selected colour-magnitude range each star has a different probability of being observed in the chosen colour and magnitude range depending on distance, age, and metallicity. We compute the stellar population weight as the inverse of this probability and colour-code the points in Fig. 5 by the stellar population weight. We apply the Bayesian hierarchical regression model described in Section 2 to the data for each mono-age stellar population. This method is employed for two cases, namely with and without the stellar population weight, referred to as the stellar population bias “weighted" and “unweighted" case, respectively. Fig. 5 also shows our best fit linear regression slope for the stellar population unweighted case as the median of the posterior distribution of the associated, marginalised MCMC chains. Lower and upper uncertainties correspond to the 16th and 84th quartile, respectively. Table 1 gives the best-fit values with uncertainties for the slope, $\alpha={\rm d[M/H]/d}\left|Z\right|$, intercept, $\beta$, and dispersion, $\sigma$, for the five mono-age populations as well as for the entire sample, including stars with ages greater than 11 Gyr.

The slope for each mono-age stellar population as a function of age is shown in Fig. 6. The results of the current analysis revealed that all our mono-age stellar populations show a negative vertical metallicity gradient. There appears to be a steady increase in steepness of the negative slope with increasing age. The youngest population of stars with ages between 0 and 2 Gyr show an almost flat vertical metallicity gradient, d[M/H]/d$\left|Z\right|$ = $-0.013^{0.066}_{0.070}$ dex kpc*-1*. The oldest population reach a significantly more negative value, at d[M/H]/d$\left|Z\right|$ = $-0.441^{0.122}_{0.125}$ dex kpc*-1* in the case of the stellar population unweighted case. The stellar population weighted case reveals a qualitatively similar trend of negatively steeper vertical metallicity gradient. The stellar population weighted case gives a slightly steeper slope than the unweighted case, especially for the relatively older populations, $\tau\gtrsim 4$ Gyr. As it can be seen in Fig. 5, at lower $\left|Z\right|$, the stellar population weight for the higher metallicity stars is higher than the stellar population weight for the lower metallicity stars. As a result, higher weights for higher metallicity stars lead to a steeper negative slope and higher values for the intercept (Fig. 7). Even with the added stellar population bias effect, the current results reveal an almost flat vertical metallicity gradient for the youngest population and increasingly negative vertical metallicity gradient for the older age populations. Hence, we conclude that the observed trend is robust and likely not due to the observational selection or stellar population biases.

Fig. 7 reveals the trend for our best-fit linear regression intercept values, $\beta$, of each mono-age stellar population for both the stellar population weighted and unweighted cases. The stellar population unweighted and weighted cases show a similar trend. The youngest stellar population of $0<\tau<2$ Gyr has a systematically higher metallicity than the older populations. The stars with ages between 2 and 8 Gyr show similar values for the intercept. The change in the intercept values is smaller for the stellar population weighted case. On the other hand, the oldest stellar population, $8<\tau<11$ Gyr, shows significantly lower metallicity. The intercept indicates the metallicity at the disc plane, $\left|Z\right|=0$ kpc. The stellar population weighted case shows systematically higher values of the metallicity at the plane for the stars older than $\sim 4$ Gyr. This trend corresponds to the steeper negative slope for the stellar population weighted case in Fig. 6. The steeper negative slope, driven by the higher weight of higher metallicity stars at lower $\left|Z\right|$, leads to the higher [M/H] at the disc plane for the stellar population weighted case. The similarity between the weighted and unweighted cases reinforces the robustness of our results. In fact, this trend is reminiscent to the age-metallicity relation observed for the Solar neighbourhood stars (e.g. Casagrande et al., 2011; Haywood et al., 2013; Anders et al., 2017). For example, Fig. 9 of Haywood et al. (2013) shows stars younger than 2 Gyr have higher [Fe/H] than older stars 222Note that Haywood et al. (2013) used [Fe/H]. On the other hand, we use the [M/H] value from RAVE DR5. Therefore, we use [Fe/H] when quoting the results of Haywood et al. (2013).. This systematically higher metallicity for stars younger than $\sim 2$ Gyr is also seen in fig. 4 of Anders et al. (2017). Haywood et al. (2013) also showed that stars older than 8 Gyr have lower [Fe/H] systematically, and identified them as the age-defined thick disc.

Finally, the values of the dispersion, $\sigma$, from our best-fit vertical metallicity gradient model for the mono-age stellar populations for both the stellar population weighted and unweighted cases are shown in Fig. 8. These values are almost the same for both the weighted and unweighted cases. Our results reveal that the intrinsic scatter increases steadily with age, from $\sigma\sim 0.13$ dex for the youngest population to $\sigma\sim 0.31$ dex for the oldest population. As discussed in Section 1, the larger intrinsic dispersion of the metallicity distribution of the older stellar populations can be explained by a radial mixing mechanism which brings more metal-rich (poor) stars from the inner (outer) region of the Galactic disc. The steady increase of the dispersion indicates that radial mixing has not reached saturation, which is expected to happen if the radial mixing processes mix stars born at different radii quickly. It is also important to note that our metallicity dispersion is systematically higher than what is obtained by Anders et al. (2017) for their mono-age stellar populations. This may be due to the large measurement uncertainties in [M/H] from the RAVE data, as compared to their metallicity from high signal-to-noise data from the APOGEE data. However, since there should not be a correlation between the age and metallicity errors, e.g. older stars should not have higher uncertainties in the RAVE data, we conclude the trend that we observe is not induced by the observational uncertainties or selection biases. Instead, it is likely to be due to an inherently higher dispersion in metallicity for the older population.

4 Discussion and Summary

In this paper, we measured the dependence of metallicity on vertical height for five different mono-age stellar populations by analysing combined data from the TGAS and RAVE DR5 catalogues with an advanced linear regression model. Our main results are summarised as follows:

The youngest population has an almost flat vertical metallicity gradient. Contrastingly, the older population shows a negative vertical metallicity gradient, i.e. lower [M/H] at a higher height. Furthermore, the mono-age stellar populations display an increasingly steeper negative vertical metallicity gradient with age. 2. 2.

From the intercept value from the linear regression model, which yields the metallicity at the disc plane (Z= 0 kpc), we found an age-metallicity relation. The youngest mono-age stellar population with ages between 0 and 2 Gyr has systematically higher metallicity than stars older than 2 Gyr. Across the three mono-age stellar populations within the age range of 2 to 8 Gyr, the metallicity at the disc plane remains almost constant. On the other hand, the oldest population of $8<\tau\leqslant 11$ Gyr shows significantly lower metallicity. 3. 3.

The linear regression model employed in our study enables us to analyse the intrinsic dispersion in the metallicity distribution. We find that the dispersion increases steadily with age up to 8 - 11 Gyr.

We minimised the observational selection bias induced by the various selection criteria by selecting a sample in a specific colour and magnitude range. Additionally, we applied the stellar population bias correction by weighting the likelihood function with the inverse of the probability of seeing a star having a particular age and metallicity and located at a certain distance in the selected colour and magnitude range. In a qualitative sense, the results obtained by taking into account the stellar population bias are similar to those obtained when no stellar population weighting was applied. This fact is evidence that the results of this analysis are robust, and unlikely induced by observational selection bias or stellar population bias.

The vertical metallicity gradients found in our study are qualitatively consistent with the study by Xiang et al. (2015) despite their vertical range (0 < $\left|Z\right|$ < 2 kpc) being much more extended than ours. At the Solar radius, their results reveal a steeper negative vertical metallicity gradient for the older mono-age stellar populations, which is in qualitative agreement with our results for stars with ages between 2 and 11 Gyr. Xiang et al. (2015) also found that stars with ages between 11 and 16 Gyr have a flatter gradient. Unfortunately, we cannot measure the vertical metallicity gradient of this mono-age population since we do not have a statistically significant subset of stars older than 11 Gyr in our sample. The bright magnitude limit and the narrow colour range employed in our selection, as well as the limited volume coverage of our sample, i.e. we are still rather confined to the plane, where the majority of stars are young, limit our sample to stars no older than 11 Gyr. At the youngest end, we found that stars with ages between 0 and 2 Gyr have an almost flat vertical metallicity gradient. Since Xiang et al. (2015) did not show the vertical metallicity measurement for stars younger than 2 Gyr, our results reveal for the first time the flat vertical metallicity gradient of the youngest population, and connect it with the gradually steeper negative vertical metallicity gradient of the older population.

Before discussing the various scenarios that can explain our results, it is important to consider the Galactic disc evolution scenarios that our results can reject. The fact that the older population has a distinctly negative vertical metallicity gradient and the gradient becomes steeper with age strongly disfavors the existence of a very strong vertical phase mixing of stars (like strong heating by satellite accretion). Such a process would have driven a flat vertical metallicity gradient across the mono-age stellar populations formed before the heating happened. In the light of our results, the Milky Way disc should not have undergone dramatic heating leading to strong vertical mixing. Furthermore, they seem to support the recent claims that the Milky Way has not undergone any major merger in the last 10 Gyr (Ruchti et al., 2015; Casagrande et al., 2016).

As justified in Section 1, it is natural to assume that each mono-age disc population had a negative radial metallicity gradient and no vertical metallicity gradient when they formed. Under these conditions, our observed trend challenges one of the numerical results of Kawata et al. (2017). According to their study, if a mono-age stellar population of the Galactic disc had a constant scale height at all radii, radial mixing would have driven a positive vertical metallicity gradient as discussed in Section 1. Our results invalidate constant scale height for the initial structure of the mono-age disc populations and suggest that the Galactic disc stars should not have had a star-forming disc with a constant scale height across all radii.

On the other hand, our results support another model of Kawata et al. (2017) where they consider that each mono-age stellar population formed from a flaring star-forming disc. Flaring star-forming regions have been observed previously in the Milky Way in the distribution of the HI and molecular gas (e.g. Scoville et al., 1993; Nakanishi & Sofue, 2016), and young stellar populations traced by Cepheids (Feast et al., 2014). It is not impossible to expect that all the thin disc stars formed in a flaring star-forming region (e.g. Rahimi et al., 2014; Minchev et al., 2015). Kawata et al. (2017) demonstrate that if this flaring thin disc has a negative radial metallicity gradient when the mono-age stellar populations formed, then radial mixing processes could lead to a negative vertical metallicity gradient as discussed in Section 1. If radial mixing was a slow process, then this scenario could explain our observational trend.

Based on the same flaring disc scenario, if the star-forming region had a similar level of flaring across all mono-age stellar populations, but a steeper negative d[M/H]/dR at an earlier epoch, radial mixing could drive the steeper vertical metallicity gradient for the older mono-age stellar populations. However, Anders et al. (2017) found an opposite trend in the radial metallicity gradient of the current disc, with stars with ages between 1 and 4 Gyr having d[M/H]/d$R$ $\simeq-$ 0.07 dex kpc*-1*, whereas the older population shows a flatter gradient of d[M/H]/d$R$ $\simeq-$ 0.03 dex kpc*-1*. A similar trend was observed by Casagrande et al. (2011) when they analysed radial metallicity gradients in the Solar neighbourhood using the guiding radius of stars. We still do not know whether the flatter radial metallicity gradient for the old disc is due to the radial metallicity gradient of the star-forming disc being flatter at an earlier epoch, or due to radial mixing flattening the radial metallicity gradient more prevalently for the older disc. Thus, this raises the question about what was the initial radial metallicity gradient of each generation of stars in the Galactic disc. An extensive set of stars with accurate measurements of age and metallicity spanning a large vertical and radial range should provide stronger constraints on the initial metallicity distribution of the Galactic disc at different epochs and the strength of radial mixing effects by comparing with numerical simulations. Such data will be available soon with $Gaia$ DR2 complemented by high-resolution spectroscopic survey data, such as APOGEE-2, $Gaia$-ESO and GALAH, and asteroseismic age information from the K2 campaign (e.g. Stello et al., 2015) and NASA’s TESS mission, and ultimately by ESA’S Plato mission.

Acknowledgements

IC, DK and MC acknowledges the support of the UK’s Science & Technology Facilities Council (STFC Grant ST/K000977/1 and ST/N000811/1). IC is also grateful the STFC Doctoral Training Partnerships Grant (ST/N504488/1). LC is supported by Australian Research Council Future Fellowship FT160100402.

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

Funding for RAVE has been provided by the Australian Astronomical Observatory; the Leibniz-Institut fuer Astrophysik Potsdam (AIP); the Australian National University; the Australian Research Council; the French National Research Agency; the German Research Foundation; the European Research Council (ERC-StG 240271 Galactica); the Istituto Nazionale di Astrofisica at Padova; The Johns Hopkins University; the National Science Foundation of the USA (AST-0908326); the W. M. Keck foundation; the Macquarie University; the Netherlands Research School for Astronomy; the Natural Sciences and Engineering Research Council of Canada; the Slovenian Research Agency; the Swiss National Science Foundation; the Science & Technology Facilities Council of the UK; Opticon; Strasbourg Observatory and the Universities of Groningen, Heidelberg and Sydney. The RAVE web site is at http://www.ravesurvey.org

We thank Jo Bovy for making $\tt GALPY$ publicly available (Bovy, 2015) and to Željko Ivezić, Andrew J. Connolly, Jacob T. VanderPlas and Alexander Gray for doing the same with $\tt astroML$ (Vanderplas et al., 2012). We employed plotting routines from both these packages.

Appendix A Comparison analysis with PJM2017

In this paper, we used the age estimate obtained using the isochrone fitting method $\tt Elli$ (Lin et al., 2017). To test if our results are not affected by choice of the technique employed to get age estimates, we have applied the same analysis to the age estimates provided in McMillan et al. (2017), hereafter $\tt PJM17$, who used a different isochrone fitting method and different isochrones for the TGAS-RAVE sample. First, we compare the ages derived with $\tt Elli$ and $\tt PJM17$ in Fig. 9. The $\tt PJM17$ method provides overall larger derived ages than $\tt Elli$ and the difference between these two age estimates is more significant in the older population.

We applied the same analysis described in Sec. 2 to the dataset using the $\tt PJM17$ age estimates and their uncertainties. Furthermore, we also cut in relative age error at less than 50 %. A subsequent MCMC approach is taken to characterize the slope, intercept and dispersion of the vertical metallicity profile of the mono-age populations considered in the main text. For this comparison, we also discuss the results obtained for stars older than 11 Gyr because $\tt PJM17$ provides a statistically significant sample of these stars. The last age bin contains stars older than 11 Gyr in Fig. 11, but for visualisation, we plot the results at age=14 Gyr. We also applied the distance estimated in $\tt PJM17$ rather than the distance from the TGAS parallax to show that the results are not affected by choice of distance estimate. We call this sample the “PJM2017 sample" and the sample selected in Sec. 2 the “$Elli$ sample". Fractional age uncertainty as a function of age for both samples is shown in Fig. 10. Both samples are comparable in terms of the distribution of age uncertainties against age.

Fig. 11 shows the results of the stellar population bias weighted (noted as “sp-weighted") and unweighted (noted as “unweighted") cases. Within their estimated uncertainties, the results for the gradient of the metallicity as a function of height for the unweighted case are consistent between these analyses, except the age bins of 0-2 Gyr and 8-11 Gyr where the $\tt PJM17$ sample shows rather unexpectedly steeper and flat gradient, respectively. The results from the $\tt PJM17$ sample preserve the main conclusion of this paper qualitatively, namely that older stars have a more negative vertical metallicity gradient than younger ones, with the stellar population weighted result showing it more clearly.

Interestingly, the results for the slope for the 8<$\tau$<11 Gyr age bin found from our analysis using $\tt Elli$ are similar to those obtained for the >11 Gyr age bin in the $\tt PJM17$ sample. This trend is present across all ages considered, with the results for the vertical metallicity gradient obtained for the $\tt Elli$ sample being comparable to those for the $\tt PJM17$ sample in one age bin older. In other words, we can argue that although the age estimate may suffer from the systematic uncertainty as demonstrated in Fig. 9, we can trust the relative trend of vertical metallicity gradient with respect to stellar age. Hence, this comparison shows that the qualitative trend we obtained from the analysis using the $\tt Elli$ sample is not affected by the systematic biases or errors in the isochrone age estimates.

Fig. 11 also provides the results taking into account the age uncertainties by computing the weights of the contribution to each age bin for each star following Anders et al. (2017) and assuming the age is drawn from a Gaussian distribution. We also weigh the age weights by the stellar population described in Section LABEL::method. The results are labeled as “age-sp weighted" in Fig. 11. We applied the same weighting procedure for the logarithm of the likelihood function as we did for the stellar population bias weights in Sec. 2. Fig. 11 shows that the slopes in the age-sp-weighted case are systematically flatter as stars with different mean ages with large errors contribute to adjacent age bins. Therefore, we think that the age-weighted results are underestimating the slope. Nevertheless, there is still a qualitative trend of more negative slopes in older mono-age populations. This also provides the robustness of our conclusions. Fig. 11 also presents the results of the intercept and the dispersion. The trends are generally consistent with the results of the slope, and again this demonstrates the robustness of our conclusion.