The Metallicity Gradient and Complex Formation History of the Outermost Halo of the Milky Way

Metal-poor stars have been widely recognized as "fossils" of the earliest generation of stars in the universe and as important tracers of the assembly history of our Galaxy. In particular, the most metal-poor stars are likely to be among the most ancient, possibly even true second-generation stars (e.g., Beers & Christlieb 2005; Frebel & Norris 2015; Hansen et al. 2016; Hartwig et al. 2018). Discerning the origins of different metal-poor populations in the Milky Way is crucial for understanding when, where, and how different Galactic components formed. Consequently, our understanding of the assembly history of the Milky Way has evolved significantly over the past half century, in large part through studies of metal-poor stars.

1.1. Metallicity Gradient

1.2. Galactic Dual Halo Formation History

2.1. Non-local "In Situ" Samples

Our in situ samples consist of two sets of SDSS (R ∼ 2000) giants compiled by Chen et al. (2014) and Janesh et al. (2016). The Chen et al. (2014, hereafter C14) sample comprises 15,723 red giant branch (RGB) stars from SDSS DR9 (Ahn et al. 2012), compiled to study the thick disk, inner-halo, and outer-halo of the Galaxy. The Janesh et al. (2016, hereafter J16) sample is made up of 6036 K giants from SDSS DR9, selected to study substructure in the stellar halo. Both samples cover distance ranges in excess of 100 kpc, making them good candidate data sets for studying the outermost halo. The majority of stars in C14 and J16 do not have reliable Gaia parallaxes (with <20% uncertainty) available. We use the original distances derived in C14, and calculate photometric distances for J16 following the method described in Beers et al. (2000). We also retain the original SDSS radial velocities (mean uncertainty ∼2 km s⁻¹) and proper motions for both samples. Spatial distributions for these samples are included for reference in Figure 1. The C14 sample is highly concentrated within ∼20 kpc of the Galactic center while the J16 sample is distributed more uniformly. Though the majority of both samples lies within r ∼ 50 kpc, each spans a range of >100 kpc.

Zoom In Zoom Out Reset image size

2.2. Local "Ex situ" Samples

The primary ex situ sample used for our analyses is compiled from SDSS DR15 (Aguado et al. 2019). The initial query to the SDSS catalog server³ resulted in 357,816 stars with signal-to-noise ratios (S/N) > 10 in the effective temperature range 4500 K < T_eff < 7000 K, where the SEGUE Stellar Parameter Pipeline (SSPP; Lee et al. 2008) is most reliable. Duplicate stars were removed by choosing the measurement with the highest S/N. Stars with spectra taken on plug-plates which were part of SEGUE cluster- or structure-targeting programs were removed prior to analyses.⁴

Proper motions and radial velocities were taken from Gaia DR2 (Gaia Collaboration et al. 2018) where available. When unavailable, we used the original kinematic parameters from the SDSS archive. The resulting sample contains mainly (>99%) Gaia proper motions and SDSS radial velocities. We adopt distances from the Bailer-Jones treatment of the Gaia parallaxes, and restrict the sample to stars with <20% distance uncertainty (Bailer-Jones et al. 2018). Stars with uncertainties on their radial velocities exceeding 20 km s⁻¹ were removed from the sample. The resulting sample has a mean radial-velocity uncertainty of ∼1.5 km s⁻¹. We limit our sample to a local volume within 4 kpc of the Sun, leaving us with 118,037 stars.

Our complementary ex situ sample was compiled by Huang et al. (2019) from DR1 of the SkyMapper Southern Survey (SMSS; Wolf et al. 2018). The authors provide metallicity estimates for 972,994 RGB stars and compile kinematic parameters where available. Proper motions and (Bailer-Jones) distances from Gaia DR2 are available for the majority (∼70%) of the sample. Of the 972,994 stars in this sample, 423,995 have available radial-velocity estimates, compiled by Huang et al. (2019) from a variety of catalogs, primarily Gaia DR2 and the Galactic Archaeology with HERMES survey (GALAH; R ∼ 28,000; Buder et al. 2018). The resulting sample has a mean radial-velocity uncertainty of ∼1 km s⁻¹. After limiting the sample to a local volume, we are left with a sample of 395,144 stars for which viable kinematics can be obtained.

A comparison of the g-band magnitude distributions for all samples used in our analyses is shown for reference in Figure 2. We divide our samples into two plots based on survey source, as the SDSS g-band and the SMSS g-band differ from each other⁵ and cover a different magnitude range.

Zoom In Zoom Out Reset image size

Kinematic parameters, such as U, V, and W, are derived using the galpy Python package (Bovy 2015). In this work we adopt R_⊙ = 8 kpc as the Sun's distance from the Galactic center projected onto the Galactic plane. We use v_LSR = 220 km s⁻¹ for the velocity of the local standard of rest (LSR; Kerr & Lynden-Bell 1986) and (U, V, W)_⊙ = (−9, 12, 7) km s⁻¹ (Mihalas & Binney 1981) for the motion of the Sun with respect to the LSR.

Orbital parameters such as energy, maximum height from the Galactic plane (Z_max), and r_apo are derived with a version of the potential code used in Chiba & Beers (2000), adapted to run in Python. This code adopts the analytic Stäckel potential developed by Sommer-Larsen & Zhen (1990), consisting of a flattened, oblate disk and a nearly spherical massive halo. We note that we use the Stäckel potential and not galpy's Milky Way-like potential, MWPotential2014, following the recent comparison of Kim et al. (2019) between the Stäckel and galpy potential. Kim et al. (2019) suggest that MWPotential2014 may not be the ideal choice for studies focusing on the outer-halo because many highly energetic, predominantly retrograde (outer-halo) stars are found to be unbound when using the shallower galpy potential.

We estimate uncertainties on orbital values produced by the Stäckel potential via a Markov Chain Monte Carlo sampling method. Assuming that the uncertainties on the input parameters are normally distributed about their observed values, we generate 1000 randomly sampled orbits per star, and adopt the standard deviation of the resulting distribution as our uncertainty. Since our in situ non-local samples do not utilize high-precision Gaia data, they have larger uncertainties on average than our ex situ samples. We choose not to make any cuts on uncertainty for C14 and J16 because such cuts would exclude an overly large number of stars from our analyses. Our two ex situ samples utilize Gaia data, so we are able to trim high-uncertainty stars from our data to minimize potentially spurious features.

The manner in which we prune high-uncertainty stars from our samples depends on the subsequent analyses we intend to perform on them. In Section 4.1.2 we bin stars over r_apo in steps of size 5 kpc each, so we choose to restrict our sample to stars with uncertainty on r_apo < ±10 kpc. This allows us to create a low-uncertainty sample without losing too many distant halo stars due to an overly strict cutoff. To minimize contamination from disk-system stars in our local samples, we exclude stars with r_apo < 10 kpc and Z_max < 3 kpc, leaving us with 10,078 SDSS halo stars and 6576 SMSS halo stars.

Select values for our compiled ex situ samples are available in Tables 1 and 2. The final ex situ samples' derived orbits span distances up to Z_max ∼ 60 kpc and r_apo ∼ 80 kpc.

4.1. Global Metallicity Gradient

4.1.1. In situ Results

In Figure 3 we construct MDFs for C14 (left panels) and J16 (right panels) over increasing slices of r. The MDF of the C14 sample shifts toward the metal-poor regime and its metal-poor tail noticeably increases in relative proportion as we move farther from the Galactic center, particularly beyond 45 kpc. Though this effect is not as noticeable in the J16 sample as the C14 sample, likely due to the J16 selection function,⁶ both clearly show that the dominant stellar component changes from the metal-richer populations of the metal-weak thick disk (MWTD; [Fe/H] ∼ −0.8 to −1.8) and inner-halo ([Fe/H] ∼ −1.6) to the metal-poor outer-halo ([Fe/H] ∼ −2.2) over increasing distance. We note that there exists, interestingly, a relatively strong inner-halo like population even beyond 30 kpc in both samples. This may be associated with a major merger event, and is discussed in more detail in Section 4.1.2.

Zoom In Zoom Out Reset image size

Although our in situ samples are likely to be affected by metallicity selection biases, here we are interested in the nature of their lowest-metallicity tails, with [Fe/H] < −2.0. As clearly shown by the empirical comparison of giant-branch luminosity with metallicity for Galactic globular clusters in Figure 5 of Huang et al. (2019), this dependency is minimal at the lowest abundances. Additionally, the increasingly apparent bimodality at larger distances in Figure 3 (also seen in Carollo et al. 2007, 2010) cannot be explained by this bias alone.

4.1.2. Ex situ Results

We conduct a similar analysis using the local ex situ SDSS and SMSS samples to mitigate the metallicity-distance bias, based on the suggestive evidence of a possible metallicity gradient in the outermost halo from the non-local in situ samples analyzed above.

We examine how the frequency of VMP stars and the average metallicity of our samples vary as a function of r_apo, as shown in the upper and lower panels of Figure 4, respectively. Stars are binned in steps of size 5 kpc until fewer than 10 stars are available per bin, after which all subsequent bins are combined. We note that, because of this binning choice, the final bin in each panel may be influenced by low-number statistics and suffer from high uncertainty.

Zoom In Zoom Out Reset image size

The VMP frequency of the SDSS sample shown in the top-left panel of Figure 4 climbs very slowly in the range r_apo = 10–40 kpc, after which it experiences a sharp increase, rising from ∼20% to ∼60% over the next 10 kpc. The average metallicity slowly decreases from [Fe/H] ∼ −1.4 at r_apo = 10 kpc, plateauing around [Fe/H] ∼ −1.6 at r_apo = 25–40 kpc, then dropping rapidly to [Fe/H] ∼ −2.0. We note that the statistics in the largest distance bins of the left panels appear contrary to the overall trends, but have high uncertainty compared to the majority of the preceding bins.

The VMP frequency for the SMSS photometric sample experiences a steady climb over r_apo, maxing out at ∼30%, as seen in the top-right panel of Figure 4. The average metallicity decreases from [Fe/H] ∼ −0.8 at 10 kpc to [Fe/H] ∼ −1.6 at 40 kpc. Although the changes in VMP frequency and mean metallicity are not as dramatic for the SMSS sample as for the SDSS sample, both samples display a clear metallicity gradient over r_apo. The differences in the samples are likely due to the fact that the SMSS sample does not reach as far into the halo as the SDSS sample.

The steep change in mean metallicity from [Fe/H] ∼  −1.6 to approximately −2.0 at r_apo  ≳ 40 kpc in the SDSS panels in Figure 4 could indicate a shift between the relative dominance of the inner-halo ([Fe/H] ∼  −1.6) and outer-halo ([Fe/H] ∼  −2.2), as discussed in Section 1.2. Another possibility is that this region of the metallicity distribution represents the shift between stars donated to the outer-halo by the Sequoia merger event, estimated at [Fe/H] ∼ −1.6 (Myeong et al. 2019), and stars donated by smaller mergers of more metal-poor satellites. Accordingly, we investigate the detailed Galactic halo assembly history using these local samples in the next subsection.

4.2. Detailed Accretion History

4.2.1. Metallicity Distribution

Since the behavior of the SDSS local halo stars changes at approximately 35–40 kpc (e.g., Figure 4, left panels), we further investigate the reason for this behavior by constructing MDFs for "near" (10 ≤ r_apo < 35 kpc) and "far" (r_apo ≥ 35 kpc) halo samples, shown in the left panels of Figure 5. A similar diagram is shown for the SMSS halo stars in the right panels of Figure 5. The SMSS halo stars do not present the same sharp changes at 35–40 kpc (see Figure 4, right panels), in part because they do not probe the same distance range as the SDSS sample.

Zoom In Zoom Out Reset image size

However, since the SMSS sample displayed similar general characteristics to the SDSS sample in our analyses of Figure 4, we chose to create additional near/far MDFs from the combination of both the SDSS and SMSS samples, in order to bolster the number of stars available in the r_apo ≥ 35 kpc range (see Figure 6). Gaussian distributions are fit to the combined data using the scikit-learn Gaussian Mixture Model package in Python to identify components with potentially distinct origins. The near-halo combined MDF primarily consists of a distinctive component at [Fe/H] ∼ −1.4, with a smaller, more metal-rich peak around −0.6, possibly belonging to a portion of the MWTD that was not completely removed by the cuts made in Section 3. The far-halo combined MDF also has a dominant [Fe/H] ∼ −1.4 peak, as well as a more metal-poor component at [Fe/H] ∼  −2.3. The near halo also possesses a small VMP population, fitted with a peak at [Fe/H] ∼  −2.0, but this population rises in relative significance in the far halo.

Zoom In Zoom Out Reset image size

Figure 6 shows that there may be at least two separate populations at r_apo ≥ 35 kpc. The more metal-rich peak ([Fe/H] ∼  −1.4) seen in the far-halo (lower) panel of Figure 6 could be a selection of inner-halo population stars still present at r_apo ≥ 35 kpc. The exact location of the transition zone between the inner- and outer-halo regions is uncertain. Carollo et al. (2007) place it at ∼15–20 kpc while Kim et al. (2019) give an estimate of ∼30 kpc (both use the same Stäckel potential adopted in this work), so it is possible that we could still see evidence of the inner-halo population at r_apo ≥ 35 kpc. Another possible interpretation is that this [Fe/H] ∼ −1.4 peak comprises stars accreted from the Sequoia and Gaia-Sausage mergers—this could explain the hint of bimodality seen in this component (two sub-peaks at [Fe/H] ∼ −1.3 and [Fe/H] ∼ −1.5)—while the more metal-poor peak represents stars accreted from a series of more minor mergers.

4.2.2. Prograde versus Retrograde

We can examine these trends in more detail by dividing our local ex situ samples on rotational velocity (v_ϕ) into prograde (v_ϕ > 0 km s⁻¹) and retrograde (v_ϕ < 0 km s⁻¹) components, as shown Figures 7(a) and (b), respectively. There are 6941 prograde stars and 3137 retrograde stars in the SDSS sample, and 5267 prograde stars and 1307 retrograde stars in the SMSS sample.

Zoom In Zoom Out Reset image size

The prograde SDSS sub-sample (Figure 7(a), left panels) experiences a slow climb in VMP frequency up until r_apo ∼ 35 kpc, after which the fraction of VMP stars climbs to ∼50%. Its average metallicity exhibits a steady decrease from 10 to 45 kpc. The retrograde sub-sample (Figure 7(b), left panels) exhibits no discernible relationship between VMP frequency and r_apo in the 10–45 kpc range, and its average metallicity hovers around −1.6 in the same range. Once more, the contrary behavior of the final bins of these sub-samples may be due to low-number statistics (high uncertainty).

The prograde and retrograde SMSS sub-samples (Figures 7(a) and (b), right panels) exhibit largely the same behavior as the SDSS sub-samples. The prograde stars show a strong dependence on r_apo for both VMP frequency and mean metallicity, while the same quantities for retrograde stars show no noticeable dependence on distance. The Gaia Sausage ([Fe/H] ∼  −1.3, L_z ∼ 0 kpc km s⁻¹; Belokurov et al. 2018; Myeong et al. 2018) dominates the retrograde (and prograde) signal in the near halo and the Sequoia progenitor galaxy ([Fe/H] ∼  −1.6, L_z ∼ −2000 to −3000 kpc km s⁻¹; Matsuno et al. 2019; Myeong et al. 2019) likely donated the majority of the retrograde stars in the far halo. The similar peak metallicities and overlapping metallicity ranges (see Figure 2 in Matsuno et al. 2019 and Figure 9 in Myeong et al. 2019) of these imported populations could result in this overall metallicity plateau.

While the Sequoia (M_* ∼ 10⁷ M_⊙) event may have imprinted a bulk retrograde signal onto the outer-halo, stars donated by numerous small accretion events likely contributed both prograde and retrograde stars to the outer-halo. However, it is likely more difficult to detect (low-metallicity) stars from small mergers in the retrograde outer-halo due to the overwhelming presence of Sequoia stars. In contrast, numerous minor accretions could be the predominant contributors to the prograde outer-halo, based on the metallicity gradient and strongly increasing VMP fraction seen in Figure 7. If this behavior accurately reflects the assembly history of the outermost halo, observational efforts to compile catalogs of the most ancient, most metal-poor stars may have more success targeting candidates in the prograde rather than the retrograde outermost halo.

We compiled a set of in situ (∼21,700 stars in total) and ex situ (∼16,500 stars in total) samples to confirm the presence of a metallicity gradient in the outermost halo of our Galaxy and explore the complex assembly history of the Galactic halo. The results of the in situ analyses are suggestive regardless of the metallicity-distance selection bias. In our ex situ analyses, we find clear evidence of both a negative metallicity gradient over r_apo and an increasing relative fraction of VMP stars with distance. In particular, the local SDSS sample exhibits a VMP frequency that reaches ∼60% at 50 kpc, commensurate with theoretical studies (i.e., Tissera et al. 2014).

When splitting our samples into prograde and retrograde components, we find that the retrograde appears to exhibit no metallicity-distance correlation while the prograde experiences a steady decline in [Fe/H] and a strong increase in VMP frequency with distance. This may be due to the influence of a more-massive merger (metal-richer satellite) in the retrograde direction, versus numerous minor accretions (metal-poorer counterparts) in the prograde direction. As a result, the prograde outermost halo may be the best place to search for the most metal-poor stars.

Since we have placed tight constraints on uncertainty, we may have excluded some stars in the outer-halo that could have given our analysis more significance, but we are also not as likely to detect completely spurious features even at large distances. In addition, our local, ex situ samples are not susceptible to the metallicity-distance bias that may affect our non-local, in situ results.

We note here that, recently, Conroy et al. (2019b) published an exploration of the Galactic halo using a sample of some ∼4200 giants from the H3 Spectroscopic Survey (Conroy et al. 2019a), and claimed that no metallicity gradient is detectable in their sample. However, the results of their analyses are not dissimilar from our own findings. Although they find a flat metallicity relation across the majority of the halo, they admit possible evidence for a decreasing mean metallicity beyond r ∼ 50 kpc, which may coincide with the behavior in the r_apo ≥ 35 kpc region examined in this work. The VMP component they identify as potentially originating from multiple distinct populations parallels our own hypothesis of a halo component at [Fe/H] < −2.0 consisting of numerous accretions of small mini-halos.

Have we fully explored a volume that could qualify as a comprehensive "outermost" halo, up to the outskirts of the Galaxy? The exact bounds of the outer-halo population are not yet known, and though this work shows the potential for a signature that may extend beyond the volume explored here, further efforts are required to quantify the behavior of the outer-halo MDF beyond 50 kpc. For example, improved kinematics from Gaia DR3 will allow us to expand the narrow magnitude ranges probed by our ex situ samples (see Figure 2). Near-future observations with the Large Synoptic Survey Telescope (LSST; Ivezić et al. 2019), combined with spectroscopic follow-up, as well as large spectroscopic surveys undertaken with the Dark Energy Spectroscopic Instrument (DESI Collaboration et al. 2016), the WHT Enhanced Area Velocity Explorer (Dalton et al. 2014), and the 4 m Multi-Object Spectroscopic Telescope (Helmi et al. 2019), will enable a more thorough exploration of the distant halo and the significance of its metallicity gradient. Further analyses of the VMP catalogs from the LAMOST (Cui et al. 2012) survey (e.g., the DR3 VMP catalog of Li et al. 2018 ⁷ ), should prove illuminating as well.

It would be of particular interest to examine the prograde/retrograde metallicity distributions with a sample of distant CEMP (Beers & Christlieb 2005) stars. Similarly to metal-poor stars, CEMP stars have been recognized as useful tracers of Galactic formation history (e.g., Beers & Christlieb 2005; Carollo et al. 2010, 2014; Lee et al. 2017, 2019; Yoon et al. 2018, 2019), and their individual nucleosynthetic sub-classes can provide even more information about the origins of various stellar populations. CEMP stars that display over-abundances of elements associated with the slow neutron-capture process (CEMP-s) are thought to originate in more-massive galaxies, while CEMP stars that exhibit no over-abundances of neutron-capture products (CEMP-no) are thought to originate in less-massive galaxies (e.g., Lee et al. 2017; Yoon et al. 2018, 2019). A comparison of the CEMP-s to CEMP-no ratios (which can be identified using absolute carbon abundance criteria, readily measured with medium-resolution spectroscopic data, see Yoon et al. 2016) in the prograde and retrograde outermost halo could help clarify the origins of the populations present there.

We thank Xiang Xiang Xue for access to the J16 data set, Young Sun Lee for his assistance and discussions regarding SDSS data, and the referee for their careful consideration of this work. The authors acknowledge partial support from grant PHY 14-30152, Physics Frontier Center/JINA Center for the Evolution of the Elements (JINA-CEE), awarded by the US National Science Foundation.

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

Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS website is http://www.sdss.org/.

The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington.

Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III website is http://www.sdss3.org/.

SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.

Software: astropy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018), galpy (Bovy 2015), numpy (van der Walt et al. 2011), matplotlib (Hunter 2007), scipy (Virtanen et al. 2020), scikit-learn (Pedregosa et al. 2011), STILTS (Taylor 2006).