Yule-Simpson's paradox in Galactic Archaeology

TL;DR

This paper explores how Simpson's paradox affects interpretations of chemo-kinematical relations in Galactic Archaeology, revealing that neglecting stellar age and birth radius can lead to conflicting conclusions about the Milky Way's formation history.

Contribution

It demonstrates the impact of Simpson's paradox on key Galactic relations and emphasizes the need for precise stellar ages to accurately interpret Galactic evolution.

Findings

The MVR in RAVE giants is consistently negative when considering narrow [Mg/Fe] bins.

Radial metallicity gradients cause negative age-group trends, influencing MVR interpretations.

Lithium abundance trends vary with age and birth radius, constraining chemical evolution models.

Abstract

Simpson's paradox, or Yule-Simpson effect, arises when a trend appears in different subsets of data but disappears or reverses when these subsets are combined. We describe here seven cases of this phenomenon for chemo-kinematical relations believed to constrain the Milky Way disk formation and evolution. We show that interpreting trends in relations, such as the radial and vertical chemical abundance gradients, the age-metallicity relation, and the metallicity-rotational velocity relation (MVR), can lead to conflicting conclusions about the Galaxy past if analyses marginalize over stellar age and/or birth radius. It is demonstrated that the MVR in RAVE giants is consistent with being always strongly negative, when narrow bins of [Mg/Fe] are considered. This is directly related to the negative radial metallicity gradients of stars grouped by common age (mono-age populations) due to the inside out disk formation. The effect of the asymmetric drift can then give rise to a positive MVR trend in high-[alpha/Fe] stars, with a slope dependent on a given survey's selection function and observational uncertainties. We also study the variation of lithium abundance, A(Li), with [Fe/H] of AMBRE:HARPS dwarfs. A strong reversal in the positive A(Li)-[Fe/H] trend of the total sample is found for mono-age populations, flattening for younger groups of stars. Dissecting by birth radius shows strengthening in the positive A(Li)-[Fe/H] trend, shifting to higher [Fe/H] with decreasing birth radius; these observational results suggest new constraints on chemical evolution models. This work highlights the necessity for precise age estimates for large stellar samples covering wide spatial regions.

Figures

Full text

††pagerange: Yule-Simpson’s paradox in Galactic Archaeology–LABEL:lastpage††pubyear: 2019††pubyear: 2019††pagerange: Yule-Simpson’s paradox in Galactic Archaeology–LABEL:lastpage

1 Introduction

Galactic Archaeologists use the observed stellar kinematics, chemical abundances, and derived ages as fossil records to recover the chemo-dynamical history of the Milky Way (Freeman & Bland-Hawthorn, 2002; Bland-Hawthorn et al., 2010; Rix & Bovy, 2013; Binney, 2013; Minchev, 2016; Hogg et al., 2016). Age information is by far the hardest one to obtain, yet extremely important to break degeneracies among different chemical evolution models (Miglio et al., 2017).

To infer the Galactic disk history, we rely on observed correlations, such as the mean stellar metallicity111As commonly done in the literature, throughout this work we use ’metallicity’ to mean the iron fraction [Fe/H], not [M/H] - the overall metallicity of a star, defined as the total amount of elements heavier than Helium., [Fe/H], as a function of Galactic disk radius. The negative metallicity gradient observed today suggests that the disk has formed inside out (Matteucci & Francois, 1989). The observational situation is improving rapidly: With the expansion of the disk coverage around the Sun thanks to a number of Galactic surveys, it has become possible in recent years to identify a number of new chemo-kinematical trends.

For example, an inversion of the radial metallicity gradient from negative to positive (and in [$\alpha$/Fe] gradient from positive to negative) for samples at increasingly larger distance from the disk midplane, $|z|$, has been observed in most large Milky Way surveys (e.g., SEGUE - Cheng et al. 2012, RAVE - Boeche et al. 2013b, APOGEE - Anders et al. 2014, Gaia-ESO - Recio-Blanco et al. 2014, LAMOST - Wang et al. 2019), explained by Minchev et al. (2014a) (see also Rahimi et al. 2014; Miranda et al. 2016; Ma et al. 2017; Schönrich & McMillan 2017) as the effect of inside-out disk formation and disk flaring (increasing disk thickness with Galactocentric radius) of mono-age populations (Fig. 10 in Minchev et al. 2014a). It was shown that gradient inversion can result from the variation in the stellar age distribution with radius: at large $|z|$ stars that are old, metal-poor, and [$\alpha$/Fe]-rich dominate in the inner disk (due to the inside-out formation), while younger, more metal-rich, and more [$\alpha$/Fe]-poor stars are more abundant in the outer disk (present at large $|z|$ thanks to disk flaring). This naturally results in a positive metallicity (and a negative [$\alpha$/Fe]) gradient at large $|z|$, although for groups of stars of common age (mono-age populations) the gradient is always negative.

Another example of a trend that can be linked to the Galaxy’s past is disk flaring (e.g., Kazantzidis et al. 2008; Bovy et al. 2016; Kawata et al. 2017; Wang et al. 2018; Bland-Hawthorn & Gerhard 2016). Disk flaring can be related to the dynamical effect of merging satellites, e.g., the Sagittarius dwarf (Laporte et al., 2018; Thomas et al., 2019). Using two suites of simulations of galactic disk formation in the cosmological context, Minchev et al. (2015) showed that flaring is always present for mono-age populations. When considering the total stellar population, however, it was found that the disk flaring was lost or at least significantly decreased. This was explained by a combination of the nested flares of mono-age populations, where younger groups of stars flare and dominate in density at progressively larger radii as a result of the inside-out disk formation.

As a final example, Minchev et al. (2013) (hereafter MCM13) found that in the solar vicinity resulting from their chemo-dynamical model the variation of [Fe/H] with vertical distance from the disk midplane, $|z|$, showed a strong negative gradient ($\sim-0.3$ dex/kpc) but almost flat relations for mono-age populations (see their Fig. 9, top). This phenomenon could be understood by taking into account that the fraction of metal-poor old stars increases with $|z|$, thus giving rise to a negative gradient for the total population (see results with SEGUE G-dwarfs – Schlesinger et al. 2011, RAVE – Boeche et al. 2013b, and LAMOST data – Wang et al. 2019).

What all of the above examples have in common is that a relation resulting from the total sample is weakened, erased, or even inverted, compared to the trends in subsets grouped by considering a third variable (in this case, age).

This phenomenon is known as Simpson’s, or Yule-Simpson’s, paradox (hereafter, YSP), in classical statistics (e.g., Yule 1902; Simpson 1951; Smith & Goltz 2012; Thompson 2006). It has mostly been described in social sciences, medicine, biology, and baseball batting averages, but also in quantum mechanics (Li et al., 2013; Selvitella, 2017). YSP arises when both the independent and dependent variables of a given relation depend on yet a third variable, often referred to as a "lurking" or "confounding" variable (in the above cases - age). A lurking variable is one for which data are unavailable, but it nevertheless has influence on other variables in the study. See, for example, Brase & Brase (2016) for more information on the topic.

Five cases of the YSP, although not identified as such, occurring in the field of Galactic Archaeology have been described in several works in the last five years (Minchev et al., 2013, 2014a, 2014b; Minchev et al., 2015; Minchev et al., 2018). We here elaborate on those and discuss two new cases. We will refer to a "weak" case of Simpson’s paradox when a relation is simply flattened or weakened and a "strong" case when it is completely inverted.

In §§2.2-2.5 we make use of the Milky Way chemo-dynamical model by MCM13, who used a hybrid technique to combine a classical semi-analytical chemical evolution model (Chiappini, 2009) with a simulation in the cosmological context (Martig et al., 2012). The model has been shown to comply with a range of observational constraints, as well as predict new ones (for a summary, see Minchev 2016). In §2.2 and §2.1 we present YSP cases for global Galactic disk chemo-kinematical trends. In §§2.3-2.7 we will be concerned with relations of samples confined to the solar neighborhood.

2 Examples of Simpson’s paradox

2.1 Variation of disk thickness with Galactocentric radius

As a first illustration of YSP, we consider a cosmological zoom-in hydrodynamical simulation, using initial conditions from one of the Aquarius Project haloes (Springel et al., 2008; Scannapieco et al., 2009) and a Tree-PM SPH code. The spatial and mass resolutions are 300 pc and $4.4\times 10^{5}$ M*⊙*, respectively. The stellar mass is $5.5\times 10^{10}M_{\odot}$ and the disk scale-length is $h_{d}=4$ kpc. More details about this simulation can be found in Aumer et al. (2013) (model Aq-D-5).

In the top panel of Fig. 1 color curves show the variation of disk scale-height, $h_{z}$, with galactic radius, $r$, for mono-age populations of width 1 Gyr and median values as seen in the color bar. Single exponential models were fitted to get $h_{z}$ in radial bins of width $0.5h_{d}$ and distance from the disk midplane $h_{d}$ (see details in Minchev et al. 2015). Strong disk flaring (an increase of scale-height with radius) can be seen for the oldest mono-age populations, decreasing for younger subsamples. When the total stellar population is considered, the vertical disk density at all radii requires a sum of two exponentials for a proper fit (Gilmore & Reid, 1983). Intriguingly, the resulting thin (short-dash curve) and thick (long-dash curve) disks do not show much flaring, even though they are composed of the strongly flared mono-age subsamples. This is an example of geometrical (also referred to as structural or morphological) definition of galactic thin and thick disks, as opposed to separation by age or by chemistry in the [$\alpha$/Fe]-[Fe/H] plane for which "low-" and "high-[$\alpha$/Fe] sequences", respectively, are more appropriate names (see discussion in Martig et al. 2016).

This figure highlights a simple example of a "weak" case of YSP, where the strong variation in $h_{z}$ with radius for mono-age subsamples is lost when the total population is considered. As discussed by Minchev et al. (2015), this results because the flaring of younger samples typically dominate at larger radii due to the inside-out disk formation. These "nested flares" constitute a geometric thick disk which does not flare, or in which the flaring is strongly reduced, consistent with the lack of disk flaring found in observations of external galaxies (van der Kruit & Searle, 1982; de Grijs, 1998; Comerón et al., 2011). It can be seen in the bottom panel of Fig. 1 that, indeed, the density of older mono-age groups drops much faster with radius and it flattens for younger populations.

This model of thick disk formation predicts a negative radial age gradient in geometrically-defined thick disks, which was found to be indeed the case for the Milky Way using APOGEE data (Martig et al., 2016). Flaring in Milky Way stellar populations of limited age range has been found by Kalberla et al. (2014) (Red Clump), Feast et al. (2014) (Cepheids), Carraro et al. (2015) (young clusters) and in LAMOST (Wan et al., 2017; Wang et al., 2018; Xiang et al., 2018) and APOGEE (Mackereth et al., 2017) mono-age populations. A proper Milky Way disk mass model is needed to asses the interplay among the flares of different mono-age populations on the overall scale-height variation with radius (that of the total disk mass).

Flaring in mono-age populations from a theoretical point of view appears to be unavoidable, shown to result from merger perturbations (Kazantzidis et al., 2008; Villalobos & Helmi, 2008; Martig et al., 2014), misaligned gas infall (Scannapieco et al., 2009; Roškar et al., 2010), and reorientation of the disk rotation axis (Aumer & White, 2013). A thick disk model incorporating disk flaring is attractive because flaring becomes a necessity rather than a nuisance (as often seen in numerical works in the past).

As will be shown next, inside-out disk formation combined with disk flaring can also explain the inversion of chemical abundance gradients with distance from the disk midplane.

$\bullet$ Possible erroneous interpretation in the absence of age measurements: The Milky Way disk presents no flaring and thus dynamical mechanisms known to flare disks were not at work during its evolution. This is not true if you slice by age - a better interpretation is that mono-age disks always flare but flaring is lost in the total (or wide age range) sample due to the inside-out forming disk. This is a weak YSP case.

2.2 Inversion in radial abundance gradients

As discussed in §1, inversion in the [$\alpha$/Fe] (e.g., [Mg/Fe]) gradient with distance from the disk midplane has been found in a number of Galactic surveys. We investigate the origin of this in the MCM13 model.

The top left panel of Fig. 2 shows the radial [Mg/Fe] gradient for the MCM13 model (similar to Fig. 10 by Minchev et al. 2014a). The black-dashed curve results from the total stellar population with a maximum distances from the disk midplane of $|z|=0.3$ kpc. Different colors correspond to different mono-age populations, as indicated. Although the younger mono-age groups have well-defined positive slopes, that of the total population is strongly flattened inside $r\lesssim 8$ kpc. This can be understood if we consider the relative stellar surface density as a function of radius for each mono-age population, which is indicated by the height of the colored rectangles (conveying the same information as the disk surface density variations with radius shown in the bottom panel of Fig. 1). The old stars – red and orange bins – dominate in the inner disk, causing the mean [Mg/Fe] at $r\lesssim 8$ kpc to go up. On the other hand, the younger populations dominate in the outer disk, as evident from the total mean closely following the $\rm age=5$ Gyr line.

The top right panel of Fig. 2 is similar to the left one, but for a sample at $0.8<|z|<1.5$ kpc. For this significantly higher distance from the disk midplane, the density of old stars is found to strongly dominates at small radii, causing an increase in mean [Mg/Fe] by about $0.2$ dex compared to the sample confined to the disk plane and, thus, a negative gradient.

The bottom row of Fig. 2 presents the effect of uncertainties, as expected from high-resolution spectroscopy (e.g., APOGEE - Majewski et al. 2017, WEAVE - Dalton et al. 2012, 4MOST - de Jong et al. 2012) $\rm\delta[Mg/Fe]=0.05$ dex and asteroseismic ages (e.g., CoRoT, K2 - Howell et al. 2014, TESS - Ricker et al. 2015 and PLATO - Rauer et al. 2014) $\rm\delta age=20\%$, drawn from a Gaussian distribution. Interestingly, the gradient of the total sample is only minimally affected, while a decrease of $\sim 0.05$ dex is found in the oldest mono-age group.

The above illustrates a "strong" YSP case, resulting from the effect of inside-out disk formation combined with disk flaring (see §2.1), giving rise to a different mixture of ages as a function of Galactic radius at different $|z|$ – older, [$\alpha$/Fe]-rich stars are confined to the inner disk and large $|z|$, while younger, [Mg/Fe]-poor stars dominate above the disk plane in the outer disk. Similar explanation holds for the inversion of [Fe/H] gradient from negative close to the disk plane, to positive above it.

$\bullet$ Possible erroneous interpretation in the absence of age measurements: (1) The sign shift in the [$\alpha$/Fe] and [Fe/H] radial gradients found in observations at larger distance form the disk midplane resulted from an outside-in (thick) disk formation or a dominating accreted population at large $|z|$. A better interpretation, according to our model, may be that [Mg/Fe] ([Fe/H]) depends negatively (positively) on $r$ only when one marginalizes over age. This is a strong YSP case. (2) The flattened gradient close to the disk plane is in disagreement with classical chemical evolution models. Steeper gradients are expected in mono-age populations. This is a weak YSP case.

2.3 Birth radius vs vertical velocity dispersion

We now demonstrate that disk flaring and heating as a function of radius can be seen in the local velocity distribution thanks to the effect of stellar radial migration.

In the left panel of Fig. 3 we show the variation of the vertical velocity dispersion, $\sigma_{z}$, with birth radius, $r_{birth}$, for stars found in the simulated solar neighborhood of the MCM13 model. $\sigma_{z}$ is computed as the standard deviation of the vertical velocity in each radial bin. The black-dashed curve shows the total population, giving the impression that the hottest stars in the solar vicinity today have arrived from the inner disk. Looking at this relation it may be tempting to conclude that it is radial migration of hot inner disk stars that give rise to a thick disk. However, dissecting by age (color-coded curves) an opposite relation is revealed for the oldest stellar populations showing that stars born at 2-3 kpc have half the velocity dispersion of coeval populations born at the solar radius. This predicts that the kinematically hottest stars near the Sun today were born locally or in the outer disk - a natural result of the flaring in mono-age disks (see §2.1), which is unavoidable in simulations for different reasons, the most likely possibly being perturbations from orbiting satellites.

The strong positive slope seen in old stars flattens toward redshift zero (present time), becoming slightly negative for age $\lesssim 4$ Gyr, in agreement with the expected effect of migration in quiescent disks (Minchev et al., 2012; Roškar et al., 2013). Note that the flaring induced by outward migration of stars with larger vertical actions and inward migration of outer-disk-born-stars with small vertical actions ($\sim 50\%$ increase in $h_{z}$ in $\sim 4$ disk scale-lengths, see Minchev et al. 2012) is trumped by the impact of mergers and the effect of a cosmological context in general, in which case flaring results from locally born stars, i.e., non-migrators.

The $r_{birth}-\sigma_{z}$ relation was first presented in order to understand why the velocity dispersion increase with [Mg/Fe] in RAVE Giants and SEGUE G-dwarfs showed an abrupt drop at the highest [Mg/Fe] end (Minchev et al., 2014b). It was found that the observations could be explained by kinematically cool, old (thus low metallicity and high [Mg/Fe]) stars arriving to the solar neighborhood from low Galactic radii ($\sim 2-4$ kpc). Indeed, larger changes in guiding radii (or angular momenta) are expected for the coldest subpopulations in a given age group (e.g., Vera-Ciro et al. 2014; Daniel & Wyse 2018).

Minchev et al. (2018) recently presented a largely model-independent method for estimating stellar birth radii based on age and [Fe/H] measurements, using AMBRE:HARPS (de Laverny et al., 2013; De Pascale et al., 2014; Hayden et al., 2017) and HARPS-GTO data (Adibekyan et al., 2012; Delgado Mena et al., 2017) combined with stellar ages computed with the StarHorse code (Santiago et al., 2016; Queiroz et al., 2018), as described in Anders et al. (2018). This opened the possibility to test the $r_{birth}-\sigma_{z}$ relation in observations.

In the right panel of Fig. 3 we show $\sigma_{z}$ vs $r_{birth}$ using HARPS-GTO data for six mono-age populations with median values indicated in the color bar and a bin width $\rm\Delta age=2$ Gyr. Also shown are typical error bars, corresponding to two standard deviations of 1000 realizations in a bootstrapping calculation. As in Minchev et al. (2018), we imposed the following quality criteria on ages and abundances: $\rm\delta[Mg/Fe]<0.07$ dex, $\rm\delta age/age<0.25$ or $\rm\delta age<1$ Gyr, resulting in a sample of 603 stars.

The trends in the data seen in the right panel of Fig. 3 are remarkably similar to the MCM13 model expectation shown in the left panel. The result that stars born at larger radii are kinematically hotter than those born at smaller radii is remarkable as it provides, for the first time, observational evidence that outward migration cools the local disk, rather than heating it, as has been proposed in the past (e.g., Schönrich & Binney 2009b; Roškar et al. 2013). This finding supports work showing that migration does not contribute to thick disk formation (e.g., Minchev et al. 2012; Minchev et al. 2014a; Vera-Ciro et al. 2014; Grand et al. 2016, see also discussion in §1), except in redistributing stars heated by external mechanisms (Quinn et al., 1993; Villalobos & Helmi, 2008) or turbulent gas clouds at high redshift (Bournaud et al., 2009; Forbes et al., 2012).

$\bullet$ Possible erroneous interpretation in the absence of age measurements: The trend in the total population of the $r_{birth}-\sigma_{z}$ relation suggests that migration from the inner disk created the (local) Milky Way thick disk. This is not true if we slice by age, finding exactly the opposite trend for old and intermediate ages in both model and data. A better interpretation is that hotter stars arrived from outer radii due to disk flaring. The outward migrators kinematically cool the local disk. This is a strong YSP case.

2.4 The vertical metallicity gradient

The left panel of Fig. 4 presents the mean [Fe/H] as a function of distance from the disk midplane, $|z|$, for stars in the simulated solar vicinity of the MCM13 model. A negative trend in the total population (black-dashed curve) is seen, consistent with SEGUE G-dwarf (Schlesinger et al., 2011; Rix & Bovy, 2013) and RAVE data (Boeche et al., 2013b). In contrast, mono-age populations (color-coded curves) show much weaker variations - some decrease for old and an increase for the two youngest age bins.

In the right panel of Fig. 4 we show the effect of observational uncertainties, considering $\rm\delta age=20\%$ and $\rm\delta[Fe/H]=0.1$ dex and a Gaussian probability distributions (as in Fig. 2). Similarly to the radial [Mg/Fe] gradients (§2.2), the observational error does not seem to affect much the total sample, however, notable negative trends appear for mono-age populations older than $\sim 4$ Gyr (cyan through red curves). A remarkable drop of $\approx 0.3$ and $\approx 0.4$ dex at high $|z|$ is seen for the green and cyan curves, respectively. Recently, Ciucǎ et al. (2018) found similar negative trends in the older populations in RAVE DR5 (Kunder et al., 2017) data, which, given the results of this section, may be caused by observational uncertainties.

The explanation for the negative vertical metallicity gradient in the total population is again related to the interplay among stellar ages, chemical enrichment, and stellar density, similarly to the radial [Mg/Fe] gradient inversion with $|z|$ (§2.2). Younger metal-rich mono-age groups dominate closer to the disk midplane, while older metal-poor ones are more abundant at larger $|z|$. The fraction of younger to older stars varies with Galactic radius. This both causes a variation in the vertical metallicity gradient (see, e.g., Hayden et al. 2014 using APOGEE data) and an inversion of [Fe/H] and [Mg/Fe] gradients with increasing $|z|$ discussed in §2.2. Depending on the observational biases of different dataset, different slopes in [Fe/H] (or [M/H]) with $|z|$ would be measured in the total population as summarized by Schlesinger et al. (2014), which may result in different interpretations.

Ciucǎ et al. (2018) reasoned that the mild negative gradient seen in the oldest population suggested mono-age disks must have been born flared, which after radial mixing would bring hotter, metal-poor stars from the outer disk to the solar neighborhood, for a given mono-age population. The effect will be the same, however, even if mono-age disks were heated subsequently by external perturbations, i.e., the same scenario producing the nested flares in Fig. 1. The slightly positive gradient in the youngest populations then signifies the opposite effect - migrators from the inner disk weakly heat the local disk, while those arriving from the outer disk cool it, mostly cancelling the overall contribution (as seen in Fig. 5 by Minchev et al. 2012 and in the negative slope of the youngest stars in Fig. 3). The latter effect results only in disk evolution lacking external perturbations.

This section presents a weak YSP case.

2.5 The metallicity-velocity relation (MVR)

A correlation between the observed total metallicity, [M/H], and rotational velocity in the Galactic plane, $V_{\phi}$, of stars in the solar neighborhood was already reported by Carney et al. (1990). Haywood (2008) suggested that this metallicity-velocity relation (MVR), can serve as a test for the amount of stars migrated into the solar neighborhood. He showed that metal-poor solar neighborhood stars at low [$\alpha$/Fe] possess large angular momenta. This meant that they were visitors from the outer disk, currently near the Sun on their pericenters. Schönrich & Binney (2009a) argued that the presence of a negative correlation in the MVR for stars in the chemically defined thin disk (or more correctly, the low-[$\alpha$/Fe] sequence) is an indication of migration. Using a compilation of data for stars in the solar neighborhood, Navarro et al. (2011) found a flat MVR for low-[$\alpha$/Fe] stars, concluding that migration in the solar neighborhood was unimportant. Lee et al. (2011) found a negative relation for low-[$\alpha$/Fe] stars, in contrast to Navarro et al. (2011). A positive slope in the MVR for high-[$\alpha$/Fe] stars was described by Spagna et al. (2010), Lee et al. (2011), Kordopatis et al. (2011), and Adibekyan et al. (2013). More recently, Allende Prieto et al. (2016) used the APOGEE-TGAS sample which utilizes the Gaia DR1 data (Gaia Collaboration et al., 2016), confirmed the positive and negative slopes of the high- and low-[$\alpha$/Fe] stellar populations in the local disk, respectively.

It will be shown below that the positive slope in the MVR of high-[$\alpha$/Fe] stars spanning a wide range in age is just another strong YSP case resulting from the combination of MVRs of mono-age (or mono-[Mg/Fe]222[$\alpha$/Fe] is a good proxy for age for samples confined to a small radial bin (Haywood et al., 2013; Bergemann et al., 2014), but not for data sets spanning large Galactic radii (Minchev et al., 2017).) populations with well-defined negative slopes.333These results were first presented at the RAVE collaboration meeting in Ljubliana, Slovenia in 2014, using RAVE-DR4 (Kordopatis et al., 2013).

We use a subsample from RAVE DR5 (Kunder et al., 2017) with high-quality chemical and kinematical information. We selected a sample of giant stars close to a random magnitude-limited sample as done in Minchev et al. (2014b). We excluded giants with $\rm log(g)<0.5$ to avoid possible effects due to the boundaries of the learning grid used for the automated parameterization, and considered stars in the temperature range $\rm 4000<T_{eff}<5500$ K (thus avoiding horizontal branch stars; see Boeche et al. 2013a). Our final sample consists of stars with signal to noise S/N$>70$ and a distance from the Sun $d<1$ kpc.

The top-left panel of Fig. 5 presents the MVR for our RAVE-TGAS sample. The total population (black-dashed curve) shows a distinct S-shape rotated by $\sim 90^{\circ}$. A positive trend is found in the range ${\rm-1<[Fe/H]<0}$, while a negative one can be seen at ${\rm[Fe/H]>0}$ and ${\rm[Fe/H]<-1}$ dex. When we consider narrow bins in [Mg/Fe] (color-coded curves), however, the MVRs of low-[Mg/Fe] groups show strong negative slopes, flattening and inverting as [Mg/Fe] increases. A V-shape is seen in the highest [Mg/Fe] bin.

The top-middle panel of Fig. 5 is similar to the left one, but showing the MCM13 model, including uncertainties similar to what we expect for the data: $\rm\delta[Mg/Fe]=0.12$ dex, $\rm\delta[Fe/H]=0.1$ dex, and $\rm\delta V_{\phi}=2$ km/s. A remarkable match is seen to the RAVE-TGAS data for both the total population and the subsamples, notably, the gradual inversion from negative to positive with increasing [Mg/Fe] and the V-shape of the highest [Mg/Fe] bin.

Finally, in the top-right panel of Fig. 5 we present the model excluding the uncertainties, which shows a drastic change in shape for individual sub-populations. The MVR trend is always negative, slightly decreasing in steepness for the highest [Mg/Fe] bin. The slope measures $\rm dV_{\phi}/d[Fe/H]=-127$ km/s/dex for stars with median [Mg/Fe] of 0.13 dex (cyan curve), flattening to $\rm dV_{\phi}/d[Fe/H]=-50$ km/s/dex for the highest [Mg/Fe] bin with median value of 0.34 dex (maroon curve).

2.5.1 Effect of observational uncertainties

To understand the effect of observational uncertainties seen in the model, in the bottom row of Fig. 5 we show the [Fe/H] distributions of the total samples and the corresponding mono-[Mg/Fe] populations. In both data and model we see similar shifts in the peaks to lower [Fe/H] with increasing [Mg/Fe], as expected.

The small pink circles in the top-right panel of Fig. 5 indicate the approximate [Fe/H] values at which the density peaks for each mono-[Mg/Fe] population. We see that those act as pivot points around which the MVR flattens when neighboring bins are allowed to communicate, i.e., when uncertainties are convolved with the model (middle column). The only exception is the highest-[Mg/Fe] bin, because it is significantly lower in density than the intermediate bins, and it lacks a clear peak.

In a range of [Fe/H] over which the densities of two adjacent mono-[Mg/Fe] populations are comparable, both MVRs are deflected similarly. For example, the two highest-[Mg/Fe] bins in the range $-1.0<$[Fe/H]$<-0.5$ dex show similar densities, consequently, the maroon curve shifts upwards and red curve downwards as error is implemented. More generally, the MVRs of density tails are affected by the neighboring bin’s peaks in the corresponding direction: since the slope is always negative, metal-poor tails are shifted toward larger $V_{\phi}$ and vice versa. It should be noted that relatively similar densities in the observations are not to be taken at face value due to the selection function of a particular survey. For example, the RAVE survey is known to favor low-[$\alpha$/Fe] stars (Boeche et al., 2013b), which is evident in the distributions of mono-[Mg/Fe] populations in the bottom left panel of Fig. 5.

The negative MVR slope at [Fe/H]$<-0.9$ dex is minimally affected, as these are the stars with the highest [Mg/Fe] values and the uncertainties are not large enough to allow other [Mg/Fe] bins to reach that end of the metallicity distribution. This distinct shape is present in both data and model, as well as seen in other works, although not acknowledged. For example, such an inversion can be found (1) in the bottom panel of Fig. 7 by Lee et al. (2011), if one focussed on the mean values, rather than the line fit, and (2) in the top panel of Fig. 4 by Adibekyan et al. (2013) at $\rm[Fe/H]\approx-1.1$ dex.

Similarly to the effect of errors, using wider [Mg/Fe] ranges would result in overlapping between the neighboring bins (all of which show clearly negative MVRs), approaching the black-dashed line in the right panel of Fig. 5. If we split into the high- and low-[Mg/Fe] sequences, as often done in the literature, we would end up with a positive and a negative MVR, respectively.

In summary, the genuinely negative MVR can become positive if (1) a sample spans a large range of [Mg/Fe] or (2) uncertainties are not sufficiently small.

2.5.2 What causes the invariantly negative MVR?

The negative trends for all mono-[Mg/Fe] populations (right panel of Fig. 5) are easy to understand. While a star at the solar radius, $r_{0}$, moving on a perfectly circular orbit rotates with the value of the circular velocity (here $V_{0}=220$ km/s), stars with guiding radii inside/outside $r_{0}$ will rotate slower/faster than that. Due to the exponential stellar surface density decrease with radius, however, more low-angular momentum stars are present in the solar neighborhood (this is the asymmetric drift effect, see Binney & Tremaine 2008). Consequently, the mean tangential velocity, $V_{\phi}$, of the total population is always below $V_{0}$. Since our sample is selected from a narrow radial range, the angular momentum is given by $L\approx V_{\phi}r_{0}$ and in terms of the stellar guiding radius, $r_{g}$, as $L\approx V_{0}r_{g}$. Therefore, $V_{\phi}\approx(V_{0}/r_{0})r_{g}$, i.e., $V_{\phi}$ can be thought of as a proxy for $r_{g}$. We will not find the youngest stars with guiding radii far away from $r_{0}$, as those will be on orbits too close to circular to be able to reach our sample. The total population shows a peak at solar [Fe/H] because these are the youngest stars near the solar radius, thus having the most circular orbits.

Splitting the model sample in narrow bins of [Mg/Fe] in the top-right panel of Fig. 5 showed that the MVR is always negative. This simply reflects the radial negative metallicity gradient of a mono-age, or in this case, a mono-[Mg/Fe] population. Since $V_{\phi}\sim r_{g}$, for a given [Mg/Fe] bin, a negative metallicity gradient is reflected into a negative MVR gradient.

In summary, the positive MVR gradient in the high-[$\alpha$/Fe] stars measured in a number of works (Carney et al., 1990; Rocha-Pinto et al., 2006; Spagna et al., 2010; Lee et al., 2011; Adibekyan et al., 2013; Allende Prieto et al., 2016) is simply the result of the asymmetric drift and the negative radial metallicity gradient of mono-age populations.

$\bullet$ Possible erroneous interpretation in the absence of age measurements: The MVR slope of a mixed-age population tells us something about the Milky Way migration history. In fact the MVR slope is strongly dependent on the sample selection function and the age range of the sample. When dissected by age (or [$\alpha$/Fe]), one finds that all groups have negative slopes closely related to the radial metallicity gradients of mono-age populations. This is a strong YSP case.

2.6 The age-metallicity relation (AMR)

The age-metallicity relation (AMR) in the solar vicinity is believed to give us a measure of how fast the iron content of stars has increased with cosmic time. It has been found somewhat confusing that the AMR is quite flat, or even reversed, for stars with $\rm age\lesssim 10$ Gyr (e.g., Edvardsson et al. 1993; Lin et al. 2018; Vickers & Smith 2018). This has been taken sometimes as evidence that stars at the solar radius have not enriched in iron during that period, which is in disagreement with predictions from Milky Way chemical evolution models, or thought to result from the target selection of a given sample. While it is widely accepted now that radial migration is responsible for the scatter in [Fe/H] at the same age (Sellwood & Binney, 2002; Roškar et al., 2008; Schönrich & Binney, 2009a), it has only recently been suggested that stars migrating from different radii present well-defined AMRs but conspire to flatten the relation resulting from the total population (Minchev et al., 2018).

In Fig. 6 we plot the AMR for the local HARPS-GTO data, which was used in the right panel of Fig. 3. The black-dashed curve shows the relation for the total sample. The color-coded lines represent five mono-$r_{birth}$ populations, as color-coded, of width 2 kpc and median values given in the color bar. Birth radii were estimated from the stellar age and [Fe/H], as described by Minchev et al. (2018). This is not an unbiased result since the $r_{birth}$ estimate is based on age and [Fe/H]. Nevertheless, this figure reflects how the functional form of the ISM metallicity evolution with radius and cosmic time (also estimated by Minchev et al. 2018) determines the shape of the mono-$r_{birth}$ AMRs.

In Fig. 6 we measure a gradient of 0.02 dex/Gyr for $\rm age<8.5$ Gyr in the total population (black-dashed curve). However, the slope in the AMR of stars born in the range $6<r_{birth}<8$ kpc (blue curve) is 0.04 dex/Gyr and is similar for the other $r_{birth}$ bins. This figure provides an explanation for the flatness of the local AMR, when radial migration is not taken into account. The slope of the total population will depend crucially on selection effects of a particular survey and the quality cuts applied. In this particular case we find quite a steep one compared to other data sets. For example, from Fig. 7 by Lin et al. (2018) we estimated $\sim 0.013$ dex/Gyr for $\rm age<10$ Gyr, while Vickers & Smith (2018) find a double peak at $\sim 3$ and $\sim 8$ Gyr. While the shape indeed depends on the particular target selection when the total sample is considered Casagrande et al. (2016), such effects will be significantly reduced if we split in mono-$r_{birth}$ populations, relying on the assumption that stars with a particular value of [Fe/H] and age were born at roughly the same radius. Therefore, we do not expect the AMRs of mono-$r_{birth}$ populations to deviate much in other data sets from what we show here.

In summary, the chemical enrichment and the radial migration history are so closely entwined in the data that they can only be inferred together, as showed by Minchev et al. (2018) and Frankel et al. (2018).

$\bullet$ Possible erroneous interpretation in the absence $r_{birth}$ measurements: The metallicity at the solar radius has not increased in the last 8-10 Gyr of evolution. This is not true if mono-$r_{birth}$ populations are considered, showing well defined slopes in the AMR of stars born at a given radial bin. This is a weak YSP case.

2.7 Lithium abundance variation with [Fe/H]

The left panel of Fig. 7 shows the abundance of lithium, A(Li), vs [Fe/H] for our local HARPS data used in §2.3 and §2.6. The difference here is that the lithium is estimated from the AMBRE:HARPS sample (as detailed by Guiglion et al. 2016a, b) and crossmatched with the ages derived from HARPS-GTO stellar parameters (Anders et al. 2018, using the StarHorse code), resulting in 326 dwarf stars. The black-dashed curve shows the mean resulting from the total population and the color-coded lines represent six mono-age populations, as in Fig. 3. While the relation from the total sample has a positive slope, mono-age populations show negative trends, except for the two youngest age bins. The variation of mean A(Li) with [Fe/H] for mono-age groups, to our knowledge, has not been shown before in observations, although it could be argued would result from the scatter plots shown in Figures 5 and 7 by Ramírez et al. (2012) and should correspond to bins of common stellar mass, shown in several works (e.g., Nissen & Schuster 2012; Ramírez et al. 2012; Delgado Mena et al. 2015; Bensby & Lind 2018). The flattening we find in A(Li)-[Fe/H] for younger groups of stars is very similar to the flattening seen with increasing stellar mass in Fig. 6 by Bensby & Lind (2018), as expected from the correspondence between stellar mass and age. This is, once again, a clear strong case of YSP, akin to the $r_{birth}-\sigma_{z}$ relation of §2.3 and the MVR of §2.5, which all show a dramatic reversal of the trend in the total local population.

Studies are usually concerned with the upper envelope (not shown here) in the A(Li)-[Fe/H] relation of dwarf stars, reasoning that this is closer to the original lithium abundance in the ISM, since depletion over the lifetime of stars takes place (e.g., Rebolo et al. 1988; Romano et al. 1999; Lambert & Reddy 2004; Guiglion et al. 2016b; Cescutti & Molaro 2019; Guiglion et al. 2019). This upper envelope, which will be traced by extending the metal-poor tails of the mono-age populations in Fig. 7, is found to delineate stars from old to young when moving from low to high [Fe/H]. Further work is needed to understand if this is simply due to the inability of young stars to reach the solar neighborhood on apo- and peri-centers as their orbits are close to circular, and the insufficient time for migration.

Splitting the same data by birth radius (right panel of Fig. 7), we find that mono-$r_{birth}$ populations exhibit mostly positive trends as in the total sample, but with significantly steeper slopes (a weak YSP case, similar to the AMR relation, see Fig. 6). This general shift from high to low [Fe/H] for inner to outer $r_{birth}$ bins is also seen in the chemical evolution models by Prantzos et al. (2017) and Chiappini (2009) (shown by Guiglion et al. 2019). In those models, however, the variation of A(Li) with [Fe/H] is shown for the gas, i.e., at the time of stellar birth. Since lithium is depleted throughout stars’ lifetime, we need to consider the upper envelope for each $r_{birth}$ bin. The above models provide a good match to the trends seen in the right panel of Fig. 7 assuming the scatter around the mean for all mono-$r_{birth}$ groups is similar.

Delgado Mena et al. (2015) described a lithium decrease at $\rm[Fe/H]>0$ using HARPS data, a result subsequently confirmed by Guiglion et al. (2016b) using a larger dataset. Recently Guiglion et al. (2019) reasoned that the observed decline in lithium abundance at super-solar metallicity can be explained as the effect radial migration, as stars in the metal-rich tail are expected to have been born roughly in the range $2<r<5$ kpc (see Fig. 8 by Minchev et al. 2018). The right panel of Fig. 7 shows this is indeed the case, using our estimated birth radii for HARPS data.

It should be noted that our results are not dependent on any chemical evolution or dynamical modeling, except for a simple assumption on keeping the mono-$r_{birth}$ distributions physically meaningful (see for details Minchev et al. 2018). It is clear that, depending on the birth radius bin, the upper envelope will have different shape, notably, shifting to higher [Fe/H] and having a faster drop at the [Fe/H]-poor tail, as $r_{birth}$ decreases. This opens up the possibility to study the production/destruction of lithium as a function of birth radius, that will provide stronger constraints on chemical evolution models.

$\bullet$ Possible erroneous interpretation in the absence of age and $r_{birth}$ measurements: In the total sample the upper envelope in the A(Li)-[Fe/H] plane will be biased toward older stars at the higher and (especially) lower [Fe/H] tails (see Fig. 2 by Minchev et al. 2018), in which case the inferred ISM evolution of lithium will be incorrect. This problem can be resolved with the help of age and $r_{birth}$ estimates. Dissecting the A(Li)-[Fe/H] relation by age presents a strong YSP case and that by $r_{birth}$ - a weak one.

3 Conclusions

We showed in this work that Yule-Simpson’s paradox (YSP) is omnipresent in the field of Galactic Archaeology. Chemo-kinematical relations believed to constrain the Milky Way past can look completely different, depending on whether a subset of stars with common age or a mixture of ages is considered (and similarly for birth radius). In other words, a statistical relationship between two variables does not necessarily represent a cause-and-effect relationship. We can summarize our results as follows:

$\bullet$ We described seven YSP cases, which we classified as "weak" if a relation was flattened or weakened, or "strong" if it was completely reversed. Weak cases were found in the variation of galactic disk thickness with Galactocentric radius (§2.1), the vertical metallicity gradient (§2.4), the age-metallicity relation (AMR, §2.6), and the A(Li)-[Fe/H] relation when split by birth radius (§2.7). Strong cases appeared in the inversion of radial abundance gradients with distance from the disk midplane (§2.2), the $r_{birth}-\sigma_{z}$ relation (§2.3), the metallicity-rotational velocity relation (MVR, §2.5), and the A(Li)-[Fe/H] relation when split by age (§2.7).

$\bullet$ YSP is found in both global (§2.1 and §2.2) and local (samples confined to the solar neighborhood, §§2.3-2.7) chemo-kinematical trends.

$\bullet$ The shape of the MVR of stars in the solar vicinity has been suggested as a discriminant for the amount of radial migration suffered by the local disk, where a flat (inverse) relation has been argued to signify the unimportance (importance) of migration. Using RAVE data and the MCM13 chemodynamical model, we showed that the slope in the MVR for any narrow [Mg/Fe] subpopulation (used as a proxy for age) is always negative. An inversion from a negative to a positive slope of high-[$\alpha$/Fe] stars can be caused by the increasingly lower angular momentum of older, kinematically hotter stars at [Fe/H]$\lesssim-0.5$ dex with a slope measurement dependent on the selection function of a particular survey, as well as the observational uncertainties. The persistent MVR negative slope for any [$\alpha$/Fe] bin simply reflects the Milky Way negative metallicity gradient of mono-age populations due to the inside-out disk formation, combined with the effect of the asymmetric drift.

$\bullet$ Using estimates of stellar birth radii for AMBRE:HARPS data, we showed that the A(Li)-[Fe/H] relation in the solar neighborhood looks different for different mono-$r_{birth}$ populations. This can provide new constraints on chemical evolution models by requiring to match the upper envelope of the relation in different $r_{birth}$ bins, as well as their combination at a given final radius expected from the effect of radial migration.

Depending on the observational biases of different dataset, different trends can be measured in the total population of a given relation, thus resulting in various interpretations. Dissecting by age, however, decreases strongly such biases, as coeval stars will be affected by dynamical processes similarly, for a given radius. Moreover, one need not worry what age a particular stellar sample is biased toward, as long as the relation of interest does not involve the stellar density. For example, in the case of the vertical metallicity gradient, assuming stars with the same metallicity and age are born at the same radius (i.e., no significant intrinsic abundance scatter at the time of star formation, cf, Spitoni et al. 2018), the only effect we need to worry about is radial migration.

We demonstrated that splitting a stellar sample by age or birth radius is of utmost importance for recovering the Galaxy evolution. Given that our birth radius estimate depends crucially on precise ages measurements, we conclude that ages are urgently needed to make progress in the field of Galactic Archaeology. We expect to get those form the K2 (Howell et al., 2014), TESS (Ricker et al., 2015) and PLATO (Rauer et al., 2014) asteroseismic missions in the near future, as well as isochrone age determination from spectroscopic surveys, such as GALAH (De Silva et al., 2015), Gaia-ESO (Gilmore et al., 2012), APOGEE, WEAVE, SDSS-V’s Milky Way Mapper (Kollmeier et al., 2017), and 4MOST, combined with Gaia astrometry.

Simpson’s paradox undoubtedly exists in other Galactic Archaeology relations, as well as other areas of astrophysics. There may be other lurking variables in addition to age and birth radius we discussed here, which must be accounted for in statistical experiments, in order to avoid non-meaningful results. We hope this work brings awareness to different communities and helps advance our understanding of the Milky Way formation and evolution.