The chemodynamics of prograde and retrograde Milky Way stars

Context: The accretion history of the Milky Way is still unknown, despite the recent discovery of stellar systems that stand out in terms of their energy-angular momentum space, such as Gaia-Enceladus-Sausage. In particular, it is still unclear how these groups are linked and to what extent they are well-mixed. Aims: We investigate the similarities and differences in the properties between the prograde and retrograde (counter-rotating) stars and set those results in context by using the properties of Gaia-Enceladus-Sausage, Thamnos/Sequoia, and other suggested accreted populations. Methods: We used the stellar metallicities of the major large spectroscopic surveys (APOGEE, Gaia-ESO, GALAH, LAMOST, RAVE, SEGUE) in combination with astrometric and photometric data from Gaia's second data-release. We investigated the presence of radial and vertical metallicity gradients as well as the possible correlations between the azimuthal velocity, $v_\phi,$ and metallicity, [M/H], as qualitative indicators of the presence of mixed populations. Results: We find that a handful of super metal-rich stars exist on retrograde orbits at various distances from the Galactic center and the Galactic plane. We also find that the counter-rotating stars appear to be a well-mixed population, exhibiting radial and vertical metallicity gradients on the order of $\sim$ -0.04 dex/kpc and -0.06 dex/kpc, respectively, with little (if any) variation when different regions of the Galaxy are probed. The prograde stars show a $v_\phi$-[M/H] relation that flattens -- and, perhaps, even reverses as a function of distance from the plane. Retrograde samples selected to roughly probe Thamnos and Gaia-Enceladus-Sausage appear to be different populations yet they also appear to be quite linked, as they follow the same trend in terms of the eccentricity versus metallicity space.


Introduction
The role that accretion events have played in the evolution of the Milky Way and the quest to find possible remnants that are at the origin of the old disc have been central topics in Galactic archaeology for over than five decades (e.g. Eggen et al. 1962;Searle & Zinn 1978;Gilmore & Reid 1983;Chiba & Beers 2000;Gilmore et al. 2002;Wyse et al. 2006;Kordopatis et al. 2011). The advent of the Tycho-Gaia Astrometric Solution (TGAS) catalogue of the European space mission, Gaia (Gaia Collaboration et al. 2016b,a), and, in particular, the second Gaia data release (GDR2, Gaia Collaboration et al. 2018a), have enabled us to measure with much greater accuracy the positions and the 3D velocities of millions of stars in a volume several kiloparsecs wide. Such studies have shed an unprecedented light on the questions cited above.
On the one hand, the discovery of ripples in the Galactic disc has shown that its morphology and characteristics, at all radii, continue to be impacted by external factors (Gaia Collaboration et al. 2018c;Antoja et al. 2018;Laporte et al. 2018Laporte et al. , 2019. On the other hand, the discovery of kinematic groups outside the disc, such as the Gaia-Sausage Myeong et al. 2018b), Gaia-Enceladus , Sequoia (Myeong et al. 2019), and Thamnos (Koppelman et al. 2019) are all believed to be remnants of past accretions; however, it is still unclear whether or not they are distinct features and to what extent they might have contributed at the formation of the thick disc or the inner halo (e.g. Haywood et al. 2018;Fernández-Alvar et al. 2019;Belokurov et al. 2020).
In that respect, retrograde stars in the Milky Way hold a key place in our understanding of the assembly history of our Galaxy because there is no clear mechanism that could form them exclusively in situ. Counter-rotating stars in the halo have been discussed since Majewski (1992); Carney et al. (1996); Carollo et al. (2007); Nissen & Schuster (2010); Majewski et al. (2012). In particular, Majewski (1992) identified, via an investigation of proper motions and a multi-colour analysis of a sample of few hundred stars towards the north Galactic pole, a retrograde rotation among stars reaching 5 kpc from the plane. Majewski (1992) reports that they exhibit no radial metallicity gradient, while also stating that this population may be younger, on average, than the dynamically hot metal-poor stars that are closer to the plane. This result was later confirmed by Carney et al. (1996) using a kinematically biased sample of 1500 stars and by Carollo et al. (2007) using calibration data from 20 000 stars from the Sloan Extension for Galactic Understanding and Exploration (SEGUE, Yanny et al. 2009). Nissen & Schuster (2010) carried out a spectroscopic observation at high resolution of a sample of 94 kinematically selected dwarf stars and found that the low-[α/Fe] sequence stars identified in the [α/Fe] vs [Fe/H] plane were mostly counterrotating targets. According to the authors, some of the low-[α/Fe] sequence stars may have originated from ω Centauri's Article number, page 1 of 16 arXiv:2009.07192v2 [astro-ph.GA] 22 Sep 2020 globular cluster progenitor (the latter also being on a retrograde orbit; see e.g. Dinescu et al. 2002), a hypothesis that has also been suggested by Majewski et al. (2012), using low-resolution (R ∼ 2600) spectra of ∼ 3 000 stars, and by Myeong et al. (2018a), using a catalogue of ∼ 62 000 halo stars from a crossmatch of Gaia DR1, the SDSS DR 9 (Ahn et al. 2012), and LAM-OST DR 3 (Luo et al. 2015) catalogues.
Using the exquisite GDR2 data, Gaia Collaboration et al. (2018b) found that a non-negligible amount of the halo stars in the extended Solar neighbourhood that are somewhat less bounded than the Sun are on retrograde orbits. Compared with cosmological simulations, they conclude that this rate of counter-rotating stars, despite being rare, is not unexpected. Helmi et al. (2018) suggested that they originate from the merger with Gaia-Enceladus, the latter in a slightly counter-rotating orbit with a mass ratio of 4:1 (see also Gallart et al. 2019;Feuillet et al. 2020).
While the metallicity extent, chemical structure as well as the number of subpopulations that constitute these retrograde stars have already been investigated in previous papers in actionangular momentum space (e.g. Helmi et al. 2018;Myeong et al. 2018bMyeong et al. , 2019Koppelman et al. 2019;Feuillet et al. 2020;Naidu et al. 2020), the investigation of the correlation between the stellar azimuthal velocity, v φ , and the stellar metallicity, [M/H], or iron abundance, [Fe/H], can provide an additional insight in this endeavour (e.g. Kordopatis et al. 2013a, for the identification of thick disc stars and characterisation of its properties). The level of correlation between two or more parameters can also highlight the mechanisms that have formed and shaped a stellar population that is being considered.
While this was already visible in previous datasets (see e.g. Carney et al. 1996), Spagna et al. (2010); Kordopatis et al. (2011);Lee et al. (2011) isolated the correlation between kinematics and metallicity for thick disc stars and measured it to be on the order of ∂v φ /∂[M/H] ≈ +50 km s −1 dex −1 . A mild negative correlation has also been measured for the chemically identified thin disc stars in Recio-Blanco et al. (2014); Kordopatis et al. (2017); Allende Prieto et al. (2016). Whereas for the thin disc, this correlation is well-understood as the effect of blurring and churning in the disc (e.g. Sellwood & Binney 2002;Schönrich & Binney 2009;Minchev & Famaey 2010), multiple explanations can be found in the literature concerning the origin of the positive correlation for the thick disc stars. Indeed, it has been suggested that it could either be the signature of the collapse of a primitive gas cloud (e.g. Kordopatis et al. 2017), a signature of inside-out formation and gas re-distribution in the primitive disc (Schönrich & McMillan 2017), or, as recently suggested by Minchev et al. (2019), a correlation resulting from the superposition of mono-[α/Fe] sub-populations with negative slopes (as in the thin disc). In the latter case, the measured positive correlation is due to the combination of several populations of different ages with different relative weights and proportions (the so-called Yule-Simpson paradox).
In this paper, we aim to study the chemokinematics of the stars in the Gaia-sphere, with a particular focus on the retrograde targets, in order to identify trends that could shed some light on the formation origins of the thick disc and the origin of those retrograde stars. Section 2 describes the dataset used in this analysis. Sections 3 and 4 show the spatial metallicity gradients and the correlations between the azimuthal velocity and metallicity for the prograde and retrograde stars, whereas Sect. 5 discusses selections probing Gaia-Enceladus-Sausage stars and other accreted populations, selected in the phase-space. Finally, Sect. 6 presents our conclusions.
We applied the quality-flag present in the catalogue to remove all the stars for which the spectroscopic parameters are too far from the isochrones to have reliable distances and ages (i.e. those that do not have the value of flag equal to zero; see Sanders & Das 2018, for more details). We also removed the RAVE-on (Casey et al. 2017) entries in order to avoid duplicates with RAVE-DR5. As far as the T eff range is concerned, we only kept the stars between 3500 K and 6800 K in order to avoid too cool or too hot stars for which spectra parameterisation is intrinsically difficult to obtain and its results uncertain, thus introducing a potential bias 1 .
In addition to the filters above, further cuts were required in order to make sure our sample contained only single stars. To achieve this, we cross-matched the sample of stars that fulfil the above criteria with the GDR2 archive, based on the GDR2 sourceid. We extracted the renormalised unit weight error (RUWE) of the targets and discarded the stars with a RUWE greater than 1.2 (suggesting that Gaia's astrometric solution has not converged appropriately and that the considered stars are potential binaries) 2 , an astrometric excess noise greater than 1, and, finally, a parallax relative uncertainty (σ / ) greater than 0.1. The latter filter ensures that what dominates the stellar distance estimation is the parallax measurement and not the prior adopted in Sanders & Das (2018). We note that the adopted distances do not take into consideration the zero-point offset reported for GDR2 parallaxes (e.g. Lindegren et al. 2018;Schönrich et al. 2019), as it has a dependence on the position on the sky, the magnitude, and the colour of the star, which complicates an application that does not involve the introduction of further biases (e.g. Gaia Collaboration et al. 2018a; Arenou et al. 2018). That said, the global effect of this zero point on the velocities and gradients is discussed in the relevant sections below and in greater detail in Appendix A.
We also discarded the stars that have an uncertainty in metallicity greater than 0.2 dex, a v φ uncertainty greater than 50 km s −1 , and a Galactocentric cartesian X position that suggests they are located past the Galactic center (in order to avoid probing entirely different regions of the Galaxy at a given R). Our final working sample was obtained by removing the intersurvey repeats, using the duplicate keyword in the Sanders & Das (2018) table, which preferentially selects the stars in the order: APOGEE, GALAH, GES, RAVE, LAMOST, and SEGUE (for a comparison between the metallicity, the distance, and the azimuthal velocity for the repeated inter-survey stars, see Appendix B). Eventually, we ended up with 2 419 655 unique stars, out of the 4 906 746 entries in the initial catalogue. The relative fractions of each catalogue as a function of the Galactocentric radius, R, and absolute distance from the Galactic plane, |Z|, are shown in Fig. 1 Fig. 1. Relative fraction of targets within each survey for the stars closer than 5 kpc from the Galactic plane after application of the quality cuts described in Sect. 2, as a function of the Galactocentric radius, R (top), and as a function of absolute distance from the Galactic plane, |Z| (bottom). uncertainties in distance, metallicity and v φ , split by individual surveys, are shown in Fig. 2. As anticipated, this figure shows that low-resolution surveys tend to have larger uncertainties in metallicity, and, to some extent, v φ , yet the uncertainties in distance are mostly driven by the apparent magnitude of the targets.
The stellar orbits were computed using the galpy code (Bovy 2014) with the MWPotential2014 and the actionangle formalism for axisymmetric potentials using Binney (2012)'s Staeckel approximation. To be compatible with Sanders & Das (2018) velocities and priors, we adopted the Solar peculiar velocity from Schönrich et al. (2010), the velocity of the local standard of rest as V LSR = 240 km s −1 , and the Sun's position equal to (R, Z) = (8.2, 0.015) kpc.
Once we had our final dataset in hand, we verified again that the metallicities of the stars belonging to the different surveys are roughly on the same scale. This was achieved by visually inspecting that the metallicity distributions close to the Sun's position (|R − R| < 0.2 kpc, |Z| < 0.2 kpc) peaked at the same value. The metallicity distributions, shown in Fig. 3, exhibit a good agreement given the distinct selection functions and the discrepancies already identified in Fig. B.1.  Figure 4 shows the Kiel diagram of the prograde (v φ > 0 km s −1 , 2 397 183 stars) and retrograde (v φ < 0 km s −1 , 22 472 stars) samples. Simply by comparing the two diagrams, especially at the turn-off and red giant branch (RGB) regions, we can already notice that the age range of the retrograde stars is smaller than that of the prograde stars, yet the large width of the turn-off as well as the width of the RGB, suggests that retrograde stars encompass a range of ages. Figures 5 and 6 show the radial and vertical gradients for the two populations for selected ranges of distances from the plane and Galactocentric radii. The associated measured gradients are reported in Tables 1 and 2, respectively, and are discussed in the following subsections.

Metallicity gradients for the prograde stars
We find that the radial metallicity gradient of the prograde stars exhibits two regimes. First, in the inner disc (5 ≤ R ≤ 8 kpc), a globally null gradient is present, at all |Z| (0.01 ± 0.002 dex kpc −1 closest to the plane 3 , −0.005 ± 0.005 dex kpc −1 in the Z = [1 − 2] kpc bin).
For Galactocentric distances greater than 8 kpc, the metallicity gradient flattens as one moves further from the plane, regardless of whether we correct for the zero-point offset on the parallaxes (see also Tables A.1 and A.2). For the distances derived without the zero-point correction, the metallicity gradient eventually reverts from a negative one to a slightly positive one at distances where the thick disc dominates (i.e. |Z| > 2 kpc). This result, first reported in Boeche et al. (2013) using RAVE-DR4 data (Kordopatis et al. 2013a), is not due to small number statistics (as approximately 5 · 10 4 targets are still available at those distances) nor to the non-homogeneity of the considered samples, as the metallicity gradient can also be measured when using LAMOST, APOGEE, RAVE-DR5, and GALAH separately (see top-left panel of Fig. 7). The inversion of the gradient (or the flattening) can be interpreted as a thick disc that is more centrally concentrated and more metal-poor than the thin disc combined with a thin disc that exhibits a flare at large radii. This result has also been suggested in studies in which the discs have been defined chemically ([α/Fe]-high population for the thick disc and Similar to the radial metallicity gradients that exhibit a vertical dependency, we find that the vertical metallicity gradients of the prograde stars also present a radial dependence, even when . Normalised metallicity distributions (obtained using a kernel density estimation with an Epanechnikov kernel and a smooting parameter of 0.06 dex) for stars close to the Sun (|R − R| < 0.2 kpc, |Z| < 0.2 kpc), colour-coded by survey. The plot has been truncated at -0.7, in order to better visualise the peak of the distributions. A good agreement is found between the different surveys in the sense that they are all peaking at the same value (∼ 0). We note that the SEGUE distribution contains less than 100 stars (most of them K dwarfs), resulting to a noisier shape.
individually considering the stars belonging to each survey (see bottom left panel of Fig. 7). The gradients we measure for all of the selected stars range from −0.24 dex kpc −1 at the inner disc to −0.17 dex kpc −1 at the outer disc. The vertical gradients that are derived for the range 5.3 ≤ R ≤ 9.2 kpc are similar to the gradients previously found in the literature, on the order of −0.2 dex kpc −1 to −0.27 dex kpc −1 (see Nandakumar et al. 2017, for a comparison between the different surveys), compatible with a mixture of two populations of different scale-heights, namely, h z,thin ∼ 300 pc and h z,thick ∼ 1000 pc (e.g. Gilmore & Reid 1983).

Metallicity gradients for the retrograde stars
As far as the trends of the retrograde stars are concerned, they are strikingly different than the ones found for the prograde stars. The radial gradients are always significantly negative, on the order of ∼ −0.04 dex kpc −1 , with no clear indication of a dependency with |Z|, at least up to |Z| ∼ 2 kpc, (even when the surveys are considered individually).
As far as the vertical gradients are concerned, these are much shallower than the ones derived for the prograde stars, yet they are still significant, on the order of ∼ −0.05 dex kpc −1 (see also bottom-right panel of Fig. 7 for values derived from each survey separately). A mild dependence on the radius might exist, when considering all of the stars simultaneously, in the sense that the outer radial bin seems to have a slightly flatter vertical gradient than the one measured at the inner Galaxy, regardless of the zeropoint parallax correction.
The absence, or small dependence, of the vertical metallicity gradient variations on the distance from the Galactic center and of the radial metallicity gradient on the distance from the Galactic plane is compatible with either a well-mixed population or a single population (coming from e.g. a single accretion) Article number, page 4 of 16  Radial gradients for the prograde stars (in black) and the retrograde stars (in red). Solid, dashed, and dotted lines correspond to selections for 0.2 ≤ |Z| < 1 kpc, 1 ≤ |Z| < 2 kpc, and 2 ≤ |Z| < 4.5 kpc, respectively. Error bars are computed as σ/ √ N .
that would have a pre-existent metallicity gradient (e.g. Abadi et al. 2003). Further investigation is therefore needed in order to better understand this retrograde population. In the next section, we scrutinise the correlations between v φ and the metallicity.

Correlations between rotational velocity and metallicity
overlaps by half a step with the previous one. The prograde stars are shown on the right-hand side panels and the retrograde stars on the left-hand side panels. We note that these trends do not change (although they become more noisy due to fewer stars) when only one survey is taken into account at once (which is in 6. Vertical gradients for the prograde stars (in black) and the retrograde stars (in red). Solid, dashed, and dotted lines correspond to selections for 7.2 ≤ R < 8.2 kpc, 5.2 ≤ R < 7.2 kpc, and 9.2 ≤ R < 11.2 kpc, respectively. agreement with Nandakumar et al. 2017, who state that the selection functions of the considered surveys do not significantly alter the measured Galactic gradients), nor when the zero-point offset in parallax is corrected for.
We also note that the usual way of plotting the chemokinematic correlations that can be found in the literature is the opposite than the one presented in these figures: commonly, it is v φ that is marginalised over metallicity-bins. However, by doing so, the tails of the velocity distribution are systematically missed or ignored, which is not what we are looking to do in this study. For this reason, the plots that follow are not obtained in the "usual" way, that is,

Trends for the prograde stars
Regarding the prograde stars, we can find the expected and known trends in the disc: close to the plane (purple colours on Article number, page 5 of 16 the right panels of Figs. 8 to 10), stars show a negative correlation down to v φ ∼ 180 km s −1 (the typical velocity-lag of the thick disc). Then they show a positive correlation for lower angular momenta.
The negative correlation, typically found for the chemically defined thin disc stars (e.g. Lee et al. 2011;Recio-Blanco et al. 2014;Allende Prieto et al. 2016;Kordopatis et al. 2017) is due to the blurring of the stellar orbits: older thin disc stars have a higher eccentricity than young stars via more Lindblad resonances with the spiral arms; this eventually allows them to visit radii that are far from their birth radius (Lynden-Bell & Kalnajs 1972). As a consequence, old outer thin disc stars can reach the solar neighbourhood at their pericentre, hence with velocities lower than the local standard of rest (LSR), and the inner thin disc stars can reach the solar neighbourhood at their apocentre, hence with velocities higher than the LSR. Because of the radial metallicity gradient in the thin disc, stars that are more metalrich than the locally born stars tend to move faster than the LSR, and the stars that are more metal-poor than the locally born stars will tend to move slower than the LSR. Comparing the metallicities at which the inflexion happens in the Figs. 8, 9, and 10, we can notice that it shifts from super-solar metallicities at the inner disc to sub-solar metallicities at our outer disc sample, which is as expected for a disc with a negative radial metallicity gradient.
The positive correlation that we find for stars with v φ 160 − 180 km s −1 (typically found for the thick disc stars, see Spagna et al. 2010;Kordopatis et al. 2011;Lee et al. 2011;Recio-Blanco et al. 2014;Kordopatis et al. 2017;Re Fiorentin et al. 2019) is often used as an argument against radial migration in the thick disc. We note, however, that Minchev et al. (2019) suggested that this correlation may be due to the superposition of stars of different ages, each mono-age population itself having a negative correlation (the so-called Yule-Simpson effect). The precision of the ages available for our sample does not allow us to either support or reject this statement.
That said, we find that the trends become flatter as we move farther from the plane. Eventually, for 6 < |Z| < 10 kpc (or even |Z| > 4 kpc, when the zero-point offset in parallax is taken into account), that is, where the canonical thick disc still represents a significant fraction of the stars (assuming a scale-height of 1 kpc and a local normalisation of ∼ 0.1), the trends seem to become negative for all of the stars over all the v φ -range (i.e. even at large velocities where the ratio canonical thick disc or canonical halo is high) and at all Galactocentric radii. This result is compatible with Minchev et al. (2019) and could potentially suggest that the thick disc stars far from the plane and the inner halo might be one single mono-age population. To our knowledge, this flattening trend has not been identified in any previous study.

Retrograde stars
Unlike the prograde stars, all of the retrograde stars seem to exhibit similar behaviour at any distance from the plane and any distance from the Galactic centre. The correlation is always positive, on the order of ∼ 25 − 30 km s −1 dex −1 , again suggesting that the probed population is well-mixed at all R (one would otherwise expect different trends). However, the retrograde sample contains, rather unexpectedly, super-solar metallicity (SMR) stars (see for instance Figs. A.2 and A.3 as well as Fig. 8, left panel, around v φ ∼ −325 km s −1 ). Although only ∼ 600 of them (∼ 300 when the zero-point offset in parallax is taken into account, observed mostly by the APOGEE, LAMOST, and RAVE surveys), they show an inverse (i.e. negative) correlation with metallicity, similar to the thin disc prograde stars. These SMR stars are located at all Galactic sky coordinates ( , b), with 50 per cent of them being closer than 1 kpc from the Galactic plane, and reaching distances up to |Z| ∼ 6 − 7 kpc. They are also seen in both the inner and the outer Galaxy (up to R ∼ 15 kpc). Interestingly, when the zero-point offset in parallaxes is not taken into account, 20% of the retrograde stars are located in the bulge region, that is, between 1 < R < 2.5 kpc, but the latter sample becomes prograde when correcting for the zero point (also see the plots in Figs. A.2 and A.3).
We cross-matched our retrograde sample with the APOGEE-DR16 catalogue (Ahumada et al. 2020) and out of the 1197 stars which have a reliable abundance determination 4 , we find that 18 have [M/H] > 0 (see Fig. 11) 5 . A specific investigation of the stellar spectra belonging to the APOGEE, LAMOST and RAVE surveys will be performed in future papers with the aim of confirming whether those SMR retrograde stars are indeed as metalrich as suggested by the the Sanders & Das (2018) catalogue. We show, nevertheless, in Fig. 12, the APOGEE DR16 spectra and their best-fit templates for five of those SMR retrograde targets as a preview of the study that will follow, also indicating that these stars have a good parameterisation.
Recently, Fragkoudi et al. (2020), Belokurov et al. (2020) and Grand et al. (2020) revived the idea that the thick disc may have formed after a violent gas-rich merger (e.g. Brook et al. 2004). In this context, because of the compressible nature of the gas, some of the newborn stars could have been formed on counter-rotating or radial orbits. The time and mass of this gas-         Fig. 11. Scatter plot of magnesium abundance as a function of metallicity for the retrograde stars in the APOGEE-DR16 catalogue and marginalised normalised histograms. Black, red, and green points and histograms correspond to counter-rotating stars that we labelled 'Gaia-Enceladus-Sausage', 'Thamnos', and 'other', respectively (the criteria for selecting those stars are described in Sect. 5). The sizes of the points are indicative and correspond to the estimated ages from Sanders & Das (2018); we note, however, that the majority of these stars are giants and, hence, their ages are unreliable. rich merger, according to Grand et al. (2020), could be then inferred from the properties and the fraction of the locally born retrograde stars. Figure 11 shows the magnesium abundance as a function of metallicity for all of the retrograde stars in the APOGEE DR16 sample. We can see from this figure two chemical sequences starting at [M/H] −1.5: a high-α, typically associated with stars born in situ, and a low-α, associated with accreted populations (e.g. Nissen & Schuster 2010). Whereas the majority of the counter-rotating stars are located in the low−α sequence, a non-negligible fraction of them, extending up to super-solar metallicities, are on the high-α sequence, indicating that they are possibly born locally. As we do not know the selection function of our sample, we cannot draw conclusions on their relative proportions. However, the fact that they are slightly α-enhanced, even for [M/H] > 0, corroborates to the fact that they were likely born from in situ material.
According to Helmi (2020); Naidu et al. (2020), the largest fraction of retrograde low-α stars belong to Gaia-Enceladus-Sausage, especially at large distances from the plane. In the scenario where Gaia-Enceladus-Sausage was a massive merger (with a mass-ratio of 4:1), this could imply that such a massive galaxy would have an internal metallicity gradient and puffed-up the Galactic disc that was present at the time of the merger. The low-metallicity accreted stars from such a merger would therefore be deposed first, on high eccentricity, and the metal-rich stars (formed close to the center of the accreted galaxy) would be deposited on more circular orbits due to dynamical friction (e.g. Koppelman et al. 2020). In the next section, we investigate more closely the case of Gaia-Enceladus-Sausage and try to identify differences with the other retrograde stars and groups of accreted stars, as identified recently in the literature.

Gaia-Enceladus-Sausage, Thamnos/Sequoia, ωCen, and the other counter-rotating stars
In what follows, we define three sub-samples of counter-rotating stars.
First, we labelled, as 'Gaia-Enceladus-Sausage' the targets fulfilling the quality criteria presented in Sect. 2 and having L Z < 0 kpc km s −1 and v φ > −200 km s −1 as well as e > 0.65 (16 283 out of the 22 472 counter-rotating stars, LAMOST and APOGEE targets encompassing 56.4 and 16.1 per cent of the targets, respectively). This selection is somewhat compatible with remnants of a counter-rotating merger of initial inclination of approximately 30 deg (see Helmi 2020, Sect. 4.2.1) but we stress that it ignores the prograde counter-part of the merger remnants and is likely contaminated by other accreted populations (see for example Feuillet et al. 2020). Therefore, our sample should not be directly compared with other Gaia-Enceladus or Gaia-Sausage samples in the literature.
We also selected stars probing the chemically and dynamically peculiar population identified by Koppelman et al. (2019) and dubbed 'Thamnos', the targets fulfilling the previous quality criteria and having L Z < 0 kpc km s −1 , e < 0.65, as well as v φ > −200 km s −1 (4 007 stars, LAMOST and APOGEE encompassing 52.1 and 19.2 per cent of the targets, respectively), that is, the stars having similar angular momentum and v φ as our 'Gaia-Enceladus-Sausage' sample, but having low eccentricities. Despite being still unclear if 'Thamnos' is part of the low-energy tail of Sequoia (see Myeong et al. 2019;Naidu et al. 2020), in what follows we nevertheless labelled this sample 'Thamnos' as the majority of Sequoia stars are on higher energy orbits.
Finally, all the other retrograde stars that are not 'Gaia-Enceladus-Sausage' nor 'Thamnos', with L z < 0 kpc km s −1 , are labelled in what follows as "other" (1 461 stars, LAMOST, RAVE and APOGEE encompassing 58.1, 13.6, and 13.1 per cent of the targets, respectively). It is noteworthy that this category comprises, 'Sequoia' stars as well as other sub-populations (see e.g. Naidu et al. 2020). Figures 13 and 14 show the radial and vertical metallicity gradients for those three samples. We did not separate them into different vertical or radial ranges, as we did in Sect. 3.2, since we found that the gradients for the retrograde stars are not particularly dependent on spatial cuts. We find that the 'Gaia-Enceladus-Sausage' sample is globally more metal-rich than the 'Thamnos' sample (in agreement with Koppelman et al. 2019;Mackereth et al. 2019) and that the two former samples are, in turn, more metal-rich than the 'other' counter-rotating sample in our sample.
When considering all of the surveys simultaneously, the radial and vertical metallicity gradients that we measure for 'Gaia-Enceladus-Sausage', 'Thamnos', and the 'other' sub-samples are similar, within 1σ and 2σ, respectively (see Table 3), yet this result should be taken with a grain of salt as it is not necessarily corroborated when considering LAMOST or APOGEE targets separately (maybe due to small number statistics or different volume coverage).
This lack of evidence of difference in the metallicity gradients, other than the zero-point offset, between each of the counter-rotating sub-samples, made us investigate how the metallicity changed as a function of the eccentricity. Results are shown in Fig. 15. Strikingly, we find that 'Gaia-Enceladus-Sausage' and 'Thamnos' samples follow the same trend, as opposed to the 'other' subsample, which appears to show no variation of metallicity as a function of the eccentricity (a result that is also found when separating the targets per survey). This  Fig. 12. APOGEE DR16 spectra (in a selected arbitrary wavelength range, in black) and best-fit templates (in red) for five of the super-solar metallicity retrograde stars in our sample. The Gaia sourceid is reported at the top-right corner and the atmospheric parameters of the stars, as well as their azimuthal velocity, v φ , at the top-left corner.   trend seems to favour the fact that the 'Thamnos' and 'Gaia-Enceladus-Sausage' samples are strongly linked, yet the question remains about whether they originate from the same parent population.
The [Mg/Fe]−[M/H] plot coming from the APOGEE-DR16 data (Fig. 11) shows not only that our 'Gaia-Enceladus-Sausage' and 'Thamnos' samples are mixed populations, containing both high-α and low-α stars, but also that they exhibit different distributions in both the metallicity space and the [Mg/Fe] space (see coloured dots and histograms). The two sample Kolmogorov-Smirnov tests of the metallicity and [Mg/Fe] between the 'Gaia-Enceladus-Sausage' and 'Thamnos' samples strongly suggest that the two samples cannot have been drawn from the same underlying distribution, with p−values always lower than 10 −7 , even when decomposed into radial and vertical spatial bins.
Finally, we also explored the possibility of having a significant amount of stars belonging to the retrograde (and peculiar on many aspects) globular cluster ωCentauri ( Gratton et al. 2011;Marino et al. 2011). Figure 16 shows no clear indication of anything of the sort for either one of the three subsamples, suggesting that, globally, our sample of retrograde stars does not contain any obvious sub-population of stars stripped from the ωCen progenitor. We note, however, as remarked by Zasowski et al. (2019), that sodium abundances in APOGEE exhibit a large scatter compared to the published uncertainties, partly due to the presence of strong telluric absorptions close to the Na lines. Therefore, firmer conclusions on this topic might require a more in-depth investigation of the abundance patterns of APOGEE, which is beyond the scope of this paper.

Discussion and conclusions
The merger-tree of the Milky Way is a matter of vivid debate. Despite the recent discoveries of many over-densities in the action-energy space in the extended Solar neighbourhood, it is still rather unclear how these populations relate to each other, or to what extent they are well-mixed within the Galaxy. In this paper, we use a compiled catalogue of the metallicities and velocities of more than 4 million stars from Sanders & Das (2018) to investigate the chemical gradients and the chemodynamics of the counter-rotating stars, while also performing similar analyses for the prograde stars. This parallel analysis allowed us to highlight the differences between the two populations and also to confirm that our analysis method returned sound results.
Our main results can be summarised as follows: -We determined that stars on prograde orbits show a v φ − [M/H] relation that flattens -and perhaps even reverses as a function of distance from the plane. To our knowledge, this trend has not yet been reported in the literature. We investigated potential biases due to parallax zero-point offsets or the inhomogeneity of the survey metallicities in order to make sure that the results of our analysis are robust. Globally, a correction of the parallax zero point affects mainly the retrograde stars, increasing their velocities and, hence, decreasing the sample-size of this population. Our conclusions, however, remain robust as they do not rely on the absolute densities of the populations. The different zero points in metallicity between the surveys do not alter our conclusions either. Indeed, where, for example, the absolute values of the metallicity gradients might differ if a specific survey is used instead of another, the global trends remain unchanged. That said, we stress that additional biases might still exist due to the different quality cuts that we applied to our input catalogue. These are difficult to quantify (e.g. the effect of filtering for the stars with large uncertainty in metallicity), still, we present a brief discussion in Appendix C.
From the points above, the following picture can be drawn. The flattening (or inversion of sign, from positive to negative) of the correlation of v φ with metallicity for the prograde stars far from the plane may suggest, in agreement with other studies, that the inner halo and the thick disc far from the plane may be populations that share a common past. Indeed, it is believed that mono-age disc populations have intrinsically negative v φ − [M/H] correlations . Hence, going from a positive correlation for the thick disc stars close to the plane to a negative (or flat) one for the stars at large distances from it, where the probed population is a mixture of thick disc and halo, would indicate that the targeted population has a smaller age range.
One point that is worth mentioning here relates to recent results, for example, from Haywood et al. (2018) Belokurov et al. (2020), who revived the idea that the low-angular momentum thick disc might be made up of stars from the proto-disc that are now part of the inner halo after having been heated by a major past accretion (see also, Gilmore et al. 2002;Wyse et al. 2006;Kordopatis et al. 2013b). In particular, using the same input catalogue as we have in this work, Belokurov et al. (2020) identified a small, yet non-negligible, prograde population of intermediate metallicity ([M/H] −0.7) with low angular momentum (v φ 100 km s −1 ), which they have dubbed 'the Splash'. These stars are found to be slightly younger than the accreted ones from the last major merger (i.e. Gaia-Enceladus-Sausage) and with a different star-formation history than the latter. These aspects led these authors to suggest that the Splash stars were born locally in the protodisc and have had their orbits altered by this massive accretion (see, however, Amarante et al. 2020, for an alternative explanation). Our analysis, in selecting each time only the prograde or the retrograde stars, unavoidably dilutes the signature of the Splash with either the thick disc or the halo and Gaia-Enceladus-Sausage targets. However, investigating Figs. 8 to 10, at around 100 v φ 150 km s −1 (where Belokurov et al. 2020, suggest that there is the overlap between the Splash and the disc), we do not see any peculiarities in the trend other than the well-known change of behaviour in the v φ − [M/H] space, going from an anti-correlation to a positive correlation.
As far as the counter-rotating stars are concerned, we find that more retrograde targets are also, on average, more metalpoor. The lack of spatial changes in the metallicity gradients, or in the v φ − [M/H] correlations, within the probed volume (5.2 < R < 11.2 kpc and |Z| < 10 kpc), indicate that the counter-rotating stars are remarkably well-mixed. This is in line with, for example, Deason et al. (2018); Helmi (2020) (and with hints of this already discussed in Watkins et al. 2009;Deason et al. 2013) suggesting that the orbital turning points (i.e. shells) of Gaia-Enceladus-Sausage are at a distance of ∼ 20 kpc, that is, outside the volume that we investigate here. Our results expand upon this, as they show that none of the shells of Gaia-Enceladus-Sausage, Thamnos, Sequoia, or any other accreted population that could compose the retrograde stars, is within 10 − 12 kpc from the Sun.
Amongst the counter-rotating stars, the samples we dub 'Thamnos' and 'Gaia-Enceladus-Sausage' seem to compose two different populations that are, nevertheless, also somehow linked. Indeed, we find that the 'Thamnos' sample is more metalpoor than the 'Gaia-Enceladus-Sausage' one and on more circular orbits. However, it is striking to observe that they follow the same trend in the metallicity-eccentricity space. This trend cannot be explained simply by the presence of intrinsic metallicity gradients within a common progenitor of Gaia-Enceladus-Sausage and Thamnos/Sequoia. Indeed, despite the initial mass estimates of Gaia-Enceladus-Sausage found in the literature, ranging from 6·10 8 to 5·10 9 M (e.g. Helmi et al. 2018;Myeong et al. 2019;Kruijssen et al. 2019;Mackereth et al. 2019;Fattahi et al. 2019;Feuillet et al. 2020), which could lead to such internal radial metallicity gradients, it is unclear how a counter-rotating accretion with an angle of ∼ 30 deg (see Helmi 2020) would deposit its most metal-rich stars, which would otherwise be expected to be present at the innermost regions of the progenitor at the highest eccentricities.
We believe that the key to understanding the effect of the past accretions on the properties of the thick disc is the super-solar metallicity counter-rotating population that we identified in our sample. These stars, if proven to really exist, could have been formed from local gas at the moment of the accretion with the Gaia-Enceladus-Sausage progenitor (see Maiolino et al. 2017;Gallagher et al. 2019, for star formation activity within galactic outflows). This population would therefore offer us an undeniable sample of locally born retrograde stars in order to precisely date and weigh this merger (see, Grand et al. 2020). The Gaia data release 3, which will become public at the end of 2020, will first confirm whether these peculiar stars are indeed counter-rotating. Then, future large spectroscopic surveys such as WEAVE (Dalton et al. 2018) or 4MOST (de Jong et al. 2019) will allow us to verify the chemical composition of these stars, investigate them in more dimensions of the chemical space, and, eventually, to better understand their link with the rest of the galactic populations.

Acknowledgments
We thank the anonymous referee for their useful feedback that helped improving the quality of this paper. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/ gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research made use of Astropy, 6 a community-developed core

Appendix B: Comparison between repeats in several surveys
We show in Fig. B.1, the comparison between the metallicity, the distance, and the azimuthal velocity, as derived by Sanders & Das (2018), for the repeated inter-survey stars. A very good agreement is obtained for distances and v φ (largest median offsets are 10 ± 20 pc and 0.5 ± 2.7 km s −1 , respectively), whereas metallicity offsets are smaller or equal to 0.10±0.16 dex (largest median offset found is between RAVE and GALAH). Since we keep only the stars that exhibit a very good astrometry and a small uncertainty in metallicity (see Sect. 2), this offset in [M/H] mostly translates the offset in metallicity derived spectroscopically by each survey. Eventually, we kept only one entry for those repeated stars, with the following order of preference: APOGEE, GALAH, GES, RAVE, LAMOST, and SEGUE.