The Lives of Stars: Insights from the TGAS–RAVE–LAMOST Data Set

The study of the history of the Milky Way Galaxy, galactic archaeology, relies on fossil evidence to draw conclusions about the formation processes that resulted in the Milky Way we see today. Some of this evidence can be seen directly in stellar number counts, as is the case for stellar streams formed by the infall of both large galaxies (e.g., the Sagittarius dwarf spheroidal; Ibata et al. 1994; Hernitschek et al. 2017) and small systems (e.g., the tidal tails around the globular cluster Palomar-5; Odenkirchen et al. 2001; Erkal et al. 2017). Such evidence provides direct observational constraints on the extent to which the Milky Way is formed in situ (Eggen et al. 1962), and which portion is formed via accretion (Searle & Zinn 1978; Helmi et al. 2017).

Evidence of formation mechanisms can be seen indirectly in the form of the Galactic warp (Djorgovski & Sosin 1989; Schönrich & Dehnen 2017) and flare (Kent et al. 1991; Feast et al. 2014), which could be caused either by accretion or by internal disk instabilities. Internal disk instabilities could also give rise to features such as the Galactic bar (Combes & Sanders 1981; Weinberg 1992) and the spiral arms (which are perhaps long-lived or perhaps transient; for discussion, see, for example, Sellwood 2011 or Martinez-Medina et al. 2017).

In the era of photometric surveys, it was difficult to draw sweeping conclusions about these broad features, but the more one looked, the more evidence pointed to a fairly chaotic formation (see for example the lumpy Milky Way topography in Newberg et al. 2002, or the complex Magellanic system in Belokurov & Koposov 2016). Estimates of kinematics could be made based on two-dimensional astrometric proper motions; however, until recently, these have been defined by long-baseline observations cross-matching the most modern photometric observations with 20th century photographic plate surveys—see for example the UCAC (Zacharias et al. 2013) and PPMXL (Roeser et al. 2010) catalogs—which generally achieve precisions of ∼3 mas yr⁻¹, or 15 km s⁻¹ kpc⁻¹ (an inconvenient limit for discussing features several kiloparsecs away). Other problems also occur in the confusion of cross-matching large numbers of surveys with varying specifications, leading to some discrepancies between proper motion catalogs (Vickers et al. 2015), systematics within individual catalogs (Wu et al. 2011; Vickers et al. 2016), and, for some, errors clustered around the plate-scale, indicating the sub-optimal condition of cross-matching plates with CCD data (Pearl et al. 2017).

In the early 2000s, a variety of large-scale spectroscopic surveys began producing data:

• 1.  
  the Radial Velocity Experiment released almost 25,000 spectra in 2006 (RAVE; currently ∼520,000 spectra available; Steinmetz et al. 2006; Kunder et al. 2017);
• 2.  
  the Sloan Digital Sky Survey (SDSS; York et al. 2000) project the "Sloan Extension for Galactic Understanding and Exploration" (SEGUE; Yanny et al. 2009) released spectra of 240,000 stars in 2009 (SEGUE-2 added another ∼120,000);
• 3.  
  the Large Area Multi-Object Spectroscopic Telescope (LAMOST; Luo et al. 2015) released 1.5 million spectra in 2013 (currently ∼7.5 million spectra are available); and
• 4.  
  the Apache Point Observatory Galactic Evolution Experiment (APOGEE; Holtzman et al. 2015) released spectra of 150,000 stars in 2015, and APOGEE-2 is currently surveying a further 300,000.

With this outpouring of data, the parameter spaces of radial velocity and abundance became available for large numbers of stars. With radial velocities, overdensities could be tested for coherence at distances impossible with extant proper motions (a nice example of this is seen in the 34 kpc distant Cetus Polar stream of Newberg et al. 2009), and the potential of the Milky Way could be tested out to tens of kiloparsecs (as done with blue horizontal branch stars by Xue et al. 2008).

With abundances, work could be done by comparing elemental abundance ratios to provide indirect insights into the ages of populations. The ratio of iron and alpha-element abundances particularly holds information about star formation rates of populations (see for example Wheeler et al. 1989 and references therein) and is useful for exploring differences between spatially mixed populations, such as those of the canonical thin and thick disks (Gilmore & Reid 1983), which are thought to have different formation histories and thus chemistry (see Fuhrmann 1998, and later Bovy et al. 2012 and Bovy & Rix 2013). Alpha elements, being primarily produced in swift, Type-II supernovae, are generally thought to be indicative of regions of high star formation, while iron is injected into the interstellar medium (ISM) by Type-Ia supernovae (SNe Ia), which occur in all disk regions at a fairly constant rate. So stars with high alpha abundances relative to iron, such as thick-disk stars, were probably born in a quickly evolving environment of high star formation; and the alpha abundance relative to iron gradually falls for all subsequent generations of stars born into a more quiescent Milky Way that is being slowly enriched in iron by SNe Ia. More complicated "chemical tagging" (see Freeman & Bland-Hawthorn 2002), utilizing a wider variety of elemental abundances arising from different nuclear processes, can provide even more information, and is explored early in F and G disk stars in Bensby et al. (2003, 2005) or more recently in, e.g., Ting et al. (2016) in APOGEE or Kos et al. (2018) in GALAH (Galactic Archaeology with HERMES, the High Efficiency and Resolution Multi-element Spectrograph; De Silva et al. 2015). Recently, chemodynamical models of Milky Way-like galaxies have been crucial to understanding these abundance observations in the context of cosmological galaxy evolution with effects such as inside-out formation scenarios (Pilkington et al. 2012), migration (Minchev et al. 2014; Kubryk et al. 2015a, 2015b), spiral arm influence (Grand et al. 2016), and realistic galactic environments (Ma et al. 2016).

However, while many proxies for age existed—such as color (with extremely blue colors generally singling out stars with very short lifespans), color–apparent magnitude (with coherent populations following isochrones in this space), and alpha–iron ratios (making certain assumptions about the formation mechanisms leading to those ratios)—precise ages via isochrone matching remained out of reach for most stars until accurate distances could be obtained. Age estimates may be made with atmospheric parameters alone, but the precision is generally low in large spectroscopic surveys.

Some other specialized methods for measuring age exist, such as elemental abundance ratios in giants (Martig et al. 2016) or astroseismic age estimation for giants (Casagrande et al. 2016), but robust isochrone-matched age determination for the full range of stellar types had to wait until more accurate distance information became available.

With the Hipparcos mission (Perryman et al. 1997; van Leeuwen 2007) providing parallaxes for 2.5 million stars in the form of the Tycho-2 astrometric catalog (Høg et al. 2000), the way was open to directly calculate the ages and kinematics for huge numbers of stars, as was done by the Geneva–Copenhagen Survey (Nordström et al. 2004) for about 14,000 stars. Similar endeavors were undertaken using the high-accuracy radial velocity Planetary Searcher (Mayor et al. 2003) for about 1000 stars (see for example Haywood et al. 2013).

The Gaia space telescope (Gaia Collaboration et al. 2016b), launched in 2013, will drastically improve the current astrometric precision on parallaxes (with microarcsecond parallaxes expected out to 20th mag) and proper motions (estimated to be two to three orders of magnitude more precise than current large, ground-based surveys). This will greatly increase the sample of stars for which six-dimensional phase space information, and ages, can be accurately calculated by cross-matching with existing spectroscopic surveys. The mission will also provide ∼150 million radial velocities, atmospheric parameters for ∼5 million stars, and detailed abundances for ∼2 million stars (according to the science performance page¹ ; see also Recio-Blanco et al. 2016).

This paper presents an analysis of the kinematic and abundance trends as a function of age for a large sample of stars in a 1 kpc wide sphere around the Sun, taking advantage of the first Gaia data release, the fourth LAMOST data release, and the fifth RAVE data release.

We will consider the behavior of stars in different age groups with a specific interest in stellar migration. Stellar migration may generally be thought of as consisting of two mechanisms: (1) the gradual heating of stars off of circular orbits and onto eccentric orbits, extending the radial range they inhabit but conserving angular momentum (so they travel quickly at perigalacticon and slowly at apogalacticon), which causes populations to be observed over larger ranges of radii over time, an effect called "blurring"; and (2) the exchange of angular momentum between individual stars and non-axisymmetric elements of the potential, such as the spiral arms or bar, allowing stars to change their orbital radii while maintaining near-circular orbits, an effect known as "churning." Both of these mechanisms of stellar migration tend to flatten metallicity gradients. This nomenclature is taken from Schönrich & Binney (2009).

Stellar mixing, or the net movement of stars away from their birth radii, has been widely investigated in simulations in recent years. Mixing was investigated as an effect of angular momentum exchange over the corotation radii of spiral arms (by e.g., Sellwood & Binney 2002; Roškar et al. 2012; Grand et al. 2014) and was found to be responsible for large changes in individual stars' angular momenta, sometimes as large as 50%. Roškar et al. (2008) looked into the effects of this process on observed age and abundance and found that this type of migration typically flattened the age–metallicity relationship. Minchev & Famaey (2010) expanded on this with the addition of a bar component and found that bars may affect the change in angular momentum in the vicinity of the corotation radius and the outer Lindblad resonance and also that bars may enhance the rate at which metallicity gradients flatten (Minchev et al. 2011). Further, Vera-Ciro et al. (2014), using a disk with spiral arms embedded in a dark matter halo, showed that churning was most efficient for stars with low vertical velocity dispersions.

In Section 2 we outline the data used and describe corrections and cuts needed to clean them. In Section 3 we describe the process for estimating the ages and age errors for the stars. Section 4 details the coordinates and orbit calculations. In Section 5 we focus on the dynamics, discussing: the vertex deviation angle, the rate of increase of velocity dispersion and velocity dispersion ratios, the Oort constants and rotation curve as a function of stellar age, and some strange behavior of the youngest stars in our sample. Section 6 folds in abundance information to investigate the chemical evolution, looking at: inside-out formation, the relative rates of enrichment for the inner disk and outer disk, blurring, and possible observational evidence for outward stellar churning. In Section 7 we conclude. An appendix is devoted to dealing with the selection functions of spectroscopic surveys; and another details some offsets and corrections between the LAMOST and RAVE data sets.

The astrometric data used are the Tycho-Gaia Astrometric Solution (TGAS; Gaia Collaboration 2016a). This data set is a combination of the Tycho-2 astrometric catalog and the astrometric data collected in the first year of the Gaia satellite's operations. From these data we make use of the parallaxes (∼0.3 mas uncertainties), the proper motions (∼1 mas yr⁻¹ uncertainties), and the mean G-band magnitude (≲0.03 mag uncertainties for our objects).

It has been noted by many that the transformation from parallax space into distance space d is nonlinear and therefore can dramatically skew the error. Bailer-Jones (2015) showed that this effect will become pronounced after fractional parallax errors of 20%. Schönrich & Aumer (2017) have noted that the distance estimates of TGAS are relatively consistent with the predictions of their internally consistent galactic model up to a fractional parallax error of 20%. To correct for this it is necessary to adopt a prior based on estimated density distributions, as in Astraatmadja & Bailer-Jones (2016), or on density and velocity distributions, as in Schönrich & Aumer (2017). We prefer to work here without a prior since we will be limiting ourselves to high-precision data anyway, to get the most accurate age estimations. In light of this, we adopt 20% parallax errors as a hard cut. In the future, however, as Bailer-Jones (2015) notes, roughly 80% of the Gaia catalog will be beyond this 20% error threshold, so it is worth considering, now, the types of priors that will be acceptable to use to correct for this.

The spectroscopic data used come from two sources: LAMOST DR4v1 and RAVE DR5.

The LAMOST spectrographs collect R ∼ 1800 spectra in the wavelength range from 370 to 900 nm. As of DR4, the LAMOST data set contains over 4 million spectra above δ ∼ −10°, with a particular spatial concentration in the direction of the Galactic anticenter. Of these, 70% have signal-to-noise ratios (S/Ns) above 20 in the g band. The median errors are δT_eff ∼ 139.5 K, δ log(g) ∼0.46 dex, δ[Fe/H] ∼ 0.17 dex.

We cross-match these LAMOST spectroscopic data with the TGAS data within 10'' and we find that about 75% and 97% of the matches are within 1'' and 4'', respectively. We compare the (B − V) Tycho-2 colors with the LAMOST effective temperatures and find that the cross-matches out to 4'' are not of a worse overall qualitative match than the <1'' sample, so we choose 4'' as the matching radius.

We further discard all objects with failed error estimates in T_eff, log(g), or [Fe/H] (meaning pipeline failures, we require errors >0), and objects with S/Ns in g, r, or i less than 20. A final cut requiring positive parallaxes leaves us with ∼73,000 objects whose median errors are δT_eff ∼ 83.6 K, δ log(g) ∼ 0.40 dex, δ[Fe/H] ∼ 0.10 dex. However, it should be noted that the pipeline errors from the LAMOST data products are known to be flawed (see for example Schönrich & Aumer 2017). We change the radial velocity errors to be in line with values expected from comparisons of the LAMOST kinematics to models (7.2 km s⁻¹, Schönrich & Aumer 2017). We leave the gravity, temperature, and [Fe/H] error values at their pipeline values, which are likely overestimated (based on comparison with RAVE measurements and error distributions for a population with observations existing in both surveys).

The RAVE spectrograph observes at a higher resolution (R ∼ 7500) in a smaller wavelength range (841–879.5 nm) around the calcium triplet lines. These observations are arranged relatively evenly in pencil beams south of δ ∼ 0°, excepting in the plane of the Galaxy. The RAVE data product boasts more precise pipeline error medians than LAMOST—δT_eff ∼ 74 K, δ log(g) ∼ 0.16 dex, δ[M/H] ∼ 0.10 dex.

Rather than cross-matching the TGAS data to the RAVE data ourselves, we utilize the helpfully included cross-match that is packaged with RAVE DR5. We place similar restrictions on the data set: a positive parallax, total parallax error <20%, S/N >20, non-zero errors on T_eff, log(g), and [M/H] measurements. We additionally require the flags "c1," "c2," and "c3" to be "n" for normal. This yields a sample of about 67,000 objects with atmospheric parameter errors around δT_eff ∼ 76 K, δ log(g) ∼ 0.14 dex, δ[M/H] ∼ 0.10 dex.

By using data from these two surveys, we essentially double our spectroscopic sample that overlaps with the TGAS data, and also gain almost-full-sky coverage rather than the half-sky coverage offered by either survey independently (see Figure 1). However, it should be noted that this creates an extremely complicated selection function, with RAVE generally observing brighter, cooler objects more oriented toward the center of the Galaxy than LAMOST. Additional selection complexity is added with respect to the Gaia to Tycho-2 matching, and the subsequent matching of that catalog to LAMOST and RAVE. We have endeavored to address this in as comprehensive a way as possible in Appendix A.

Zoom In Zoom Out Reset image size

Duplicate observations exist in the data; they are caused by overlaps between the RAVE and LAMOST targets, and also by intentionally repeated observations. We remove these duplicates from our sample by selecting only the objects with the lowest errors in atmospheric temperature, gravity, and metallicity (preferentially in that order) within a 7'' circle on the sky.

We note that the RAVE [M/H] measurement is not wholly interchangeable with the LAMOST [Fe/H] estimate and also that the radial velocity estimates are systemically offset between the two surveys. These issues are addressed in Appendix B.

We estimate the ages of our stars using a Bayesian method similar to that outlined in Jørgensen & Lindegren (2005).

An isochrone grid is constructed using the Padova isochrones² (see Bressan et al. 2012; Chen et al. 2014, 2015; Tang et al. 2014). The grid is evenly spaced in [M/H] every tenth of a dex from −1.5 to 0.3 dex. The grid is non-uniformly spaced in age: every 0.01 Gyr from 0.03 to 0.1 Gyr; every 0.05 Gyr from 0.1 Gyr to 0.5 Gyr; every 0.1 Gyr from 0.5 to 1 Gyr; every 0.25 Gyr from 1 to 4 Gyr; every 0.5 Gyr from 4 to 13 Gyr. We space the grid in this uneven manner because isochrone morphology changes with age less at older ages than at younger ages.

For each star in our sample, we use the parallax, parallax error, and observational g-band Gaia magnitude (G) from the TGAS data set. If the object is in the RAVE data set, we use the stellar parameter pipeline values of [M/H], log(g), and T_eff; if the object is in LAMOST, we use the estimated [M/H] from Appendix B and the pipeline values of log(g) and T_eff.

The isochrones are converted into parallax space based on the individual star for which the age is currently being estimated by considering the observed G magnitude of the star and the absolute G magnitude of the isochrone:

$$\begin{eqnarray}&&{\pi }_{{\rm{iso}}}=100\times {10}^{({G}_{{\rm{abs}}}-{G}_{{\rm{obs}}})/5}.\end{eqnarray} \tag{ 1 }$$

In the paper, when considering the errors on the parallaxes, we do not, as is sometimes seen, add in a systematic error of 0.3 mas. This is following the recommendation of Brown (2017).

Each isochrone consists of n points, such that i = n − 1 sequential line segments may be constructed. χ² distances to each of these i line segments are calculated:

$$\begin{eqnarray}&&{\chi }_{i}^{2}=\displaystyle \sum _{j}{\left(\displaystyle \frac{{q}_{j}^{\mathrm{obs}}-{q}_{i,j}^{\mathrm{iso}}}{{\sigma }_{j}^{\mathrm{obs}}}\right)}^{2},\end{eqnarray} \tag{ 2 }$$

where j indicates the parameters (π, [M/H], log(g), and T_eff); q^obs and σ^obs are the measurement and measurement errors of the parameter; and ${q}_{i}^{\mathrm{iso}}$ corresponds to the parameter values of the closest point on the ith line segment connecting points i and i + 1 in the isochrone. The χ² values are then converted to likelihoods as

$$\begin{eqnarray}&&{L}_{i}=\displaystyle \prod _{j}\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{j}}\times \exp \left(-\displaystyle \frac{{\chi }_{i}^{2}}{2}\right).\end{eqnarray} \tag{ 3 }$$

We further weight the likelihoods by an initial mass function of the line segment i of the isochrone, so

$$\begin{eqnarray}&&{L}_{{\rm{iso}}}={\int }_{M}{L}_{i}\times {\rm{IMF}}(M)\,{dM},\end{eqnarray} \tag{ 4 }$$

such that

$$\begin{eqnarray}&&{L}_{\mathrm{iso}}=\displaystyle \sum _{i}{L}_{i}^{{\prime} },\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}&&{L}_{i}^{{\rm{{\prime} }}}={L}_{i}\times {\rm{IMF}}\left(\displaystyle \frac{{M}_{i}+{M}_{i+1}}{2}\right)\times ({M}_{i+1}-{M}_{i}).\end{eqnarray} \tag{ 6 }$$

Here, M_i+1 and M_i are the initial masses of the isochrone points bounding isochrone line segment i and

$$\begin{eqnarray}&&\mathrm{IMF}(M)\propto {M}^{-\alpha }\end{eqnarray} \tag{ 7 }$$

is the Kroupa (2001) initial mass function. α in this case is

$$\begin{eqnarray}{\alpha }({M})=\left\{\begin{array}{ll}0.3 & {\rm{f}}{\rm{o}}{\rm{r}}\,{M}\lt 0.08\,{{M}}_{\odot }\\ 1.3 & {\rm{f}}{\rm{o}}{\rm{r}}\,0.08\,{{M}}_{\odot }\lt {M}\lt 0.5\,{{M}}_{\odot }\\ 2.3 & {\rm{f}}{\rm{o}}{\rm{r}}\,{M}\gt 0.5\,{{M}}_{\odot }.\end{array}\right.\end{eqnarray} \tag{ 8 }$$

Marginalizing over the metallicity, we find the posterior probability of a star being age τ to be

$$\begin{eqnarray}&&G(\tau )={\int }_{{\rm{all}}[{\rm{M}}/{\rm{H}}]}{L}_{{\rm{iso}}}(\tau ,[{\rm{M}}/{\rm{H}}])\,d[{\rm{M}}/{\rm{H}}].\end{eqnarray} \tag{ 9 }$$

The 68% confidence intervals on the ages are estimated, as suggested in Jørgensen & Lindegren (2005), by finding the extent of the G(τ) function greater than 0.61 times the maximum value of G(τ).

At this point, we implement some additional data quality cuts. In our initial age distribution, we have excesses at the edge cases of our data; the reason for this is apparent from Figure 2. This Hess diagram depicts the observational sample, along with the "reddest" and "bluest" isochrones in our grid. It is apparent that a non-negligible portion of the observed stars lie outside the search space, for example in the crook redward of the old isochrone turnoff and the high log(g) observations below the main sequence. This is not especially surprising, considering that the surface gravity errors of σlog(g) ∼ 0.3 dex and temperature errors of σT_eff ∼ 100 K would naturally scatter the observations off their true locus. The stars outside this search space will be placed into the edge-case bins, overpopulating them. In light of this, we further reject any star that is 3σ inconsistent with the nearest point on the isochrone grid. A sample of about 105,000 objects remains. Rejecting the youngest (and oldest) edge-case age bins leaves ∼98,000 (∼87,000) stars.

Zoom In Zoom Out Reset image size

Further, we find that the weighting scheme implemented in Appendix A produces weights that follow a near-normal distribution in log-space. The weights vary by a few orders of magnitude and we wish to prevent highly weighted stars from completely overpowering the rest of the sample. Based on these two observations, we cut the data such that the allowed weighting is within three standard deviations of the mean of the calculated weight values (in logarithmic space). This leaves ∼55,000 stars.

We briefly explore the effects of extinction on our sample by applying the full column density extinctions of Schlegel et al. (1998), with the extinction coefficient of A_G = 2.55 × E(B − V) from Belokurov et al. (2017), to our data and recalculating the ages. We find that τ_observed − τ_dereddened = 0.08 ± 1.89 Gyr. Since the systematic difference is so small, we feel it safe to neglect extinction in our results.

The coordinates used in this work are all right-handed. The Sun is located at (X, Y, Z) = (−8.27, 0, 0) kpc and rotates at (V_R, V_ϕ, V_Z) = (0, − 236, 0) km s⁻¹ (where V_R is positive outward)—or equivalently (U, V, W) = (0, 236, 0) km s⁻¹. The solar motion with respect to the local standard of rest is (U, V, W) = (13.0, 12.24, 7.24) km s⁻¹. These values are all taken from Schönrich (2012) and Schönrich & Aumer (2017).

The observables are transformed into XYZ space using the naive distance estimation of d =  1/π and into UVW space using the method described in Johnson & Soderblom (1987).³ The average values of, and errors on, these derived quantities are calculated by performing 100 random samples of the data scattered about the errors in radial velocity, proper motion, and parallax. In this paper when we refer to velocities or positions, it is to the average values of these Monte Carlo samples.

Orbits are calculated using the galpy code described in Bovy (2015).⁴ The Milky Way potential used is the included MWPotential2014 (see Section 3.5 of Bovy 2015), which consists of three components: a power-law bulge that is exponentially cut off (McMillan 2011), a Miyamoto–Nagai disk (Miyamoto & Nagai 1975), and an NFW halo (Navarro et al. 1997). The default potential is modified to be consistent with our adopted location, rotation speed, and peculiar solar velocity relative to the local standard of rest. Orbits are calculated on time steps such that the total change in energy is negligible over the course of the orbit (the relative change in energy is less than 1 × 10⁻⁹ for all but three of our orbits).

5.1. Vertex Angle and Velocity Dispersion

In an axisymmetric potential, it is expected that stars will be born on near-circular orbits; as they are heated onto more eccentric, epicyclic orbits, the velocity plane of V_ϕ versus V_r will be symmetric about V_r. The Milky Way potential, however, is not axisymmetric, so the effects of the bar and the spiral arms can be quite influential in the structure of the velocity space, and the velocity plane is consequently asymmetric.

For example, the well known Hercules stream may be stars trapped in orbits around Lagrange points of the bar (Pérez-Villegas et al. 2017) or, alternatively, stars in the outer Lindblad resonance (Hunt et al. 2018). The outer Lindblad resonance of the bar is also predicted to create observable deviations and bimodalities in the velocity plane (Dehnen 2000).

Besides the bar, the relative locations of the spiral arms are thought to affect the orientation of the ridgeline of the velocity plane (Antoja et al. 2011). The deviation of the velocity vertex (Dehnen & Binney 1998, Equation (23)),

$$\begin{eqnarray}&&{{\ell }}_{v}\equiv \displaystyle \frac{1}{2}\arctan \left(\displaystyle \frac{2{\sigma }_{{xy}}^{2}}{{\sigma }_{{xx}}^{2}-{\sigma }_{{yy}}^{2}}\right),\end{eqnarray} \tag{ 10 }$$

where ${\sigma }_{{ij}}^{2}$ are elements of the velocity tensor,

$$\begin{eqnarray}&&\langle ({V}_{i}-\langle {V}_{i}\rangle )\cdot ({V}_{j}-\langle {V}_{j}\rangle )\rangle ,\end{eqnarray} \tag{ 11 }$$

describes this global tilt in the velocity plane.

Stars are generally thought to be born on cool, near-circular orbits with very little velocity dispersion relative to their siblings. Over time, they will be heated by interactions with molecular clouds (see, for example, Aumer et al. 2016; Gustafsson et al. 2016) and scattered by interactions with the non-axisymmetric potential, and so their velocity dispersions will increase over time. Observationally, this means that younger stars will tend to be more concentrated in the velocity plane, while older stars will tend to be more diffuse.

In Figure 3 we plot the velocity dispersions, ratios thereof, the deviation of the vertex angle, and the rotational velocity curve. The velocity dispersions have the observational errors subtracted, and the vertex angle is corrected for the observational covariances as in Smith et al. (2009). We see the characteristic rise in all components of the velocity dispersion as a function of age, an expected result of heating over time (or, possibly, of the cooling of the source gas of the disk stars over time, as explored in Aumer et al. 2016). As the dispersions rise, we see the characteristic fall in the rotation speed, V_ϕ, that is associated with the asymmetric drift. V_R (not shown) does not change significantly with age, the average being −3.9 km s⁻¹, comparable to Williams et al. (2013).

Zoom In Zoom Out Reset image size

The dispersion ratios indicate the relative importance of various heating mechanisms. The radial and azimuthal velocity dispersions are thought to be mostly increased by the spiral arms (Sellwood & Carlberg 1984), and so we expect them to rise proportionally—as is seen by the flat black curve in the bottom right panel. The vertical velocity dispersions are thought to arise from interactions with molecular clouds in tandem with the spiral arms (Lacey 1984; Jenkins 1992). It is interesting to see that the blue and red curves are largest for the oldest stars, because this implies one or more of the following: that vertical heating is more effective than planar heating; that older stars were born into vertically hotter orbits than today's young stars; or that mergers were more effective at heating the disk in earlier epochs.

The vertex angle as a function of age has an interesting history. It was noted by Parenago (1958) that the vertex angle of stars changed sharply around a spectral class of F7. Dehnen & Binney (1998) note that this change should occur at the color where the main-sequence lifetime of a star is the age of the Galactic disk, because redward of this break, stars of all ages are being observed. In our data, we are able to look at the vertex angle directly as a function of age, without the need for color as a proxy for age. We find that the younger stars—less than ∼3 Gyr old—have a stronger deviation from symmetry in the velocity plane, and that this deviation drops from 3 Gyr until about 6.5 Gyr of age. Older than 6.5 Gyr, the vertex angle remains constant, but non-zero. This could imply that mixing and blurring play a more important role in the kinematics of stars after 6.5 Gyr of age, while the spiral arms, which heavily influence the deviation of the vertex angle, are more influential when younger than 6.5 Gyr. We will look at this further in Section 6.1.

5.2. Ages and Dispersions

In Figure 3, it can be seen that the very youngest age bin reverses the trend of positive ∂σ_V/∂τ and negative ∂ℓ_V/∂τ. Such a reversal has also been noticed in the bluest stars of Aumer & Binney (2009) and the most metal-rich stars of Anguiano et al. (2018). An uptick in the σV_ϕ distribution for the bluest stars is also seen in Dehnen & Binney (1998), who posit that the very youngest stars are probably not yet well mixed and travel in moving groups on orbits similar to their progenitor clouds.

When we attempt to fit our dispersions to a power law of the form

$$\begin{eqnarray}&&\sigma (x)={\sigma }_{10}{\left(\displaystyle \frac{x+{x}_{1}}{10\mathrm{Gyr}+{x}_{1}}\right)}^{\beta },\end{eqnarray} \tag{ 12 }$$

as in Aumer & Binney (2009), we obtain poor fits because of this uptick.

In Figure 4 we plot the individual and total dispersions using our isochrone grid ages as bins, instead of the large bins used in Figure 3, so that we can more clearly investigate this uptick.

Zoom In Zoom Out Reset image size

We see a jump in dispersion for stars in our old-age edge, which is possibly due to thick-disk or halo contamination. The number of counter-rotating stars (which should most probably be halo stars, and should constitute roughly half of the halo contaminants in the sample) is negligible, about 0.15%; so we expect that this uptick is most likely due to thick-disk stars in our sample falling predominantly into this bin.

For stars younger than 1 Gyr, the velocity dispersions are higher than expected, the dispersions then flatten out until about 2 Gyr, at which point they start to increase with age. This goes against the canonical scenario of stars being predominantly born on cool, circular orbits and then being heated onto more eccentric ones. It could indicate that the younger stars are an unmixed population traveling in moving groups along the paths of their parent clouds (as proposed by Dehnen & Binney 1998).

When we consider only stars that have relative age errors less than 100%, and also implement an iterative, 3σ clipping procedure, we find that a large proportion of the stars younger than ∼1.25 Gyr are removed. The only well populated bin below this cut is at 0.04 Gyr, which is adjacent to the removed edge bin of 0.03 Gyr. Since the majority of the stars contributing to the anomalous behavior in the young age bin are removed by more stringent quality cuts, or abut the edge of our age grid, we treat them skeptically. With this in mind, we fit the velocity dispersions of stars with less than 100% age errors with velocity outliers removed in the age range of 1.25–12.5 Gyr and find

$$\begin{eqnarray}{\beta }_{\mathrm{Tot}},{\sigma }_{\mathrm{Tot},10} & = & (0.36,43.1\,\mathrm{km}\ {{\rm{s}}}^{-1}),\\ {\beta }_{R},{\sigma }_{R,10} & = & (0.33,34.9\,\mathrm{km}\ {{\rm{s}}}^{-1}),\\ {\beta }_{\phi },{\sigma }_{\phi ,10} & = & (0.32,20.2\,\mathrm{km}\ {{\rm{s}}}^{-1}),\\ {\beta }_{Z},{\sigma }_{Z,10} & = & (0.62,16.0\,\mathrm{km}\ {{\rm{s}}}^{-1}).\end{eqnarray} \tag{ 13 }$$

These fits are performed accounting for the observational errors, and should reflect the underlying, intrinsic dispersions of the population. This β_Tot value is similar to the values in Binney et al. (2000) and Aumer & Binney (2009), who found preferred values of 0.33 and 0.35, respectively. The individual components are well within the ranges found in the simulations of Aumer et al. (2016). Holmberg et al. (2009)—who similarly exclude their youngest age bins (which have higher velocity dispersions than the slightly older bins) and oldest age bins—find similar, but slightly higher, values for all components, with the exception of β_Z.

5.3. The Oort Constants

Other commonly used representations of disk dynamics are the so-called Oort constants (Olling & Dehnen 2003; Oort 1927):

$$\begin{eqnarray}2A & = & \displaystyle \frac{\langle {V}_{C}\rangle }{{R}_{0}}-\langle {V}_{C,R}\rangle -\displaystyle \frac{\langle {V}_{R,\phi }\rangle }{{R}_{0}},\\ 2B & =\,- & \displaystyle \frac{\langle {V}_{C}\rangle }{{R}_{0}}-\langle {V}_{C,R}\rangle +\displaystyle \frac{\langle {V}_{R,\phi }\rangle }{{R}_{0}},\\ 2C & =\,- & \displaystyle \frac{\langle {V}_{R}\rangle }{{R}_{0}}+\langle {V}_{R,R}\rangle -\displaystyle \frac{\langle {V}_{C,\phi }\rangle }{{R}_{0}},\\ 2K & = & \displaystyle \frac{\langle {V}_{R}\rangle }{{R}_{0}}+\langle {V}_{R,R}\rangle +\displaystyle \frac{\langle {V}_{C,\phi }\rangle }{{R}_{0}}.\end{eqnarray} \tag{ 14 }$$

The notation $\langle {V}_{R}\rangle$ is used for the mean value of the radial velocity of the population. Since our right-handed coordinate system means that V_ϕ is negative for disk rotation, here we adopt the notation V_C = −V_ϕ for azimuthal velocity. The additional subscripts ,ϕ and ,R indicate derivatives with respect to Galactocentric angle ϕ and radius R. These constants contain information about the motion of the Galaxy, such as the angular rotation speed and the slope of the velocity curve. Since we have spectroscopy (radial velocities) for our data, we may estimate these values directly, but for comparison with the literature, we will also calculate these Oort constants in different age bins.

Figure 5 shows our findings for the calculated Oort constants, as well as associated velocity derivatives ∂V_C/∂R and ∂V_R/∂R, as a function of time. Our average values for young, intermediate-age, and old populations are shown in Table 1.

Zoom In Zoom Out Reset image size

Table 1.  Oort Parameters as a Function of Age (all in units of km s⁻¹ kpc⁻¹).

	Young	Intermediate	Old
A	14.3	18.8	17.0
B	−6.7	−12.7	−2.1
C	0.2	−1.7	−8.3
K	−4.5	−7.6	−5.6
$\tfrac{\partial {V}_{R}}{\partial R}$	−4.4	−9.0	−13.0
$\tfrac{\partial {V}_{C}}{\partial R}$	−7.0	−6.5	−15.8

Download table as:  ASCIITypeset image

Bovy (2017) found (A, B, C, K) = (15.3, −11.9, −3.2, −3.3) km s⁻¹ kpc⁻¹ using main-sequence stars from TGAS; Olling & Dehnen (2003) estimated (A, B, C) = (16, −17, −10) km s⁻¹ kpc⁻¹ using Tycho-2 and ACT data (Urban et al. 1998); Feast & Whitelock (1997) calculated (A, B) = (14.82, −12.37) km s⁻¹ kpc⁻¹ using Cepheids.

Our value for ∂V_R/∂R (which is equivalent to C + K) grows steeper as a function of age and our value for ∂V_C/∂R (which is the same as –A – B) is steepest for the oldest age bin. We find these slopes to be negative for all ages, and most negative for the oldest stars. Our values for ∂V_R/∂R are in agreement with the literature (−6.5 km s⁻¹ kpc⁻¹ according to Bovy 2017 or −4 to −10 km s⁻¹ kpc⁻¹ from Siebert et al. 2011). Our ∂V_C/∂R values are larger than what are generally found (−3.4: Bovy 2017, −1 to −13: Siebert et al. 2011, −2.4: Feast & Whitelock 1997, 1: Olling & Dehnen 2003, −4: Huang et al. 2016; in units of km s⁻¹ kpc⁻¹).

6.1. Age–Metallicity Relation

In Figure 6 we plot the age–metallicity distribution of our sample of stars. The age distribution of our sample peaks at around 3 Gyr, when factoring in the spreads in age probability. Feuillet et al. (2016) found an age distribution peaked at around 3–4 Gyr using giants in APOGEE and Hipparcos; Casagrande et al. (2016) found astroseismic ages for giants observed by Kepler (Borucki et al. 2010) to have a double peak at around 2 and 4 Gyr; Nordström et al. (2004) found ages peaked around 2 Gyr for F and G dwarfs in the smaller volume of the Geneva–Copenhagen Survey sample.

Zoom In Zoom Out Reset image size

To test whether our ages are biased by degeneracy between metallicity and the estimated ages, we can utilize a secondary effect of the "blurring," mechanism which increases the velocity dispersion of similarly aged stars over time. Since the gradient ∂[M/H]/∂R of the ISM is negative, stars born at smaller radii are typically more enriched than stars born at larger radii at the same time. So, when looking at a collection of stars in a small volume, around the same radius, at the same age, the stars rotating more rapidly tend to be more metal-poor, and the stars rotating more slowly tend to be more metal-rich.

Figure 7 examines the angular momentum (L_Z)–velocity plane as a function of both age and enrichment. It can be seen that these two attributes have very different effects on our sample. In terms of age, we can easily see the effects of heating over time: epicyclic heating causes the velocity dispersion to increase, and simultaneously the tilt of the velocity vertex becomes less and less apparent.

Zoom In Zoom Out Reset image size

In terms of metallicity, we can see the effect of ∂[M/H]/∂R of the disk ISM on the velocities of similarly aged stars in our sample. The lower metallicity stars (on orbits with larger guiding center radii) tend to be rotating faster than the higher metallicity stars observed in our volume, which is caused by blurring and the conservation of angular momentum. Hopefully Figure 7 puts to rest some of the worry about possible age–[M/H] degeneracies, because the two spaces are seen to have distinct and differing effects on our sample.

Antoja et al. (2017) noted that, for their thin-disk objects, the more quickly orbiting stars had generally lower metallicities than the stars that were moving with lower V_ϕ speeds, as is expected from blurring. They also see that the stars in their data that were traveling inward (with negative V_R) in the Galaxy had a systemically higher rotational velocity than stars that were traveling outward in the Galaxy. This is a natural consequence of the deviation of the vertex angle. This is all shown nicely in Figure 7 with metal-rich objects moving outward and at lower L_Z, and vice versa for metal-poor objects.

We also see the presence of numerous moving groups in the data. Hercules, for example, while standing out most prominently in the older, higher metallicity populations, is seen faintly in most of the panels. Being present at a wide range of ages and metallicities implies that Hercules is a resonant feature. It is noticeably weaker for younger stars with low metallicity. However, comparing age distributions is difficult due to survey selection effects, such as the volume under consideration or the type of stars being measured (see, for example, Section 3 of Aumer et al. 2016).

6.2. Inside-out Enrichment

If the rotation curve is relatively flat, we may think of L_Z as a proxy for the Galactocentric guiding center radius of the orbit, which should be the birth radius in the absence of churning. In Figure 8 we show the relationship between the angular momentum of stars in our sample and the guiding center radius estimated from the stars' orbits. While not linear, the two are fairly well and monotonically correlated. For the rest of the paper we will use L_Z as a proxy for this guiding center radius, removing any assumptions implied by the potential while retaining the basic information about a star's birth radius.

Zoom In Zoom Out Reset image size

Examining trends in enrichment of populations as a function of age can lead to insights into the formation processes of the Galaxy. In Figure 9 we plot the angular momentum as a function of metallicity and age. We notice that older stars are rotating with less angular momentum than younger stars at the same metallicity. We also see that, for a given age, higher metallicity stars are orbiting with less angular momentum than lower metallicity stars (as seen in, for example, Lee et al. 2011). Since the angular momentum is broadly representative of the birth radius of a star (in the absence of churning), this means that (1) older stars of a given metallicity were born more interior in the Galaxy than younger stars of the same metallicity, and (2) more metal-rich stars of a given age were born more interior to more depleted stars of the same age. These are both expectations of an inside-out formation scenario.

Zoom In Zoom Out Reset image size

In the oldest age bin, at the metal-poor end, there is a downtick to lower angular momenta. This could possibly be explained by contaminating thick-disk stars.

We may also switch the axes to look at metallicity as a function of angular momentum and age (Figure 10). In this figure, using L_Z as a proxy for radius, we see that at all radii, metallicity increases with age, i.e., ∂[M/H]/∂τ is positive. We also see that for any given age ∂[M/H]/∂R is negative. However, these slopes are different for different ages, with there being more of a difference in metallicity between the inner portions of the disk and the outer portions of the disk for old stars; and less of a difference for young stars ($| \partial [{\rm{M}}/{\rm{H}}]/\partial R{| }_{\tau =10}\gt | \partial [{\rm{M}}/{\rm{H}}]/\partial R{| }_{\tau =1}$). Another way to say this is that (∂[M/H]/∂τ)_{outer disk} > (∂[M/H]/∂τ)_{inner disk}, in Figure 10, for example, the inner disk enriches at a rate of about 0.015 dex Gyr⁻¹, while the outer disk enriches at a rate of about 0.027 dex Gyr⁻¹ (between the 2.08 and 10.42 Gyr old populations). This again supports inside-out formation, with the interior regions reaching high levels of enrichment on a smaller timescale than the outer regions.

Zoom In Zoom Out Reset image size

6.3. Migration

In a figure similar to Figure 9, Antoja et al. (2017) found a peculiar behavior where the most metal-rich stars (>+0.5 dex) in their sample were rotating more swiftly than their slightly more metal-poor neighbors (from +0.1 dex to +0.5 dex). A possible interpretation of this is that these stars have had their angular momenta increased by way of churning. This is because, for such metal-rich objects to have traveled to our observational space, they will have needed to sacrifice a lot of angular velocity to retain the same angular momentum, as is the expectation if epicyclic heating (blurring) is the only mechanism at play.

While we see such behavior in our uncut sample, the combination of quality cuts we perform excises these stars from Figure 9 (the most efficacious being the removal of all data that do not fall within 3σ of their errors from an isochrone). As pointed out by R. Schönrich (2018, private communication), these stars could simply be solar-metallicity stars that have been scattered into the high-metallicity tail of the distribution due to large observational errors. In that case our 3σ cut would remove these because their metallicities would not be consistent with our adopted isochrones.

In a separate study of migration from the inner disk, Hayden et al. (2018; see also Kordopatis et al. 2015) calculated orbital parameters for Gaia-ESO (Gilmore et al. 2012) stars and found that a portion of those with iron abundances greater than +0.1 dex had orbital speeds inconsistent with being born in the solar neighborhood or farther interior, as implied by that level of enrichment.⁵ With a similar test, we find that 266 of 1053 (about 25%) of stars in our sample with metallicities over +0.2 dex are rotating faster than the local rotation speed, indicating that they have guiding center radii outside the solar radius. The contradiction between kinematics implying that these stars have guiding center radii external to the solar neighborhood and abundances implying interior birth environments could be rectified by a churning scenario pulling these stars outward.

Interestingly, when we collect a sample of these metal-rich, quickly orbiting stars (top panel of Figure 11), we find their age distribution to be more weighted to younger ages than slowly rotating metal-rich stars. If these are bona fide churned stars, the younger-skewed age distribution could be indicative of churning being a more effective migratory mechanism than blurring for young populations or on short timescales.

Zoom In Zoom Out Reset image size

Also in Figure 11, we plot the age–metallicity relationship for the residual between stars with less angular momentum than the local standard of rest (L_Z ≈ 1952 kpc km s⁻¹) and stars with more angular momentum (the residual between two, two-dimensional histograms, which each integrate to 1). In general, the boundary between the two populations is characterized by a line that gets gradually more metal-rich as time passes, which would represent the enrichment of the solar neighborhood in the absence of churning (with a slope of about ∂[M/H]/∂R ∼ 0.04 dex Gyr⁻¹). The stars that are rotating with less angular momentum than the local standard of rest tend to be more metal-rich than this local "enrichment rate" line since they tend to have been born at smaller radii than the solar circle; and stars that are orbiting with more angular momentum tend to be more metal-poor than this line. This is expected from the conservation of angular momentum as stars are heated onto eccentric orbits.

In this paper we have analyzed data from the TGAS–RAVE–LAMOST cross-matched and co-added data set. With directly measured parallaxes, instead of photometrically or spectroscopically inferred parallaxes, we are able to estimate the ages of these stars with greater precision than previous studies using bulk RAVE and LAMOST data alone.

Using a sample of ∼55,000 high-confidence measurements, which pass a variety of quality cuts, we have calculated the velocity dispersion as a function of age in all three components. Our values are similar to those found by Binney et al. (2000) using Hipparcos data, Aumer & Binney (2009) using Hipparcos data combined with radial velocities from the Geneva–Copenhagen Survey, and Holmberg et al. (2009) again with Geneva–Copenhagen Survey and Hipparcos data (see Equation (13)).

However, to obtain good fits it was necessary to omit the data in our very youngest and very oldest age bins. We expect that the oldest age bin suffers from an amount of thick-disk contamination, so the fact that it has more dispersion than expected is reasonable. Since our observational volume is relatively large, it is perhaps not surprising our sample includes thick-disk stars. Other samples covering similar volumes have also seen an increasing dispersion (e.g., the Hipparcos giant sample of Feuillet et al. 2016), whereas some more local samples have not (e.g., Geneva–Copenhagen Survey study by Holmberg et al. 2009). Another possible explanation for these different findings is that the different studies have different age error magnitudes. Martig et al. (2014) showed that age errors of 30% can effectively smooth and erase jumps in the age–velocity relation, as they see in their simulations with regard to an old and kinematically hot population.

The younger stars, which have an unexpectedly high velocity dispersion, are puzzling, though. We note that Holmberg et al. (2009) also omit their youngest (≲1.5 Gyr) stars from their velocity dispersion fit, and we can see that these stars have similar or even higher velocity dispersions than their slightly older neighbors. Dehnen & Binney (1998) also note that their very bluest stars (indicative of the youngest age group) do not obey the expected velocity dispersion relation. This could be caused by the star-forming regions being on orbits that are not kinematically well-mixed (as noted in Dehnen & Binney 1998). There could also be contamination in these bins from, for instance: blue horizontal branch stars or blue stragglers, which may appear young; or multiple stars, which would have additional components to their velocity dispersions and are more common in high-mass, shorter-lived populations.

When we implement more severe quality cuts on the size of the relative age errors and iteratively remove velocity outliers, we find that the biggest contribution to this rising velocity dispersion is coming from the two youngest age bins, which could indicate that this higher dispersion is just an edge effect.

Another measure of the kinematics, the deviation angle of the velocity vertex, which is formed by local potential non-axisymmetries and is gradually erased for older populations by changing potential (e.g., as the solar neighborhood changes position with respect to the spiral arms over time) and gradual heating, shows a counterintuitive reversal for the very youngest stars as well. Again, this youngest age bin is possibly contaminated by incorrectly aged stars in the edge bins, and not reliable. Curiously, this reversal is also seen by Aumer & Binney (2009) for their very bluest stars, and by Anguiano et al. (2018) for their most metal-rich and least alpha-enhanced stars. In our data, the deviation angle of the velocity vertex stops changing noticeably for ages above about 6.5 Gyr. This could imply a timescale for a local population to become well mixed and no longer show appreciable signs of the potential they were born in.

We have examined the relationship between rotational velocity (as a proxy for orbital radius) and metallicity. There is clear evidence for all radii to be enriching over time, with interior radii always being more enriched than outer radii, lending further support to hypotheses of inside-out formation. However, we also note that the amount of enrichment over time is greater in the outer regions of the disk than in the inner portions, about ∼0.225 dex over the range of our data compared to ∼0.125 dex. This would imply that the central regions formed and enriched quickly, while the outer regions have had a slower enrichment history starting at a later time.

We note that, for our sample of stars, the dividing line between stars rotating more slowly than the local standard of rest (L_Z ≈ 1952 kpc km s⁻¹) and stars rotating more quickly in the age–metallicity diagram—which could be indicative of the enrichment rate of the solar neighborhood if the main method of migration in the Galaxy is blurring—is at ∂[M/H]/∂τ ∼ 0.04 dex Gyr⁻¹. This is fiducially similar to observations by Feuillet et al. (2016) and Bernard et al. (2018), as well as to model parameters in Schönrich & Binney (2009).

We have also examined the slope of the rotation curve and the radial velocity curve. We find them both to be falling as a function of radius, with the oldest stars exhibiting the steepest slopes of descent.

In the data of Antoja et al. (2017), the highest metallicity objects could be seen to be rotating more swiftly than would be expected if they had been transported to the local volume by blurring. We do not see such behavior in our data in Figure 9. However, following the logic of Hayden et al. (2018), we count, in our data, 266 stars with metallicities greater than +0.2 dex that are rotating faster than the local rotation speed. This is out of a total of 1053 stars at this enrichment level or higher—a rate of ∼25%, which is in good agreement with the lower bound on such stars of 20% found in that study.

It is intriguing that the stars that could not have been born in the solar neighborhood (having metallicities greater than +0.2 dex), yet have solar neighborhood dynamics, have a younger age distribution than other stars at the same metallicity, with a surplus at less than 2 Gyr of age). We expect that these stars were moved to orbits with new angular momenta via churning, which is thought to be most effective on cooler, and hence younger, populations. While churning should affect all populations at various points in time, an excess is apparent in the young population. We suppose that blurring, which takes time to build up the radial extent of orbits and tends to flatten metallicity gradients, could spread out this excess on longer timescales.

We thank the developers and maintainers of the following software libraries, which were used in this work: Topcat (Taylor 2005), NumPy (van der Walt et al. 2011), SciPy (Jones et al. 2001), AstroPy (Astropy Collaboration et al. 2013), astroML (VanderPlas et al. 2014), galpy (Bovy 2015), padova,⁶ pytess,⁷ matplotlib (Hunter 2007), IPython (Pérez & Granger 2007), and Python.

We thank Emma Small for invaluable insight, expertise and guidance with regards to isochrones. We also thank Jerry Sellwood, Wyn Evans, Ralph Schönrich, and Victor Debattista for helpful discussions on dynamics.

J.J.V. gratefully acknowledges the support of a LAMOST fellowship, support from the Chinese Academy of Sciences President's International Fellowship Initiative, as well as NSFC grant 11650110439. M.C.S. acknowledges financial support from the CAS One Hundred Talent Fund and from NSFC grants 11673083 and 11333003. This work was also supported by the National Key Basic Research Program of China 2014CB845700.

Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.

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 (SPP 1177 and SFB 881); 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 website is at https://www.rave-survey.org.

This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. M.C.S. is partly supported by the National Key Basic Research and Development Program of China (No. 2018YFA0404501).

The selection function of our data set is non-trivial. The Gaia selection is biased by the sky coverage offered by its orbital position and angle; this manifests as large stripes on the sky where there are too few observations to reliably report data. LAMOST sky coverage is focused in the direction of the anticenter, and certain fields, such as the Kepler field (De Cat et al. 2015), enjoy even more dense coverage. Adding to the confusion is the fact that PIs associated with the project can upload personal target lists that are then factored into the survey. RAVE is more uniform in its sky coverage, and is claimed to be kinematically and chemically unbiased from the expectations (Wojno et al. 2017), but it still has just one hemisphere of coverage and magnitude limits differing from LAMOST and Gaia.

To correct for this complex convolution of selection functions, we must define the goal of our correction: we are interested in a sample that represents the underlying population. Since we are investigating age, metallicity, and velocities, we cannot weight by these factors; instead we will weight by sky coverage (to account for, for example, the heavily sampled Kepler field overpowering the sparsely sampled central regions) and color–magnitude space (to account for, e.g., stellar types that are selected against by parallax requirements, such as giants, or for possible color biases, such as possible searches for blue horizontal branch stars being factored into the LAMOST survey). To calculate our weightings in these spaces, we will compare to synthetic data produced by Galaxia⁸ (Sharma et al. 2011) in accordance with the Besançon Galaxy Model (Robin et al. 2003).

To create an acceptable Galaxia model, we use the default parameter file, and change the r_max value (the maximum observational radius of the model) to 0.55 kpc, and the apparent magnitude limits to 5 and 15. Note that the model is extinction-added to be consistent with the observational space using the internal three-dimensional dust values of the model.

The weighting first accounts for the on-sky density of observations by calculating a Voronoi cell for each star; a Voronoi cell is the area that is closer to the data point defining the cell than to any other data point (see Figure 12). So the inverse of the area of this cell is the density of observations at this individual observation; weighting by the cosine of the latitude accounts for the transformation of the Voronoi cell from cartesian to spherical space. Our observational coverage weight is

$$\begin{eqnarray}&&{W}_{\mathrm{Cell}}=\displaystyle \frac{\cos (b)}{{\rho }_{\mathrm{Cell}}}.\end{eqnarray} \tag{ 15 }$$

Zoom In Zoom Out Reset image size

Next, we account for the expected on-sky density of the observation, normalized (over the model), using our simulated galaxy:

$$\begin{eqnarray}&&{W}_{\mathrm{Sky}}={\rho }_{\mathrm{Model}}(l,b).\end{eqnarray} \tag{ 16 }$$

Finally, we correct for the expected population at a specific (l, b) coordinate by collecting the 1000 nearest points from the model at that (l, b) and construct a 10 × 10 × 10 color–magnitude–distance grid. The relative density in the cell of the observation compared to those of the other cells is the weight of the observation:

$$\begin{eqnarray}&&{W}_{\mathrm{CMD}}={\rho }_{\mathrm{Model}}(J,J-H,l,b,d).\end{eqnarray} \tag{ 17 }$$

For the color–magnitude comparison, we use J magnitude and J − H color from the Two Micron All Sky Survey (Skrutskie et al. 2006) (Figure 13 compares the color–magnitude diagrams of our data and our model). Errors as a function of magnitude are estimated for J, H, and d from the TGAS observational cross-match, and these errors are then added in to the Galaxia mock observations. The weights are then normalized to sum to one before being combined into a total weighting for each star:

$$\begin{eqnarray}&&W={W}_{{\rm{Cell}}}\times {W}_{{\rm{Sky}}}\times {W}_{{\rm{CMD}}}.\end{eqnarray} \tag{ 18 }$$

Zoom In Zoom Out Reset image size

This total weighting is then normalized such that the weights sum to one. This correction does not create an exact duplication of the Galaxia model in any space, because some stars are removed owing to extremely high (or low) weightings, while other stars have measurements that are inconsistent with the model and preclude them from being weighted.

B.1. Metallicity

In the data description, we have referred to the iron abundance [Fe/H] for the LAMOST data, and to the metallicity [M/H] for the RAVE data. The LAMOST pipeline only calculates [Fe/H] and the RAVE pipeline only calculates [M/H] for most stars (and [Fe/H] for a smaller subset). These measures are not equivalent, and if not brought onto similar scales they will produce systematic differences in the estimated ages. To correct this we apply a very basic shift to the LAMOST data such that

$$\begin{eqnarray}[{\rm{M}}/{\rm{H}}]=\left\{\begin{array}{ll}0.86[{\rm{F}}{\rm{e}}/{\rm{H}}]-0.46 & {\rm{f}}{\rm{o}}{\rm{r}}\,[{\rm{F}}{\rm{e}}/{\rm{H}}]\lt -0.7\\ 1.22[{\rm{F}}{\rm{e}}/{\rm{H}}]-0.26 & {\rm{f}}{\rm{o}}{\rm{r}}\,[{\rm{F}}{\rm{e}}/{\rm{H}}]\geqslant -\mathrm{0.7.}\end{array}\right.\end{eqnarray} \tag{ 19 }$$

This shift is obtained by comparing overlapping observations in RAVE and LAMOST (we find about 4100 overlaps that match to within 1''). We then fit a piecewise function to the RAVE [M/H] parameter as a function of the LAMOST [Fe/H] parameter (see Figure 14). This shifts the two surveys onto more equivalent scales, which will hopefully alleviate systematic biases between the two surveys when calculating ages. Rather than using [Fe/H] as a proxy for the full [M/H] content of a star when comparing to isochrones, we are shifting the LAMOST observations, based on the RAVE observations, to a new, hopefully more accurate estimate for the full [M/H] content of these stars. This also allows us to examine the relative age estimates for the RAVE and LAMOST samples together, rather than individually.

Zoom In Zoom Out Reset image size

We have alpha abundance measurements for a subset of our RAVE data from the pipeline of Boeche et al. (2011, see also Kunder et al. 2017), but do not incorporate them into the age estimation. We have investigated the relative effects of alpha abundances on the age estimates, and find them to be similar to the effect of metallicity—increasing the alpha abundances moves the star into a generally cooler, redder parameter space—albeit to a lesser extent. Since we do not have alpha abundances for all of our objects, and since we wish to consider the LAMOST and RAVE data concurrently, we neglect the effects of alpha abundances.

B.2. Velocity

It has been noted by Tian et al. (2015) that LAMOST spectra have a systematic shift of about 5.7 km s⁻¹ when compared to radial velocities from APOGEE (a high-resolution, high S/N spectroscopic survey of red giants). It was also noted by Schönrich & Aumer (2017) that the LAMOST radial velocities appeared to be underestimated by about 5 km s⁻¹. We note a similar offset in our LAMOST–RAVE overlapping data and apply the correction

$$\begin{eqnarray}&&{{\rm{R}}{\rm{V}}}_{{\rm{L}}{\rm{A}}{\rm{M}}{\rm{O}}{\rm{S}}{\rm{T}},{\rm{C}}}={{\rm{R}}{\rm{V}}}_{{\rm{L}}{\rm{A}}{\rm{M}}{\rm{O}}{\rm{S}}{\rm{T}},{\rm{P}}}+5.7\,{{\rm{k}}{\rm{m}}{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 20 }$$

where C and P stand for corrected value and pipeline values, respectively. As noted earlier, we also replace the LAMOST pipeline velocity errors with a uniform error of 7.2 km s⁻¹, which is the standard error found in Schönrich & Aumer (2017).