Exploring the Galaxy’s halo and very metal-weak thick disc with SkyMapper and Gaia DR2 Free

ABSTRACT

In this work, we combine spectroscopic information from the SkyMapper survey for Extremely Metal-Poor stars and astrometry from Gaia DR2 to investigate the kinematics of a sample of 475 stars with a metallicity range of |$-6.5 \le \rm [Fe/H] \le -2.05$| dex. Exploiting the action map, we identify 16 and 40 stars dynamically consistent with the Gaia Sausage and Gaia Sequoia accretion events, respectively. The most metal poor of these candidates have metallicities of |$\rm [Fe/H]=-3.31\, \mathrm{ and }\, -3.74$|⁠, respectively, helping to define the low-metallicity tail of the progenitors involved in the accretion events. We also find, consistent with other studies, that ∼21 per cent of the sample have orbits that remain confined to within 3 kpc of the Galactic plane, that is, |Z_max| ≤ 3 kpc. Of particular interest is a subsample (∼11 per cent of the total) of low |Z_max| stars with low eccentricities and prograde motions. The lowest metallicity of these stars has [Fe/H] = –4.30 and the subsample is best interpreted as the very low-metallicity tail of the metal-weak thick disc population. The low |Z_max|, low eccentricity stars with retrograde orbits are likely accreted, while the low |Z_max|, high eccentricity pro- and retrograde stars are plausibly associated with the Gaia Sausage system. We find that a small fraction of our sample (∼4 per cent of the total) is likely escaping from the Galaxy, and postulate that these stars have gained energy from gravitational interactions that occur when infalling dwarf galaxies are tidally disrupted.

1 INTRODUCTION

In the last decade, the astronomical community has experienced a renewal of interest in the properties of low-metallicity stars, particularly those with |$\rm [Fe/H]\lt -2\, dex$|⁠.¹ Motivated by the successful surveys of Beers, Preston & Shectman (1992) and Christlieb et al. (2008), in recent years numerous spectroscopic (e.g. SDSS, SEGUE, LAMOST, and APOGEE, York et al. 2000; Yanny et al. 2009; Cui et al. 2012; Zhao et al. 2012; Majewski et al. 2017) and photometric (e.g. Pristine and SkyMapper, Starkenburg et al. 2017; Wolf et al. 2018) surveys have been commissioned, scanning extensive sky areas for these very rare and key objects. We refer to Da Costa et al. (2019, their section 1) for a more complete list of spectro-photometric surveys targeting low-metallicity stars. Not surprisingly, the underlying scientific motive is the understanding of the formation of our Galaxy, as well as other galaxies in the Universe.

Specifically, the lowest metallicity stars observable at the present-day formed from gas enriched with the nucleosynthetic products from first generation metal-free stars, the so-called Population-III stars. Studies of abundances and abundance ratios in ultra- (UMP) and extremely metal-poor (EMP) stars can then yield constraints on the properties of the Pop III stars, such as their masses, and on star formation processes at the earliest times (e.g. Frebel & Norris 2015). Moreover, the kinematics of these stars can also provide much information on the events that occurred during the formation of the Milky Way (MW), which is believed to include both star formation in situ and the accretion of lower mass galaxies. Indeed, together the abundances and kinematics of the lowest metallicity stars offer a distinct perspective on the earliest stages of the formation and evolution of the MW, and by implication, of galaxies in general.

In terms of the formation of the MW, the most common scenario predicts that the most metal-poor stars will be found mainly in the Galactic halo and bulge, as these components likely formed in the earliest stages of the MW’s evolution (e.g. White & Springel 2000; Brook et al. 2007; Tumlinson 2010; El-Badry et al. 2018). In such a scenario relatively few, if any, very metal-poor stars are expected to lie in the MW disc as it formed at a later epoch after the settling into the plane of gas enriched by multiple generations of star formation (e.g. Bland-Hawthorn & Gerhard 2016). However, recent kinematic results from surveys for the most metal-poor stars have cast doubt on this scenario, altering our understanding of the formation of the MW. For example, the recent studies of Sestito et al. (2019, 2020b), Di Matteo et al. (2020), and Venn et al. (2020) have revealed a new scenario where |$\sim \, 20\, {{\ \rm per\ cent}}$| of very metal-poor stars have orbits that are confined to within 3 kpc of the MW plane; evidently the majority of these stars are not Galactic halo objects despite their low metallicities.

In particular, Sestito et al. (2019) compiled a catalogue of 42 UMP ([Fe/H] ≤ –4.0) stars from the literature and analysed their orbital properties making use of Gaia DR2 proper motions (Gaia Collaboration et al. 2018). They found that 11 out of 42 stars have prograde orbits that are confined to within 3 kpc of the MW disc. Moreover, two of these MW-planar stars are found to be on nearly circular prograde orbits, and one is the star with the lowest overall metal content currently known (Caffau et al. 2011). In the same fashion, Di Matteo et al. (2020) investigated the kinematics of a sample of coincidentally the same number of low-metallicity stars drawn from the ESO Large Program ‘First stars – First nucleosynthesis’ (Cayrel et al. 2004). Their analysis also finds that |$\sim 20{{\ \rm per\ cent}}$| of the stars show disc-like kinematics. They went on to consider a larger sample of stars covering a wider metallicity range and found consistent results. Di Matteo et al. (2020) then postulated the existence of an ‘UMP thick disc’ that is an extension to low metallicities of the Galaxy’s thick disc population.

Sestito et al. (2020b) carried out a similar kinematic analysis on a substantially larger sample, consisting of 1027 very metal-poor stars with [Fe/H] ≤ –2.5 selected from the Pristine (Starkenburg et al. 2017; Aguado et al. 2019) and LAMOST (Cui et al. 2012; Li et al. 2018) surveys. Again they find that almost one-third of the stars in the sample have orbits that do not deviate significantly from the disc plane of the Galaxy. They suggest that this implies that a significant fraction of the MW’s metal-poor stars formed with the MW (thick) disc. Moreover, they note that as a consequence, the history of the disc must have been sufficiently quiescent that (presumably old) metal-poor stars were able to retain their disc-like orbits to the present day (Sestito et al. 2020b).

Venn et al. (2020) have also investigated the kinematics of metal-poor stars using a sample of 115 objects chosen from the Pristine survey (Starkenburg et al. 2017) that have been observed at high dispersion. They find 16, out of 70, metal-poor stars whose orbits are confined to the vicinity of the Galactic plane, together with small numbers of stars that may have unbound orbits. They also identify stars whose orbital characteristics/actions are consistent with an origin in the Gaia Enceladus (Helmi et al. 2018) accretion event.

These somewhat unexpected results support the idea that the metallicity distribution of the Galaxy’s thick disc does indeed possess a low-metallicity tail, as first advocated by Norris, Bessell & Pickles (1985) and Morrison, Flynn & Freeman (1990). Moreover, the proposed low-metallicity tail would extend to lower metallicities than those authors suggested (see also Chiba & Beers 2000; Beers et al. 2014).

The origin(s) of these metal-poor thick disc stars is, however, still uncertain, though the implications of their existence for the formation and evolution of the MW, and disc galaxies in general, are likely significant. A number of different possibilities have been discussed (e.g. Sestito et al. 2019, 2020b; Di Matteo et al. 2020) including that the stars were accreted from small satellites once the MW disc had already formed, or that they represent low-metallicity stars formed in the gas-rich building-blocks that came together to form the main body of the Galaxy’s disc (see also the theoretical simulations presented in Sestito et al. 2020a).

In this work, we conduct a similar study to those mentioned above by exploiting the metallicity determinations from the SkyMapper Survey for extremely metal-poor stars (see Da Costa et al. 2019), together with Gaia DR2 astrometry (Gaia Collaboration et al. 2018), to investigate the dynamics of 475 very metal-poor ([Fe/H] < –2) stars in the southern sky. The wide extension in metallicity space, together with the relatively large number of stars, gives us a detailed view of the kinematic properties of these objects. We also consider the potential connection of any of the stars in our sample with the MW accretion events, such as those designated Gaia Enceladus, Gaia Sausage, and Gaia Sequoia that have been recently discovered in large-scale analyses of Gaia DR2 data (e.g. Helmi et al. 2018; Belokurov et al. 2018; Myeong et al. 2019; Mackereth et al. 2019). Such a connection has also been pursued in Monty et al. (2020).

The paper is organized as follows: in Sections 2 and 3 we present the data set and the orbit determination procedure, respectively, while in Sections 4 and 5, we present and discuss our results. Specifically, in Section 5.3, we discuss the small number of stars in our sample that appear not to be bound to the Galaxy. The final section (Section 6) summarizes our findings.

2 DATA

The data set used in this work consists of 475 stars with metallicities ranging from |$\rm [Fe/H]=-2.08\, \mathrm{ to}\, \lt -6.5$| dex. It is composed as follows:

• 114 giant stars with |$-6.2 \le \rm [Fe/H]_{1D,\, LTE} \le -2.25$| dex. Of these stars 113 come from Yong et al. (in preparation), while the remaining star is the most-iron poor star for which iron has been detected: SMSS J160540.18–144323.1 with |$\rm [Fe/H]_{1D,\, LTE} = -6.2 \pm 0.2$| (Nordlander et al. 2019). These stars originate with the EMP candidates discussed in Da Costa et al. (2019) and all have been observed at high resolution, principally with the MIKE spectrograph (Bernstein et al. 2003) at the 6.5-m Magellan (Clay) telescope. We shall refer to these stars as the HiRes data set.

• 45 stars observed with the FEROS high-resolution spectrograph (Kaufer et al. 1999) at the MPG/ESO 2.2-m telescope at La Silla. Again, these stars originated from the Da Costa et al. (2019) sample. We removed from the analysis all the stars with |$\rm [Fe/H]\gt -2$|⁠, and the stars in common with HiRes data set. The final count of stars belonging to this subsample is 38 and we label it as the FEROS data set.

• 122 stars from Jacobson et al. (2015) which have −3.97 ≤ [Fe/H] ≤ −1.31 dex. These stars originated in the SkyMapper commissioning-era survey (see Da Costa et al. 2019), and were also observed at high-dispersion with the MIKE spectrograph at Magellan. As for the FEROS sample, we removed seven stars with |$\rm [Fe/H]_{1D,\, LTE} \gt -2$| and the single star in common with the HiRes data set. However, as discussed in Section 3.2, there appears to be an issue with the radial velocities for the stars observed by Jacobson et al. (2015) during one specific Magellan/MIKE run, namely 2013 May 28–June 01. As a result, we have removed the stars observed in that run that lack a radial velocity from Gaia DR2 and which had not been already discarded. The final subsample used here is then composed of 91 stars and we refer to it as the Jacobson+15 subsample.

• 17 stars from Marino et al. (2019) with metallicity |$-3.26\rm \lt [Fe/H]_{1D,\, LTE}\lt -1.71$| dex. The spectra of these stars were obtained with the Keck HIRES high-resolution spectrograph (Vogt et al. 1994). After the removal of two stars present in the Jacobson+15 subsample, and two stars with |$\rm [Fe/H]_{1D,\, LTE} \gt -2$|⁠, we retain 13 stars. This subsample is referred to as the Marino+19 data set.

• 362 giant star candidates from Da Costa et al. (2019) with either |$\rm [Fe/H]_{\rm fitter}$|²<−3.0, or |$-3.0 \le \rm [Fe/H]_{\rm fitter} \le -2.5$| and g_SkyMapper < 13.7 mag. The radial velocities from the low-resolution spectra lack sufficient precision for our analysis, so the list of stars was cross-matched with Gaia DR2 to obtain radial velocities. A total of 195 stars were retained after the cross-match. These stars are referred to as the LowRes data set.

• 24 UMP giant stars ([Fe/H] ≤ −4) from Sestito et al. (2019), included to increase the number of UMP stars in the full sample and to provide a consistency check on our procedures. We have specifically selected only known giants from their sample for consistency with the SkyMapper-derived samples, which are giant dominated. We refer to Sestito et al. (2019) for a detailed description of the data set but we note it includes the star SMSS J031300.36−670839.3, which has [Fe/H]_{3D, NLTE} < −6.5 (Keller et al. 2014; Bessell et al. 2015; Nordlander et al. 2017). This data set is referred to as the Sestito+19 subsample.

Unless otherwise noted, the uncertainty in [Fe/H] values derived from high dispersion spectroscopy is taken as ±0.10, while for the stars in the LowRes data set, the uncertainty is ±0.3, and the values are quantized at 0.25 dex intervals. Fig. 1 then shows the metallicity distribution of each data set, computed using kernel density estimation with a Gaussian kernel and a bandwidth parameter of 0.5; the number of stars belonging to each set is reported in the top left corner of the panel. Each distribution has been normalized by the number of stars in the sample. Fig. 1 also shows the distribution for the total sample formed by summing the individual distributions. As is apparent, the sample spans a wide range in metallicity, with a peak around |$\rm [Fe/H]\sim -2.8$|⁠, consistent with the observed metallicity distribution function of the full SkyMapper EMP sample discussed in Da Costa et al. (2019).

Fig. 2 shows the position of the analysed stars, both in Galactic latitude and longitude and in the Cartesian Galactocentric reference frame, with each star colour-coded according to its metallicity. Since |$\sim \, 90{{\ \rm per\ cent}}$| of the stars come from the SkyMapper survey, the data set is affected by the same selection biases as discussed in Da Costa et al. (2019). Specifically, the SkyMapper survey avoids regions of the sky with significant stellar crowding, while the selection process for candidates restricts the sample to stars with E(B − V) < 0.25 mag. The net result is a lack of candidates near the Galactic plane and in the Galactic Bulge (see Fig. 4 and fig. 14 in Da Costa et al. 2019) as is evident in the inset in the middle panel of Fig. 2. Indeed, the majority of the stars lie inside the solar circle in the (⁠|$\rm {X}_{G},\rm {Y}_{G}$|⁠) plane, although at a variety of heights above and below the plane; the star nearest the Galactic Centre (GC) in the sample has a Galactocentric radius of 1.8 ± 0.8 kpc.

3 DERIVING THE KINEMATICS OF THE SAMPLE

To compute the orbit of a star the full 6D information for the position and velocity is needed. Specifically, we need right ascension (α), declination (δ), distance from the Sun (d), proper motions in right ascension and declination (μ_αcos δ, μ_δ), and the heliocentric radial velocity (v_r). Gaia DR2 provides the optimal source for these parameters, noting that strictly Gaia provides a measurement of parallax, not distance, and that v_r is only available for the brightest stars.

As recently discussed in Bailer-Jones et al. (2018), for example, simply inferring the distance from the parallax measurement alone can lead to unreliable results. To overcome this problem, Bailer-Jones et al. (2018) combined parallax measurements with a realistic prior for the distance as a function of Galactic longitude and latitude, to generate distance estimates.

Sestito et al. (2019) introduced an alternate approach to determining distances that combines the exquisite astrometry and photometry provided by Gaia DR2 with theoretical isochrones in a Bayesian analysis to infer the distance, as well as the physical properties surface gravity |$(\log \, g)$| and effective temperature (T_eff). An advantage of their technique is that it allows the breaking of the potential degeneracy between dwarf and giant star distances at a fixed T_eff. In our case however, by deliberate choice of the colour-range used to define the underlying sample of low-metallicity candidates in the SkyMapper EMP survey (see Da Costa et al. 2019), our data set consists entirely of giants,³ so that any dwarf/giant distance ambiguity does not arise. It further allows us to exploit the effective temperatures and metallicities of our stars, which are known from either the high-resolution analyses or from the spectrophotometric fits to the low-resolution spectra, to derive absolute magnitudes via the use of red giant branch (RGB) isochrones, particularly for those stars that lack a reliable parallax determination. This approach has the underlying assumption that all the stars lie on the RGB, whereas the distribution of temperatures and gravities in Da Costa et al. (2019) suggests a small fraction (5–10 per cent) of the total sample are red horizontal branch or early-asymptotic giant branch stars. Such stars are more luminous than RGB stars at the same effective temperature and thus the distance determinations based on the RGB locus will be smaller than the true distances. While this will result in some individually incorrect orbital parameters, the overall results are unaffected given the dominance of RGB stars in the sample.

Details of our distance determinations are discussed in the next section, but when we compare our derived distances with those in Sestito et al. (2019) for the stars in common, we find excellent agreement. This is illustrated in Fig. 3.

3.1 Distance determination

Our approach to determining distances for the stars in our sample is twofold. First, for the stars with unreliable Gaia DR2 parallax determinations, which we take here as those with σ_π/|π| ≥ 0.15, we adopted the following approach, which relies on the assumption that the stars in our sample, being very metal-poor, can be safely assumed to be old (age ≥ 10 Gyr).⁴ With this assumption we can then use the known T_eff and [Fe/H] values together with RGB isochrones of different metallicity to infer the absolute magnitudes and thus the distance. Specifically, we have used a set of Yonsei–Yale RGB isochrones⁵ (Y², Demarque et al. 2004) for an age 12 Gyr, [α/Fe] = +0.3, and metallicities corresponding to [Fe/H] = –3.5, –2.5, and –1.9 to infer the V-band absolute magnitude (M_V) for each star.

In practice, to find the absolute magnitude corresponding to a given star’s metallicity and T_eff, we interpolated in M_V across the isochrones at the T_eff value. Since the isochrones use visual magnitudes, we first calculated the appropriate V magnitude for each star from the Gaia G values using the coefficients provided by the Gaia documentation.⁶ Reddening values from Schlegel, Finkbeiner & Davis (1998) were adopted, corrected according to the recipe in Wolf et al. (2018). For stars with metallicities between −4.5 and −3.5, the (M_V) value is a linear extrapolation, while for the small number of stars with [Fe/H] ≤−4.5, which come primarily from the Sestito et al. (2019) subsample, the M_V inferred for [Fe/H] = −4.5 was used. The uncertainties in the distances were then determined by assuming an uncertainty of |$100\, \rm K$| in T_eff and |$0.1\, \rm dex$| in metallicity (⁠|$0.3\, \rm dex$| for stars in LowRes subsample) and then propagating these values into the distance determination.

Second, for the stars with nominally reliable Gaia DR2 parallax determinations, that is, those with σ_π/|π| < 0.15, we compared the Bailer-Jones et al. (2018) distances with the distances inferred from the RGB isochrones. This is shown in Fig. 4. While most stars do scatter about the 1:1 line, there are sizeable differences between the two estimates for ∼25 per cent of the stars, most commonly with the RGB-based distance being larger than the Bailer-Jones et al. (2018) value, indicating that the parallax may have been overestimated, or that the RGB-based distance is incorrect.

We have not sought to investigate the origin of the discrepancy for each individual case, noting that we include uncertainties in T_eff and [Fe/H] when estimating the uncertainty in the RGB-based distance. There is, however, a potential systematic uncertainty introduced by RGB isochrone based approach. In particular, as discussed by Joyce & Chaboyer (2018), the location of theoretical RGBs in the Hertzsprung–Russell diagram is sensitive to the adopted value of the mixing length parameter α_MLT. The value of α_MLT employed in any particular isochrone set (e.g. α_MLT = 1.7 for the Y² isochrones) is usually determined by requiring a fit to the solar values, but, as demonstrated in Joyce & Chaboyer (2018), at low metallicities the location of the RGB computed with a solar-calibrated α_MLT is more luminous by ∼0.3 mag at a fixed T_eff than a comparison with globular cluster RGB observations would suggest: a ∼10 per cent smaller value of α_MLT is required for consistency with the observations. It is possible therefore that our RGB-based distances are systematically overestimated, though the comparison shown in Fig. 4 suggests that it is not a major effect.

In practice, we have adopted the Bailer-Jones et al. (2018) distance and its uncertainty whenever the absolute value of the relative difference |$(\Delta = \frac{D_{\rm RGB}-D_{\rm BJ+18}}{D_{\rm RGB}})$| is smaller than 0.35. This is shown as the grey shaded region in Fig. 4. For the remainder, that is, for the stars outside the shaded area with |Δ| > 0.35, we adopted the distance inferred from RGB isochrones. Overall, this results in the use of the RGB isochrone distance for 357 stars, while for the remaining 118 the Bailer-Jones et al. (2018) distance is employed.⁷

The largest (heliocentric) distance of the stars in our sample is the RGB-based distance of ∼35 kpc for the high-luminosity giant (⁠|$\log g\, \approx$| 0.3) SMSS J004037.55–515025.1, which has [Fe/H] = –3.83 and is in the Jacobson+15 subsample. The smallest is the Bailer-Jones et al. (2018) distance of 0.21 kpc for the subgiant BD+44 493 (⁠|$\log g\, \approx$| 3.2 and [Fe/H] = –4.30) from the Sestito+19 subsample. Overall, the median heliocentric distance for the entire sample is ∼5 kpc, with a median ∼7 kpc for RGB-based distances, and a median of ∼3 kpc for Bailer-Jones et al. (2018) distances.

3.2 Orbital properties

To compute the orbital parameters, we used the full 6D information on the position and velocity for each star. Gaia DR2 provides coordinates and proper motions, while the distances have been obtained as discussed in the previous section. Radial velocities come from the high-dispersion spectra, when available, and from Gaia DR2 for the LowRes sample.

We note that there are a number of the stars with radial velocities from the high-dispersion spectra that also have radial velocities from Gaia DR2, and this allows us to check for anything unusual or unexpected. As mentioned in Section 2, in this comparison process we discovered an anomaly in the Jacobson et al. (2015) radial velocities for a particular Magellan/MIKE run. In that run 32 stars were observed of which eight also have radial velocities from Gaia DR2. The comparison for these eight stars shows extreme disagreement for seven stars, with values of the difference V_r(J+15) – V_r(Gaia DR2) ranging from –400 to +415 |$\rm km\, s^{-1}$|⁠. We are at a loss to explain the origin of the disagreements⁸ and have consequently excluded from the analysis the stars from this run that lack Gaia DR2 radial velocities, while using the Gaia DR2 radial velocities in the kinematic calculations for the remaining seven stars that have [Fe/H] ≤ –2.0 dex. We stress that such large disagreements are seen only for this one observing run in the Jacobson+15 sample, the radial velocities from other runs are very consistent with Gaia DR2 values when available. This is also the case for the stars in the HiRes and FEROS samples. Overall, for the 41 stars with radial velocities from our high-dispersion spectra and from Gaia DR2, the velocities agree well with a mean difference, in the sense of our velocities minus Gaia, of 1.7 |$\rm km\, s^{-1}$| and a standard deviation of 5.5 |$\rm km\, s^{-1}$|⁠. This agreement indicates that any systematic uncertainties in the radial velocity determinations are very minor compared to other contributors to uncertainties in the orbit determinations. We have always used the radial velocity and the corresponding uncertainty from the high-dispersion spectra when available; the Gaia radial velocities and their uncertainties were utilized only when there was no alternative.

The kinematics of our sample of metal-poor stars have been determined using the galpy⁹python package (Bovy 2015). The orbit of each star was obtained by direct integration backward and forward in time for 2 Gyr. This choice relies on the assumption that such a time-scale is shorter than any significant variation in the Galactic potential.

We adopted the potential identified as the best candidate among the ones studied in McMillan (2017).¹⁰ Briefly, it consists of an axisymmetric model with a bulge, thin, thick, and gaseous discs, and a Navarro–Frenk–White (Navarro, Frenk & White 1996) dark matter halo. The heliocentric distances derived in the previous section have been converted to distances from the GC in the galpy routine, specifying the galactocentric position of the Sun as |$(X, Y, Z)=(-8.21, 0, 0.0208)\, \rm kpc$|⁠, and its circular speed as |$v_0=232.8\, \rm km\, s^{-1}$|⁠. Both quantities are taken from McMillan (2017).

For each star, we determined the apogalacticon and perigalacticon |$(D_{\rm apo}, \, D_{\rm peri})$| of the orbit, the maximum vertical excursion from the Galactic plane (Z_max), the eccentricity |$\left(e=\frac{D_{\rm apo}-D_{\rm peri}}{D_{\rm apo}+D_{\rm peri}}\right)$|⁠, the energy (E), the three actions |$(J_{\rm R},\, J_{\rm \phi },\, J_{\rm Z})$|¹¹ and the velocity components |$U,\, V,\, W$| in the frame of the local standard of rest (LSR). As have others (e.g. Myeong et al. 2018), we emphasize that action space is the ideal plane in which to evaluate large samples of MW stars to identify and study possible substructures and debris from accretion events. The reason is that the actions are nearly conserved under the hypothesis that the potential is smoothly evolving (Binney & Spergel 1984).

As examples, we consider the two most iron-poor stars known: SMSS J031300.36–670839.3 (Keller et al. 2014; Bessell et al. 2015) and SMSS J160540.18–144323.1 (Nordlander et al. 2019). For the former we find an ‘outer-halo’ orbit with |$e = 0.70 \pm 0.05\, , D_{\rm peri} = 6.5 \pm 2.0,\, D_{\rm apo} = 36.6 \pm 9.8$|⁠, and |Z_max| = 34.2 ± 9.2 kpc. These parameters are in good agreement with those listed in Sestito et al. (2019). For the latter star, however, we determine an extreme ‘outer-halo’ orbit that may in fact be unbound as the derived energy E is close to zero. The inferred parameters are |$e = 0.93,\, D_{\rm peri} = 6.5 \pm 2.0, D_{\rm apo} \approx 423$|⁠, and |Z_max| ≈ 327 kpc; the latter two quantities are quite uncertain.

As discussed in detail in Section 5.3, SMSS J160540.18–144323.1 is, in fact, one of a small number of stars (30 out of 475) for which we find apparent apogalacticon distances larger than the MW virial radius, that is, larger than ∼250 kpc. For such stars, a substantial fraction of the 500 random realizations resulted in unbound orbits (i.e. D_apo = ∞), thus potentially biasing both the medians and the uncertainties derived from the orbital parameter distributions. The uncertainties for these specific stars are considered in more detail in Section 5.3.

As an independent check on the uncertainties and on the role of the adopted potential, we can compare our orbit parameters with those listed in Sestito et al. (2019, see their table 4) for the 24 UMP stars in common. The agreement is generally excellent. Specifically, defining Δ as the difference between our values and those of Sestito et al. (2019) normalized by our values, then for the 24 stars we find median Δ values of 0.04, 0.03, 0.03, and 0.08 for D_apo, D_peri, J_ϕ, and E, respectively, noting that for the energy comparison we have taken into account the different solar energy used here to that in the Sestito et al. (2019) study.

4 RESULTS

The physical properties and the computed orbital parameters of the first 10 stars are listed in Tables 1 and 2 while the complete tables are available with the online Supporting Information. For Table 1 the columns are, respectively, an index number, the Gaia DR2 and SkyMapper or other IDs, the on-sky location in degrees, the parallax and its uncertainty from Gaia DR2, the adopted distance and its uncertainty, a flag indicating whether the distance is from the RGB isochrones (value=0), or from Bailer-Jones et al. (2018) (value=1), the proper motions from Gaia DR2 and their uncertainties, the heliocentric radial velocity and its uncertainty, |$\log \, T_{\rm eff}$| and its uncertainty, the abundance [Fe/H], the reddening, and the data set from which the star originates. Similarly for Table 2, the columns are the index number (as for Table 1), the eccentricity and its uncertainty, the apo- and peri-galactic distances, the maximum deviation from the Galactic plane, the actions |$(J_{\rm R},\, J_{\rm \phi },\, J_{\rm Z})$|⁠, the energy, and the U, V, and W velocity components in the LSR frame.

Fig. 5 shows the inferred orbital parameters for all the stars with D_apo ≤ 250 kpc, with each star colour-coded according to its metallicity, as shown in the right colour bar. In particular, panels (a) and (b) show the vertical action |$(J_{\rm Z}\, [\rm kpc\, \rm km\, s^{-1}])$|⁠, indicative of the vertical excursion of the star, and the orbital energy |$(E\, \rm [km^2\, s^{-2}])$|⁠,¹² as a function of the azimuthal action |$(J_{\rm \phi }\, \rm [\rm kpc\, \rm km\, s^{-1}])$|⁠.¹³ The quantities have been normalized by the solar values computed for the McMillan2017 potential employed here: |$J_{\rm \phi , \odot }=2014.24\, \rm kpc\, \rm km\, s^{-1}$|⁠, |$J_{\rm Z, \odot }=0.302\, \rm \rm kpc\, \rm km\, s^{-1}$|⁠, and |$E_{\rm \odot }=-153507.15\, \rm km^2\, s^{-2}$|⁠. We note that if we adopt the MWPotential2014 employed by Sestito et al. (2019), and include the increased dark matter halo mass, we obtain solar values similar to those in that work. Retrograde orbits are characterized by a negative value of J_ϕ, while prograde orbits have a positive J_ϕ. We find that overall |$\sim 42{{\ \rm per\ cent}}\, (185/445)$| of our stars with D_apo ≤ 250 kpc exhibit retrograde orbits, and note that the selection of the stars for inclusion in our sample should not have any bias as regards prograde or retrograde orbits.

Regarding the uncertainties in the derived orbital quantities, these are listed for each individual star in Table 2, but as examples, we find that for stars with D_apo ≤ 20 kpc, the median errors in D_apo, D_peri, |Z_max|, and e are 0.9 kpc, 0.6 kpc, 1.0 kpc and 0.08, respectively. These increase to 11.6 kpc, 1.7 kpc, 7.7 kpc, and 0.08 for stars with 20 ≤ D_apo ≤ 50 kpc, respectively, and to 56.4 kpc, 2.7 kpc, 52.0 kpc, and 0.09 for stars with 50 ≤ D_apo ≤ 250 kpc.

In panels (c) and (d), we show the maximum height |$(|Z_{\rm max}|\, \rm [kpc])$| and eccentricity e as function of the apogalatic distance |$(D_{\rm apo}\, \rm [kpc])$|⁠. A preliminary inspection of panels (a) and (c) reveals that, despite the low metallicity of the stars in the sample, we detect a significant number of stars with small vertical excursion, in agreement with Sestito et al. (2019, 2020b) and Di Matteo et al. (2020). In particular, if we follow Sestito et al. (2020b) and adopt J_Z/J_Zsun < 1.25 × 10³, shown as the dotted horizontal line in panel (a), to characterize orbits that are confined to the disc, then ∼50 per cent of our sample meets this definition. Similarly, if we follow Sestito et al. (2019) in using |Z_max| = |$3\, \rm kpc$| (horizontal dashed–dotted line in panel c) of Fig. 5) to discriminate between ‘disc-like’ and ‘halo-like’ orbits, we find that 102, or ∼21 per cent, of the stars in our sample meet this criterion, that is, have orbits that do not deviate far from the Galactic plane.

Further, panel (d) suggests that, while stars with |$D_{\rm apo}\lesssim 25\, \rm kpc$| have an approximately uniform distribution in eccentricity, highly eccentric (e ≳ 0.5) orbits are favoured for stars with |$D_{\rm apo}\gtrsim 25\, \rm kpc$|⁠, while panels (c) and (f) show that there is an apparent dearth of stars with low values of |Z_max| beyond |$D_{\rm apo}\, \approx$| 30 kpc. These apparent effects are most probably a consequence of the criteria adopted to select SkyMapper EMP candidates, as stars with low |Z_max| and large D_apo are not likely to meet the apparent magnitude cut that underlies the sample (g_skymapper < 16 for the HiRes stars and g_skymapper < 13.7 for the LowRes stars). In particular, the bottom middle and bottom right panels of Fig. 2 show that stars with |$|Z|\, \le$| 3 kpc and Galactocentric distances beyond 10–15 kpc are rare in our sample.

The two bottom left panels again show the eccentricity versus the apogalactic distance, but separately for stars with |Z_max| in excess of 3 kpc (panel e) and those with |Z_max| less than this value (panel f). Similarly, the two bottom right panels also show eccentricity versus the apogalactic distance but this time the sample is split by metallicity: stars with |$\rm [Fe/H]\gt -3$| are shown in panel (g), while the more metal-poor stars are shown in panel (h). The similarity of panels (g) and (h) show that there is no obvious dependence of the kinematics on metallicity, at least for this sample of metal-poor stars.

To more clearly illustrate this point, we show in Fig. 6 a plot of [Fe/H] against e, the orbital eccentricity. Diagrams of this nature have long played an important role in discussions of the formation of the Galaxy. For example, in their classic paper, Eggen, Lynden-Bell & Sandage (1962) argued on the basis of an apparent correlation between ultraviolet excess (an indicator of [Fe/H]) and orbital eccentricity, that the proto-Galaxy collapsed rapidly to a planar structure with a time-scale of only a few × 10⁸ yr. Specifically, in their sample of stars, those with [Fe/H] less than –1.5, approximately, all had e ≥ 0.6 (Eggen et al. 1962). Norris et al. (1985) challenged the rapid collapse interpretation arguing that the lack of low-e metal-poor stars was a result of a kinematic bias in the selection of the Eggen et al. (1962) sample. Instead, using a sample selected without any kinematic bias, Norris et al. (1985) showed that metal-poor stars with relatively low orbital eccentricities exist, a population they identified as a metal-weak component of the thick disc. The Norris et al. (1985) result was confirmed and strengthened by Beers et al. (2014, see their fig. 10) who showed that for stars with [Fe/H] ≤ –1.5 there is no correlation between orbital eccentricity and metallicity: stars can be found with e values between ∼0.1 and 1. Our results in Fig. 6 extend the lack of any correlation to substantially lower metallicities than those in Beers et al. (2014), where there were only a few stars at or below [Fe/H] = –2.5 and none below –3.0 dex. We discuss the implications of the existence of EMP stars with low eccentricities (and low |Z_max|) in Section 5.2. Fig. 6 also shows the location of candidate members of the Gaia Sausage and Gaia Sequoia accretion events. The identification and properties of these stars are discussed in detail in Section 5.1.

Finally, as noted above, we find that 30 stars from the full sample have apparent D_apo values larger than 250 kpc, that is, larger than the virial radius of the MW. The majority of these stars possess energies that are consistent with, or larger than, zero and they likely have unbound orbits. These stars will be discussed in more detail in Section 5.3, but we note again that they are not plotted in the panels of Fig. 5 or  6.

5 DISCUSSION

In the following subsections, we discuss in detail the results for the 475 very metal-poor stars analysed. We will focus specifically on three key aspects. The first is the relation between the stars in our sample and the recently described remnants of the postulated Gaia Sequoia and Gaia Sausage accretion events (Belokurov et al. 2018; Myeong et al. 2019). For the sake of this analysis, following the hypothesis of Belokurov et al. (2018) and Myeong et al. (2019), we assume that these accretion events are distinct, but see Helmi et al. (2018) for an alternative view, particularly of Gaia Enceladus as a single ancient major merger event. The purpose of our work, however, is not to discern between the scenarios proposed to explain these structures in the Galactic halo, but rather to investigate their very low-metallicity content. The second key point is the analysis of low-metallicity stars with disc-like orbital properties that likely have a fundamental role in contributing to the understanding of the formation and evolution the MW’s disc. Finally, we discuss the properties and potential origin of the stars in our sample that are either loosely bound or not bound to the Galaxy.

5.1 Gaia Sausage and Gaia Sequoia candidate members

The exquisite data provided by Gaia DR2 has recently revealed the trace of at least two early major accretion events in the history of our Galaxy, referred to as Gaia Sausage and Gaia Sequoia (Helmi et al. 2018; Belokurov et al. 2018; Mackereth et al. 2019; Myeong et al. 2019; Koppelman et al. 2019). These discoveries are a direct consequence of the development of computational techniques and resources capable of processing very large data sets.

Here we exploit the action-space classification provided in Myeong et al. (2019, their fig. 9) to identify possible members of these accretion features within our sample of low-metallicity stars. Monty et al. (2020) have adopted a similar approach finding possible members of these systems with metallicities as low as |$\rm [Fe/H]=-3.6$| dex. The number, abundances, and abundance ratios of these stars could provide important information on the early evolution of the progenitors of the two accretion events.

The top left panel of Fig. 7 shows the action map (J_Z − J_R)/J_tot versus J_ϕ/J_tot with J_tot being the sum of the absolute value of the three actions (J_tot = J_R + J_Z + |J_ϕ|). Following the classification in Myeong et al. (2019), we highlight the loci of the Sequoia and Sausage accretion events with red and blue rectangles, respectively.

We find that out of the 475 analysed stars, 16 stars are kinematically coincident with the Sausage accretion event, while 40 stars are candidate Sequoia members. As expected from their definition and the action map (Helmi et al. 2018; Belokurov et al. 2018; Myeong et al. 2019; Yuan et al. 2020), the latter are characterized by mildly eccentric (e ∼ 0.5) retrograde orbits, while the former have highly eccentric orbits (e ∼ 0.9). Appendix  B shows some typical orbits for stars identified as possible Gaia Sequoia and Gaia Sausage members.

We remind the reader that our membership identification follows the criteria introduced in Myeong et al. (2019), and is thus entirely based on the dynamics through the use of the action map. We stress that this approach does not allow for any ‘background’ population that may be present in these regions of the action map. Consequently, we cannot straightforwardly assume that all the stars in our sample that are dynamically coincident with the Sequoia/Sausage accretion events actually belong to such remnants. In Appendix  A, we have attempted to perform a more accurate analysis through the use of a clustering algorithm approach. Briefly, the clustering analysis of our very metal-poor sample does provide independent evidence for the existence of groupings consistent with the Sequoia (group 6) and Sausage (group 8) dynamical definitions, though there are also indications that our Sequoia and Sausage samples, as defined in Fig. 7, are potentially contaminated by a ‘background’ population that might be as much as ∼ 50 per cent and ∼35 per cent, respectively. These background estimates are determined by exploiting the clustering analysis groupings discussed in Appendix  A, and the numbers of stars within the Gaia Sequoia and Gaia Sausage loci.

Panels (b)–(e) of Fig. 7 show a detailed analysis of stars identified as candidate Sequoia, shown in red, and Sausage members, shown in blue. In panel (c) we note that, by construction, Sausage stars are characterized by more radial orbits, although at low J_R, some candidate Sequoia stars seem to share the similar values of J_R as Sausage stars. The Toomre diagram in panel (d) shows that both groups are consistent with halo dynamics, and again we note that there is some degree of overlap between the two groups of stars. As regards panel (b), which shows the energy versus azimuthal action, Sausage candidates show the distinctive vertical distribution, indicative of almost null azimuthal angular momentum, while Sequoia stars are clearly highly retrograde, as expected. Comparing panel (b) with Koppelman et al. (2019, their fig. 2) we note that our accreted candidates span a wider range in energy. However, we note that the definition of Sequoia and Sausage parameters differs from work to work. Indeed, Yuan et al. (2020) identify Sausage members that lie well outside the selection box of Myeong et al. (2019) and the energy range of Koppelman et al. (2019). For the sake of our analysis, we choose to be consistent with the Myeong et al. (2019) classification, although we stress again that a number of the candidates may not actually belong to the remnants of the accretion events.

As regards abundances, we find that the most metal-poor star in our sample that is a candidate member of Sequoia (SMSS J081112.13−054237.7) has a metallicity of |$\rm [Fe/H]=-3.74$|⁠, while the most metal-poor Sausage candidate (SMSS J172604.29−590656.1) has |$\rm [Fe/H]=-3.31$| dex. Both stars come from the HiRes sample so that the abundance uncertainty is of order 0.1 (excluding any systematic uncertainties such as those arising from the neglect of 3D/NLTE effects). These values are quite consistent with the results of Monty et al. (2020). In that work, which uses dwarf stars, the lowest metallicity star plausibly associated with Sequoia, G082–023, has |$\rm [Fe/H]=-3.59 \pm 0.10$| while the most metal-poor star plausibly associated with Sausage, G064–012, has |$\rm [Fe/H]=-3.55 \pm 0.10$| (Monty et al. 2020).

Finally, for the stars in the HiRes, Jacobson+15, and Marino+19 samples, we are able to investigate the chemical patterns of the likely accreted stars. Panel (f) of Fig. 7 shows [|$\rm \alpha$|/Fe] versus [Fe/H] for the 218 stars for which [|$\rm \alpha$|/Fe] values are available. Specifically, [|$\rm \alpha$|/Fe] is computed as the unweighted mean of [Mg/Fe], [Ca/Fe], [TiI/Fe], and [TiII/Fe] where available.¹⁴ The star SMSS J160540.18–144323.1 (Nordlander et al. 2019) has been arbitrarily plotted at a metallicity of [Fe/H]=–4.3 since otherwise it would be the only star with [Fe/H]<–5 in the panel. A visual inspection reveals that, with the single exception of SMSS J111201.72−221207.7, all the Sequoia and Sausage candidates are α-enhanced, and no other trend is evident. Specifically, in our sample of very metal-poor Sequoia and Sausage candidates, we see no evidence for a ‘knee’ in the ([|$\rm \alpha$|/Fe], [Fe/H]) relation. The metallicity of the knee marks the abundance where [|$\rm \alpha$|/Fe] begins to decrease with increasing [Fe/H] as the nucleosynthetic contributions from supernovae (SNe) Ia become increasingly important. Our result, however, is not inconsistent with the results of Matsuno, Aoki & Suda (2019) and Monty et al. (2020) who find evidence of the presence of a knee at higher abundances than any of the Sequoia and Sausage candidates plotted in Fig. 7. For example, Monty et al. (2020) indicate that the knee in Sequoia is at [Fe/H] ≈ –2 while that for Gaia Sausage is at [Fe/H] ≈ –1.6, values significantly more metal-rich than any of the candidates in Fig. 7.

5.2 A very metal-weak component in the thick disc?

It has recently been shown (Sestito et al. 2019, 2020b; Di Matteo et al. 2020; Venn et al. 2020) that, in contrast to the commonly accepted view, a significant fraction (⁠|$\sim 20{{\ \rm per\ cent}}$|⁠) of very low-metallicity stars resides in the MW disc rather than in the halo. Here, we find a similar result: for the stars in our sample 102 out of 475 (⁠|$\sim 21\, {{\ \rm per\ cent}}$|⁠) exhibit disc-like dynamics in having orbits that are confined to within 3 kpc of the plane of the MW.

The straightforward conclusion would be to propose these stars as supporting the existence of an extension to yet lower metallicities of the proposed metal-weak thick disc (e.g. Chiba & Beers 2000). However, a detailed analysis is required to discern the origin of these stars. Specifically, it is of key importance to understand if they are indeed disc-like stars, or if they are, for example, halo stars whose orbital plane happens to lie in the MW disc.

We therefore explore a number of orbital parameters to shed light on the nature of these stars. Specifically, Sales et al. (2009, and references therein) have shown how the orbital eccentricity can be used to probe the formation scenario of the MW thick disc. Furthermore, as mentioned in the previous sections, the actions, and in particular the azimuthal action, provide important clues for the origin of a star. Here, we couple these two orbital properties to disentangle the origin of the very low-metallicity stars residing in the MW disc.

As noted above, we classify as ‘disc’ stars those with |$|Z_{\rm max}|\le 3\, \rm kpc$|⁠. This choice is made on the basis of the following considerations. First, Li & Zhao (2017) find that the (exponential) scale height of the MW thick disc is |$z_{\rm 0}=0.9\pm 0.1\, \rm kpc$|⁠; it then follows that the vast majority of the thick disc population should be found within ∼3 scale heights, that is, within |Z_max| = 3 kpc.

Second, Fig. 8 shows the eccentricity distribution for different values of the maximum vertical excursion, from |$2$| to |$8\, \rm kpc$|⁠. For each panel, the yellow histogram (designated as ‘disc’) shows the distribution of stars within |$N\, \rm kpc$|⁠, with N = {2, 3, 4, 5, 8}, while the grey shaded histogram (designated as ‘halo’) represents stars with |$|Z_{\rm max}|\gt N\, \rm kpc$|⁠. As is evident from the figure, for heights above the plane exceeding 3 kpc, that is, panels (c)–(e), the eccentricity distributions for the stars above and below the cut-off height become increasingly similar. On the other hand, for a cut-off value of |$|Z_{\rm max}|=3\, \rm kpc$|⁠, an apparent difference is present in the sense the e-distribution for the low |Z_max| stars has a possible excess of intermediate eccentricity stars together with a possible narrow surfeit of stars with |$e\, \approx$| 0.85, which is also evident in panel (a). However, application of both Kolmogorov–Smirnoff and Anderson–Darling tests (see e.g. Scholz & Stephens 1987) to compare the ‘disc’ and ‘halo’ distributions in panel (b) revealed that the apparent differences are not statistically significant. None the less, we adopt |$|Z_{\rm max}| = 3\, \rm kpc$| as the value of |Z_max| to discriminate between predominantly disc and predominantly halo populations. For completeness, we also note that if we choose |Z_max| cut-off values of 2.5 or 4 |$\rm kpc$| and repeat the analysis discussed below, the outcomes are essentially unaltered.

Finally, the third reason for adopting a value of |$|Z_{\rm max}| = 3\, \rm kpc$| for the disc-like population is that it is consistent with Sestito et al. (2019, 2020b), allowing our results to be directly compared with theirs.

Looking again in detail at panel (b) in Fig. 8, we can see that the eccentricity distribution of the disc-like stars hints at the presence of two main groups. The first group has a relatively broad distribution peaking at e ≈ 0.55 while the second population has a narrower distribution centred at e ≈ 0.85. This interpretation is confirmed by the application of Gaussian Mixture Modelling to the e-distribution for the stars with |$|Z_{\rm max}|\le 3\, \rm kpc$|⁠, a process that does not require any choice as regards histogram bin size. The best fit is for two Gaussians, one centred at e = 0.52 containing 80 per cent of the population and with a standard deviation of 0.14. The second Gaussian is centred at e = 0.89 with a narrow σ of 0.03. The two Gaussians are overplotted with blue thick lines in panel (b) of Fig. 8. We shall refer to these two groups as the ‘low-eccentricity’ and ‘high-eccentricity’ populations, respectively, and adopt e = 0.75 as the eccentricity to separate them.

We now employ J_ϕ to identify the motion of the disc stars as either prograde or retrograde. This is shown in Fig. 9, where we mark retrograde high-e and low-e population stars with dark-blue and dark-red points, respectively, while prograde high-e and low-e group stars are indicated by azure and orange circles.¹⁵ Overall we find 72 low-e stars (53 prograde and 19 retrograde) and 30 high-e stars (15 prograde and 15 retrograde).

We then analysed the orbital parameters of each of these subgroups of disc-star candidates. The results are shown in Fig. 10. The distributions shown in panels (a)–(j) are kernel density distributions computed by adopting a Gaussian kernel with a fixed value of 0.4 for the bandwidth scaling parameter, while the top three panels show the action map, the Toomre diagram, and the E versus J_ϕ plot, respectively. The disc candidates are marked with different colours, as defined in Fig. 9.

5.2.1 Low-eccentricity stars

In panel (a) of Fig. 10, it is interesting to see that the eccentricity distributions of the prograde and retrograde high-e groups are almost identical, while, conversely, there are hints of a difference in the corresponding distributions for the low-e stars.

Specifically, the prograde low-e stars (orange line) have a broader distribution while the retrograde low-e stars (red line) have a narrower distribution peaking at e ∼ 0.5–0.6. Together these e-distributions are quite consistent with the theoretical results in Sales et al. (2009), particularly as regards the e-distributions in the top-left panel of their fig. 3 (Sales et al. 2009). In that context the retrograde low-e stars can be interpreted as an accreted population, while the prograde low-e stars are likely ‘in situ’, that is, born within the thick disc of the Galaxy.

To support this interpretation we consider again the action map, here shown in the upper left of Fig. 10 with the |$|Z_{\rm max}| \le 3\, \rm kpc$| stars identified. It is evident from this panel that the majority of the retrograde low-e stars fall within the locus defining the Gaia Sequoia accretion event, consistent with these stars having an accretion origin. Comparing the corresponding panels in Figs 7 and 10 for the Toomre diagram and the Energy (E) against azimuthal action (J_ϕ) diagram, respectively, confirms the connection.

Regarding the prograde low-e stars, a substantial number of these fall in the region of the Toomre diagram usually restricted to disc stars; their rotation velocities lag that of the Sun by relatively small amounts, less than 100 |$\rm km\, s^{-1}$| in some cases. We therefore conclude that the prograde low-e stars define a very metal-weak component to the Galaxy’s thick disc. This conclusion is supported by the eccentricity distribution of the stars, which agrees well with the eccentricity distributions for (more metal-rich) thick-disc stars shown in Li & Zhao (2017, their figs 12 and 14).

These 53 low-e, low Z_max prograde stars represent ∼ 11 per cent of our total sample. Of these 53, six are included in the high-dispersion data sets and the [α/Fe] versus [Fe/H] for these stars is shown in Fig. 11. Four of the 19 low-e, low Z_max retrograde stars are also included on the plot along with the remainder of the stars in the high-dispersion data sets. For completeness as regards the [Fe/H] distributions of the samples, we also show in the upper part of the figure the [Fe/H] values for the remainder of the prograde and retrograde low-e, low Z_max samples. Detailed abundance information, such as [α/Fe], is not available for these stars that arise from the LowRes sample. Further, in order to avoid any potential systematic effects, we have chosen not to plot the [α/Fe] values for the three low-e, low Z_maxSestito+19 sample stars in Fig. 11. The stars are BD+44 493 ([Fe/H] = –4.30), 2MASS J18082002−5104378 ([Fe/H] = –4.07) and LAMOST J125346.09+075343.1 ([Fe/H] = –4.02) and all three have prograde orbits. The [Fe/H] values for these stars are taken from table 1 of Sestito et al. (2019).

It is evident from the figure that, though the sample is very limited, the four retrograde low-e stars show no obvious difference in location in this plane from the remainder of the (halo dominated) sample with high-dispersion abundance analyses. The location of the six prograde low-e stars, however, is intriguing despite the small numbers. It appears that the mean [α/Fe] for these stars is lower than that for the full sample by perhaps 0.15, and two of the stars, namely SMSS J230525.31−213807.0 which has [Fe/H] = –3.26 and which is also known as HE 2302−2154a, and SMSS J232121.57−160505.4 ([Fe/H] = –2.87, HE 2318−1621), are among the small number (of order a dozen in total) that have [α/Fe] < 0.1 in the high-dispersion samples. Such α-poor stars, also referred to as ‘Fe-enhanced’ stars (e.g. Yong et al. 2013; Jacobson et al. 2015), may reflect formation from gas enriched in SNe Ia nucleosynthetic products that, if valid, may have implications for the epoch at which these stars settled into, or formed in, the thick disc (see also Sestito et al. 2019, 2020b; Di Matteo et al. 2020). However, we consider further discussion of the element abundance distributions in these stars beyond the scope of the present paper.

The overall [Fe/H] distribution of the retrograde and prograde stars as inferred from Fig. 11 is similar to that for the full sample. The star with the lowest combination of metallicity and eccentricity in our sample of prograde, low Z_max stars is SMSS J190836.24–401623.5 that has [Fe/H] = –3.29 ± 0.10 and for which we find e = 0.29 and a high V velocity of –21 km s⁻¹. We also note that the orbital parameters derived here for the UMP star 2MASS J1808002–5104378, which is included in our Sestito+19 data set, are very similar to those found in Sestito et al. (2019). Specifically we find for this star e = 0.13 and V = –29 km s⁻¹, while Sestito et al. (2019) list values of 0.09 and –45 km s⁻¹, respectively. Moreover, as noted by Sestito et al. (2019), the ‘Caffau-star’ (Caffau et al. 2011), which is an apparently carbon normal (i.e. [C/Fe] < 0.7) dwarf (and therefore not included in our sample) with [Fe/H] ≈ –5.0, and which is the star with the lowest total metal abundance known to date, is a further example of an ultra-low-metallicity star with a disc-like orbit. Sestito et al. (2019) determine that this star has a prograde orbit with e = 0.12 that is confined to the Galactic plane, and which has a high V velocity of –24 km s⁻¹. We agree with the suggestion of Sestito et al. (2019) that these stars may have formed in a gas-rich ‘building-blocks’ of the proto-MW disc. The origin of these stars is also discussed within the context of the theoretical simulations presented in Sestito et al. (2020a). The simulations reveal the ubiquitous presence of populations of low-metallicity stars confined to the disc plane. In particular, the simulations show that the prograde planar population is accreted during the assembly phase of the disc, consistent with our interpretation and that of Sestito et al. (2019).

Panels (b)–(h) of Fig. 10 show the kernel density distributions of the other orbital parameters. These distributions do not exhibit clear differences between prograde and retrograde low-e stars, with the possible exception of D_apo, in panel (e), and D_peri, in panel (f), which mimic the differences in the eccentricity distribution evident in panel (a).

We conclude that the low Z_max prograde and retrograde low-e stars likely have different origins, with the former possibly being formed in situ in the Galaxy’s thick disc, while the latter are likely accreted from disrupted MW satellites. Whether these latter stars belong to the Sequoia main remnant, or its higher energy tails (Thamnos 1 and 2, Koppelman et al. 2019) is beyond the scope of the present work, although, as panel (a) of Fig. 10 shows, most fall within the region defining Sequoia stars.

5.2.2 High-eccentricity stars

The interpretation of the 30 low Z_max, high-e stars in our sample is less straightforward, since, as the panels of Fig. 10 reveal, nearly all their orbital parameters show similar distributions for the prograde and retrograde stars, with the only exception being the azimuthal action (by construction). We note first that the only mechanism able to explain, at least qualitatively, the occurrence of disc stars with high-eccentricity orbits, is the heating mechanism discussed in Sales et al. (2009, e.g. fig. 3). Some of the stars in this subsample could therefore be disc stars heated by accretion events. Alternatively, some could simply be halo stars that happen to have orbital planes that lie close to the MW disc plane.

However, even though these high-e stars do not satisfy the Gaia Sausage (Myeong et al. 2019) membership criteria, we find that many qualitatively share its typical orbital properties: high-eccentricity, low or no angular momentum, small Galactic pericentres, Galactic apocentres as great as ∼20–|$25\, \rm kpc$| and strong radial motions. We therefore argue that at least some stars could be associated with the Sausage accretion event. In support of this conjecture, we note that they are consistent with the Yuan et al. (2020) classification for Gaia Sausage stars, both in the action map locus, and in the energy regime.

Myeong et al. (2019) suggest that the Sausage accretion event was an almost head-on collision with the MW. Such an event generates orbits with small perigalacticons that are strongly radial, eccentric, and with roughly equal numbers of prograde and retrograde stars. These are the properties that we see for our low Z_max, high-e stars and it therefore seems reasonable to conclude that many of our low Z_max, high-e stars have their origin in the Sausage accretion event.

For seven of the 15 prograde high-e, low Z_max stars, detailed abundances are available from our high-dispersion data sets. The [α/Fe] abundance ratios are shown as a function of [Fe/H] in Fig. 11. The corresponding data for 5 of the 15 retrograde stars are also shown in the figure. The [Fe/H] values for the remainder of the stars in these groups are shown across the top of the plot. The numbers of stars with high-dispersion analyses in both high-e, low Z_max subsamples are small, but there does not appear to be any obvious difference between them and the distribution of the full sample.

5.3 The candidate unbound stars

As introduced in Section 4, we find 30 stars with potentially unbound orbits, defined by having an apparent |$D_{\rm apo} \ge 250\, {\rm kpc}$|⁠.¹⁶ One star, HE 0020–1741, comes from our subset of the Sestito+19 sample. Sestito et al. (2019) list |$D_{\rm apo} \approx 296\, \rm kpc$| for this star, the largest value in their determinations, while in the Mackereth & Bovy (2018) catalogue¹⁷ the star is classified as unbound in accord with our result. The remaining 29 stars are from the SkyMapper samples.

To test the sensitivity of the results to the adopted potential we investigated the orbits of these stars using a different choice of the potential, namely the galpyMWPotential2014. The calculations reveal that all 30 stars again have |$D_{\rm apo} \gt 250\, \rm kpc$|⁠. In addition, for this choice of potential, we find that the star SMSS J044419.01–111851.2 may also be unbound; it is on a loosely bound orbit in the McMillan2017 potential with an apogalacticon distance of |$\sim 200\, \rm kpc$|⁠.

In order to shed light on the nature of these stars we have investigated the distributions of the 500 random realizations of the orbits, together with the corresponding orbital parameters, such as the apparent apo/peri-galacticon distance ratio (D_apo/D_peri), the energy (E) and the azimuthal action (J_ϕ). Panels (a)–(d) of Fig. 12 show the information for the 30 potentially unbound stars.

Since the direct orbit integration for the observed positions and velocities resulted in unbound orbits (apparent |$D_{\rm apo}\gt 250\, \rm kpc$|⁠), it would be incorrect to adopt the uncertainties as defined in Section 3.2 for all parameters. Specifically, for these stars we opt not to give uncertainties for the apparent |$D_{\rm apo},\, e,\, Z_{\rm max}$|⁠, and J_R values since the medians for the 500 realizations and the ‘observed’ values, that is, the values from the observed properties, differ significantly. On the other hand, for the remaining orbital parameters, namely |$D_{\rm peri},\, E,\, J_{\rm \phi },\, J_{\rm Z},$| and the |$U,\, V,\, W$| velocities, the median values are consistent with the observed ones, and therefore we compute the uncertainties in these quantities as before, from the difference between the 16th and the 84th percentile of the PDFs.

Panel (a) of Fig. 12 shows the fraction of unbound realizations for each star. Each bar is colour-coded according to its percentage, from a minimum of 42.6 per cent for SMSS J095211.09–185713.7 (number 3 in the identification panel in the Fig. 12) to a maximum of 100 per cent for SMSS J090247.41–122755.1 (number 26). White- and reddish colour tones indicate stars with a fraction of unbound orbits greater than 70 per cent, which we take as a conservative value to identify likely unbound stars. Blueish colours represent stars with a lower unbound fraction. A visual inspection of panel (b) reveals that nearly all stars with |$N_{\rm unbound}\gt 70{{\ \rm per\ cent}}$| exhibit positive energies, thus confirming that they are likely to be escaping from the Galaxy. Overall, we find that 17 stars have |$N_{\rm unbound}\gt 70{{\ \rm per\ cent}}$|⁠, and of these 15 have E > 0 (the two stars with |$N_{\rm unbound}\gt 70{{\ \rm per\ cent}}$| but E ≤ 0 are numbers 22 and 23). Panel (b) also shows that four stars (numbers 1, 8, 11, and 13) have a positive energy, but with N_unbound slightly below 70 per cent. We consider these stars as also likely unbound, bringing the total number of candidate unbound stars to 21, or 4.4 per cent of the total sample. We note in particular that aside from HE 0020–1741 (number 10), three further stars in our set of 21 unbound candidates are also classified as unbound in the Mackereth & Bovy (2018) catalogue. These are the bright r-process element enhanced star SMSS J203843.18–002332.8 (number 11, RAVE J203843.2–002333, Placco et al. 2017), together with SMSS J183246.48–343434.3 (number 30) and SMSS J202221.56–121443.9 (number 13). On the other hand, one of our stars, SMSS J103622.54−713010.3 (number 12), has a bound orbit in the Mackereth & Bovy (2018) catalogue. There are no stars in common with the list of 20 ‘clean’ high-velocity star candidates with unbound probability exceeding 70 per cent in Marchetti, Rossi & Brown (2019).

The remaining nine stars have negative (bound) energies, although the values are consistent with zero within their uncertainties. Their classification is thus uncertain as they could be unbound or on loosely bound orbits. None are found in the Mackereth & Bovy (2018) catalogue.

We now speculate as to the origin of these stars, proposing three possible physical mechanisms that could provide each star with sufficient energy to escape the Galaxy.

• A star in a close binary can be expelled from the GC via an interaction with the central black hole. The clearest example of this process is the star S5-HVS1 discussed in Koposov et al. (2020).

• A star can acquire high velocity (of order of the binary’s pre-SNe explosion orbital velocity) from being in a binary when the companion explodes as an SN (e.g. Eldridge, Langer & Tout 2011).

• A star can acquire energy as part of a gravitational interaction involving the merger of a dwarf galaxy with the MW (e.g. Abadi, Navarro & Steinmetz 2009).

For the first mechanism to happen, the star has to have an origin close to the GC, which means that its D_peri should be near zero. However, as panel (c) of Fig. 12 shows, only one star (SMSS J115906.91–261050.6, number 1, D_peri = 0.26 kpc) has a perigalacticon distance within 1 kpc of the GC, while the remainder of the candidates have |$D_{\rm peri} \gt 1\, {\rm kpc}$|⁠. This suggests that the first possibility is unlikely, particularly when it is recognized that all the stars in the sample are giants and therefore unlikely to be in a sufficiently compact binary.

The second mechanism also seems unlikely because the unbound stars are all giants with, as a consequence, relatively large stellar radii. As a result, the separation between the components of any pre-SNe binary containing the star is unlikely to be sufficiently small that the orbital velocity, which underlies the ‘kick velocity’ provided when the companion becomes a SNe, would be sufficiently high that the liberated star is no longer bound to the Galaxy. Furthermore, at least for those stars where high-dispersion spectra are available, there is no evidence of any ‘pollution’ from the SNe event.

This leaves us with the third possible origin, which is plausible given that it is generally accepted that the formation of the Galactic halo is driven by the accretion and tidal disruption of dwarf galaxies (e.g. the recent discovery of remnants of accretion events: Helmi et al. 2018; Belokurov et al. 2018; Koppelman et al. 2019; Myeong et al. 2019). Specifically, we postulate that our set of unbound stars originated in the outskirts of dwarf galaxies that were accreted by the MW, gaining energy from the gravitational interaction that resulted in the disruption of the dwarfs. Given the relatively low metallicities of the unbound stars, which range from –4 (or less) to –2 in [Fe/H] (see the lowermost panel of Fig. 12), we speculate that the disrupted systems were relatively low-mass, low-metallicity systems.In this context, we note that 16 out of 21 stars have prograde orbits, while five have negative J_ϕ and thus a retrograde orbit. This likely indicates that multiple accretion events may be involved. However, it is necessary to keep in mind that the time-scale for an unbound star to reach the virial radius from the inner regions of the Galaxy is ∼ 1 Gyr. Consequently, the gravitational interactions that generated the unbound stars in our sample likely occurred relatively recently, which may argue against the proposed ‘origin in accretion events’ scenario. Detailed evaluation of the orbits of the unbound stars, individually and collectively, is required to assess the situation and to investigate their origins(s). The metallicities of other candidate unbound stars, such as those in Marchetti et al. (2019) will also provide important input (e.g. Hawkins & Wyse 2018).

We note also that there is a fourth possibility, that uncertainties in the analysis lead to incorrect orbital parameters. For example, if the distance to the star used in calculating the orbit were overestimated, this could result in unbound or nearly bound status. This is the likely explanation for the discrepancy concerning the star SMSS J103622.54–713010.3 (unbound here, bound in the Mackereth & Bovy (2018) catalogue) as our adopted distance is more than a factor of two larger than the Bailer-Jones et al. (2018) distance. The availability of improved parallaxes and stellar parameters from the forthcoming Gaia EDR3 and DR3 releases will help alleviate these discrepancies. As another possibility, we note that the total mass of the Galaxy may in fact be larger than that used in our modelling. If this is the case then, although on high-energy orbits, the stars would remain bound (e.g. Monari et al. 2018; Fritz et al. 2020).

Panel (e) of Fig. 12 shows the [α/Fe] versus [Fe/H] relation for 26 of the 30 candidate unbound stars that are in our high-dispersion samples, together with the values for the remainder of the stars from the high-dispersion samples. The estimates of [Fe/H] for the remaining four stars, from the LowRes data set, are shown at the top of the panel. It is evident from this panel that candidate unbound stars are not distinguished from the full sample as regards the overall metallicity distribution or the [α/Fe] distribution.

6 SUMMARY AND CONCLUSIONS

In this work, we have analysed a sample of 475 very metal-poor giant stars, most of which have originated from the SkyMapper search for the most metal-poor stars in our Galaxy (Da Costa et al. 2019). The data set covers a metallicity range of almost five dex (− 6.5 < [Fe/H] ≤ −2), and together with the relatively large number of stars, makes it ideal to investigate the kinematics, and ultimately the origin, of these very rare and important objects together with the implications for the formation of the MW.

We first exploited the action map for our sample together with the classification criteria of Myeong et al. (2019) to identify candidate members of the Gaia Sausage and Gaia Sequoia accretion events. We find 16 stars dynamically consistent with Gaia Sausage and 40 with Gaia Sequoia. While we cannot be certain all candidates are in fact associated with these entities, the lowest metallicities ([Fe/H] = –3.31 for Gaia Sausage and –3.74 for Gaia Sequoia) are quite consistent with the findings of Monty et al. (2020). With a single exception, all our candidate Gaia Sausage and Gaia Sequoia stars for which we have high-dispersion spectra are α-rich, similar to the general halo population. This is again consistent with the results of Monty et al. (2020).

The recent work of Sestito et al. (2019, 2020b), Di Matteo et al. (2020), and Venn et al. (2020) has revealed an unexpected significant population of very low-metallicity stars residing in the plane of the Galaxy. We find a similar result in that |$\sim 21\, {{\ \rm per\ cent}}$| of the stars in our sample have orbits that remain confined to within 3 |$\rm kpc$| of the Galactic plane. Moreover, these stars show a different eccentricity distribution compared to the stars with larger |Z_max| values, pointing towards a different origin and/or evolution compared to the (halo dominated) bulk of the sample.

Our detailed analysis of these low |Z_max| stars reveals four subpopulations as regards orbit eccentricity and prograde or retrograde motion. Of particular interest are the stars with relatively low eccentricities (e < 0.75, median ≈ 0.5) and prograde velocities. These stars, which make up ∼ 11 per cent of the total sample, have metallicities at least as low as [Fe/H] = –4.3 and are best interpreted as revealing the existence of a very low-metallicity tail to the Galaxy’s metal-weak thick disc population (e.g. Chiba & Beers 2000). On the other hand, the low-e retrograde stars that have |$|Z_{\rm max}| \le 3\, \rm kpc$| (∼ 4 per cent of the sample) are most likely an accreted population. We also find a population (∼ 6 per cent of the sample) of low |Z_max| stars that have high eccentricity orbits (median ≈ 0.88) with small pericentres and which are split equally between prograde and retrograde motion. It seems likely that many of these stars might be associated with the Gaia Sausage accretion event (Myeong et al. 2019; Yuan et al. 2020; Koppelman et al. 2019). With the possible exception of the low-e, low |Z_max| prograde stars that may have a somewhat lower mean [α/Fe] abundance ratio, none of the four subpopulations with low |Z_max| are distinguished, as regards [α/Fe] or [Fe/H], from the full set of stars for which high-dispersion-based analyses are available.

Finally, we find that a small fraction of our sample (21 stars, ∼ 4.4 per cent) are likely to be escaping from the Galaxy, that is, are on orbits that are not bound. The [Fe/H] and [α/Fe] distributions of these stars are not distinguished from those for the full sample; for example, their metallicities are spread from –4 (or less) to –2 in [Fe/H]. Our preferred interpretation for these stars is that they have acquired sufficient energy to escape from the Galaxy via the gravitational interaction that occurs when infalling dwarf galaxies are tidally disrupted by the MW.

SUPPORTING INFORMATION

suppl_data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

We thank the referee for their comments on the original version of the manuscript. This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research innovation programme (grant agreement ERC-StG 2016, no. 716082 ‘GALFOR’, PI: Milone, http://progetti.dfa.unipd.it/GALFOR). ADM is supported by an Australian Research Council Future Fellowship (FT160100206). AFM acknowledges support by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie (grant agreement no. 797100). ARC is supported in part by the Australian Research Council through a Discovery Early Career Researcher Award (DE190100656). AMA acknowledges support from the Swedish Research Council (VR 2016-03765), and the project grant ‘The New Milky Way’ (KAW 2013.0052) from the Knut and Alice Wallenberg Foundation. Parts of this research were supported by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. APM acknowledges support from MIUR through the FARE project R164RM93XW SEMPLICE (PI: Milone) and the PRIN program 2017Z2HSMF (PI: Bedin). K.L. acknowledges funds from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 852977). AF acknowledges support from NSF grant AST-1255160.

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, http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in Gaia Multilateral Agreement.

DATA AVAILABILITY

The data underlying this article are available in the article and in its online Supporting Information.

Footnotes

REFERENCES

APPENDIX A: CLUSTERING ANALYSIS

As an alternative approach, we analysed our sample of stars exploiting the scikit-learn Spectral clustering algorithm (Pedregosa et al. 2011). This analysis has the advantage of being almost completely independent from our choices (i.e. limiting |Z_max|, prograde/retrogade, etc.), while on the other hand has the limitations of being a ‘blind’ analysis. Specifically, in order for the result to be trustworthy, we should have knowledge of how the selection biases [e.g. bright stars in the solar neighbourhood, uncrowded stellar environments and the other effects discussed in Da Costa et al. (2019)] propagate into the cluster choices. Unfortunately, it is not possible to quantify these biases.

None the less, we believe that it is worth exploring this fully independent classification approach. The clustering has been performed in the 4D space with the three actions and the eccentricity (⁠|$J_{\rm R},\, J_{\rm \phi },\, J_{\rm Z},\, e$|⁠) with the following input parameters: affinity=nearest_neighbors and assign_labels=discretize.

Fig. A1 shows the results of the clustering analysis. In the top three panels, we show the same top three panels as for Fig. 10 but now for the eight identified clusters. We note that group number 1 [red markers with coordinates (0, −1) in the action map] is the group of unbound stars, and therefore is not investigated further in the lower panels.¹⁸ Panels (a1)–(e2) in Fig. A1 and panels (f1)–(j2) in Fig. A2 then show the kernel density distributions of different orbital parameters for each subgroup, with groups 2–4 in the left-hand panels and groups 5–8 in the right-hand panels. A comparison of Figs A1 and A2 with Fig. 10 suggests the following.

• Group number 1 (G1, red markers, 28 stars) is composed of stars with energies consistent or greater than 0, visible in the top right panel of Fig. A1. These stars have apparent apogalacticon distances larger than the MW virial radius, that is, 250 kpc. They are characterized by large values of the radial action, J_R, which translates in J_ϕ/J_tot ∼ 0 and (J_Z − J_R)/J_tot ∼ −1. Overall, we find that 28 stars are grouped in G1. All of these stars are in the subsample discussed in Section 5.3, although two stars in the Section 5.3 subsample, namely SMSS J095211.09–185713.7 and SMSS J002148.06–471132.1, are not classified as G1 stars, despite having apparent D_apo larger than 250 kpc. These two stars have negative (bound) energies and are classified by the clustering algorithm in group number 7 (pink points), which is composed of stars on loosely bound orbits. There is therefore essentially no discrepancy between the subsample discussed in Section 5.3 and the high-energy G1 stars identified by the clustering algorithm. We do not show the kernel density distributions of these stars in the subsequent panels for scaling reasons.

• Group number 2 (G2, blue markers, 93 stars) are a combination of prograde and retrograde halo stars. We note that they span quite a wide energy range, but we find difficult to draw further conclusions. Presumably this group is made up of a mixture of in situ and accreted halo stars.

• Group number 3 (G3, green markers, 54 stars) is a mixture of prograde and retrograde stars on loosely bound orbits, that venture far from the Galactic plane.

• Group number 4 (G4, purple markers, 61 stars) partially overlaps with low-e stars in all three of the top panels of Fig. 10. Furthermore, their eccentricity distribution peaks at e ∼ 0.4–0.6, while most of them remain roughly confined within |$5\, \rm kpc$| of the Galactic plane and 20 kpc from the GC. A possible interpretation would be to consider these stars as thick-disc stars. This hypothesis is also supported by the Toomre diagram in the top centre panel of of Fig. A1, where purple stars occupy a locus typical of thick-disc stars. By the comparison with Fig. 10 we find a partial match of this group with the low-e prograde population (orange points), identified as candidate very metal-weak thick disc stars.

• Group number 5 (G5, yellow markers, 53 stars) partially shares the location of the low-e and retrograde population (red dots in Fig. 10) as well as partially overlapping with the locus of the Sequoia remnants identified in Myeong et al. (2019). The eccentricity distribution peaks at about e ∼ 0.6, and they are confined to the inner halo (⁠|$D_{\rm apo}\le 10\, \rm kpc$|⁠). Panel (b1) of Fig. A1 shows that |$\sim 80{{\ \rm per\ cent}}$| of these stars are confined within |$5\rm \, kpc$| from the Galactic plane. Most of them do not orbit further than 10 kpc from the GC. Comparing then the top right panels of Fig. A1 with Koppelman et al. (2019, bottom right panel of their fig. 2), we see that G5 stars have a higher energy (in absolute values) than Gaia Sequoia stars, while their energy and their angular momentum suggest a possible association with the Thamnos 1/2 groups (Koppelman et al. 2019).

• Group number  6 (G6, brown markers, 37 stars) is composed of stars with very retrograde and mildly eccentric orbits that venture far from the Galactic plane and from the GC, with Z_max and D_apo peaking at ∼15–20 kpc. These stars have energies that range from ∼−1 to |$\sim -1.7\,\,[10^5\, \rm kpc\cdot km^2\, s^{-2}]$|⁠. This group is consistent with the identification of the Gaia Sequoia remnants in Koppelman et al. (2019).

• Group number 7 (G7, pink markers, 67 stars) partially overlaps with G4 both in the action map and in the Toomre diagram, while it is well defined in the energy versus J_ϕ-plane. Panel (a2) of Fig. A1 shows that the eccentrity distribution of these stars is double peaked, with the first peak at e ∼ 0.3 and a second one at e ∼ 0.6. Given the distribution of D_apo and Z_max it would seem that these stars are a mixture of halo and thick-disc stars, with lower binding energies than their G4 counterparts. As for G4, we note that there is a clear overlap between this group and the low-e prograde stars shown in Fig. 10.

• Group number 8 (G8, grey markers, 82 stars) share roughly the same location of the high-e both retrograde and prograde populations (azure and navy dots in Fig. 10) in all top three panels of Fig. A1. Furthermore, the distribution of their orbital parameters nearly overlaps with those of the combined high-e prograde and retrograde populations. Their energy and angular momentum agrees with the Gaia Sausage definition in Yuan et al. (2020). We find particularly interesting the pericentre/apocentre distributions, whose analysis suggest that most of these stars move back and forth from the GC to the Galactic outskirts, always remaining within few |$\rm kpc$| from the Galactic plane (∼60 per cent these stars are indeed confined within |$5\, \rm kpc$| from the plane). Given the observed orbital properties, and in particular the perigalacticon distances as low as ∼1 kpc, we speculate that such a remnant can be formed via a ‘head-on’ accretion event, as in the Sausage progenitor (Myeong et al. 2019).

It is clear that the two different analyses (discussed here and in Section 5) reach qualitatively the same conclusions. First, we find solid evidence for the existence of a very metal-weak component in the Galactic thick disc. Further, from the analysis of the orbital actions and by means of the action map, we have identified possible members of the Gaia Sequoia and Gaia Sausage accretion events. The analysis also suggests that the low |Z_max|, high-e population that is composed of stars with both prograde and retrograde orbits, may also be associated with the Gaia Sausage event. Both analyses also identify a consistent set of candidates that are likely not bound to the Galaxy.

APPENDIX B: EXAMPLE ORBITS

In the following, we show four typical orbits of Halo, very metal-weak thick disc, Sequoia and Sausage stars within our sample. Each orbit is colour coded according to the integration time, while the white dot represents the current position in the Galactocentric Cartesian reference frame. As discussed in Section 3.2, each orbit has been integrated in a McMillan2017 potential (McMillan 2017) backward and forward in time for 2 |$\, \rm Gyr$|⁠. The actions have been computed with the Stäckel fudge method implemented in galpy.