The Origin of the Milky Way's Halo Age Distribution

The stellar halos of galaxies play a particularly crucial role in understanding the early formation of galaxies due to the dynamical and chemical fingerprints that they carry; fingerprints that shed light on a given galaxy's assembly history and the chemical history of its accreted satellite galaxies. In the Milky Way (MW), the advent of massive photometric and spectroscopic surveys has revolutionized the understanding of our own stellar halo. With respect to the diffuse halo, Carollo et al. (2007, 2010), and Beers et al. (2012) demonstrated that it comprises at least two stellar components with different kinematics, dynamics, and chemical composition (the inner and outer halos). These findings have now been supported by a large number of authors (de Jong et al. 2010; Nissen & Schuster 2010; Deason et al. 2011; An et al. 2013, 2015; Hattori et al. 2013; Kafle et al. 2013; Das et al. 2016), and most recently by Helmi et al. (2016) using Gaia data. The MW's stellar halo is also populated by numerous stellar streams and overdensities, the product of past mergers of satellite galaxies (Bell et al. 2008; Helmi 2008), suggesting a complex superposition of stellar populations.

Another crucial parameter that can provide a comprehensive picture of the MW assembly process is the age distribution of the halo system. Recently, Santucci et al. (2015) determined the age of the underlying stellar populations in the MW's halo out to ∼25 kpc from the Galactic center by employing a small sample (∼4700) of spectroscopically confirmed blue horizontal-branch (BHB) stars selected from the Sloan Digital Sky Survey (SDSS-DR8; Aihara et al. 2011). The resulting age map showed that in the central region of the MW, there exists a region (out to 10–15 kpc) composed of a significant concentration of very old stars (∼11.5–12.5 Gyr), with younger structures or overdensities extending out to greater distances (likely associated with the Sagittarius stream and Virgo overdensity). In a later analysis, Carollo et al. (2016) produced a high-resolution age map extending out to ∼60 kpc, using a much larger sample (∼130,000) of photometrically selected BHB stars from SDSS-DR8. In addition to the aforementioned central region dominated by old stars, an inferred age gradient of −25 ± 1 Myr kpc⁻¹ was derived, and numerous known and unknown structures and overdensities were identified. The existence and claimed value of a negative age gradient was confirmed by Das et al. (2016) in an investigation based on a sample of spectroscopically identified BHB stars from SDSS-DR8.

In cosmological simulations of MW-mass galaxies, it is possible to follow the evolution of the satellite galaxies that merge as part of a system's assembly (Brook et al. 2004a; Bullock & Johnston 2005; Cooper et al. 2010, 2015; Few et al. 2012, 2014; Gómez et al. 2012; Tissera et al. 2014). In this scenario, stellar halos are predicted to form primarily through the accretion of satellite galaxies, each with different stellar masses, gas fractions, and stellar population distributions (e.g., Tissera et al. 2014; Cooper et al. 2015). An important contribution from in situ stars within the inner ∼15 kpc has been also identified in hydrodynamical simulations, which have a range of possible origins (e.g., Brook et al. 2004b; Zolotov et al. 2009; Font et al. 2011; House et al. 2011; Tissera et al. 2013; Cooper et al. 2015; Pillepich et al. 2015; Monachesi et al. 2016b).

Recent simulations of MW-mass galaxies provide information on the age distribution of stars formed both in situ and in accreted satellites, however, most such studies do not include a discussion of the overall age gradients (e.g., Brook et al. 2004b; Tissera et al. 2013; Pillepich et al. 2015). This was in part driven by the lack of any strong empirical constraints with which to compare. One early exception is the semi-cosmological, sticky-particle work of Bekki & Chiba (2001), who found an age gradient of −30 Myr kpc⁻¹ over the galactocentric radial range of 20–50 kpc in one low-resolution simulation (see also Samland 2004).

In this Letter we focus on the age structure of the stellar halo system in a set of five simulated MW-mass-sized halos from the Aquarius Project. Understanding the age structure of the stellar halos can contribute to reveal their origin and evolution. These stellar halos have been considered extensively in previous papers, including their chemical abundance patterns, density and metallicity profiles, and assembly histories (Tissera et al. 2012, 2013, 2014, 2016). Here, the age gradients and the relative contributions of the in situ and accreted components, as a function of the galactocentric radius, are explored. The aim of this analysis is to examine to what extent the simulated Aquarius halos are able to reproduce the age trends found in the MW's halo, and to investigate the connection between the age profiles and halo assembly. For these purposes, we select a subset of the level-5 runs from Scannapieco et al. (2009), specifically Aq-A, Aq-B, Aq-C, Aq-D, and Aq-G. Although none of the five realizations underwent a recent major merger, this simulated set of galaxies produce an excess of stars, and consequently, more massive stellar halos than the MW. Such issues affect other simulations as well (e.g., Font et al. 2011; Pillepich et al. 2015), and should be addressed in future work, given that the MW is an ${{\rm{L}}}^{* }$ spiral galaxy with perhaps the least massive stellar halo for its given total mass (Harmsen et al. 2017). A more efficient supernova feedback can contribute to decrease the stellar mass fraction, and form more extended disk structures (e.g., Aumer & White 2013; Pedrosa & Tissera 2015). We note, however, that the measured mass of the MW stellar halo may be underestimated by current observational methods (Sanderson et al. 2017). The analysis presented in this Letter apply a new observational constraint: this is the MW stellar age distribution, a critical piece of information for the confrontation of models and observations. It also contributes to the understanding of the origin of the MW's stellar halo. This Letter is organized as follows: Section 2 provides a short description of the simulations, while Section 3 describes the analysis of the age distributions and associated age maps. The summary and conclusions are reported in Section 4.

We analyzed a subset of five MW-mass galaxies from the Aquarius Project run with a version of GADGET-3, an optimized version of GADGET-2 (Springel 2005) that was modified to include supernova feedback and chemical evolution by Scannapieco et al. (2005, 2006). This code allows the description of a multiphase medium and the triggering of mass-loaded galactic outflows without introducing mass-dependent scaling parameters. The chemical model describes the enrichment by Type II (SNII) and Type Ia (SNIa) supernovae according to the nucleosynthesis yields of Woosley & Weaver (1995) and Thielemann et al. (1993), respectively. These simulations are explained in detail by Scannapieco et al. (2009), and have been used throughout our ongoing series of papers (Scannapieco et al. 2009, 2010; Tissera et al. 2012, 2013, 2014, 2018); here, we only provide their main characteristics. The initial conditions are consistent with a ΛCDM cosmology having the following parameters: ${{\rm{\Omega }}}_{{\rm{m}}}=0.25,{{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.75,{{\rm{\Omega }}}_{{\rm{b}}}=0.04,{\sigma }_{8}=0.9,{n}_{s}=1$ and ${H}_{0}=100\ h\ \mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$, with h = 0.73. The dark matter particle mass was 10⁶ M_⊙ h⁻¹ and the initial gas particle mass 2 × 10⁵ M_⊙ h⁻¹. The corresponding gravitational softenings ranged from 0.5 to 1 kpc h⁻¹. The halos selected have not had a major merger since z < 2.

We aligned the disks with the xy plane as described in Tissera et al. (2012) and employed their adopted dynamical decomposition methodology. We measured  = J_z/J_z,max(E) for each star, where J_z is the angular momentum component perpendicular to the disk plane and J_z,max(E) is the maximum J_z over all particles of given total energy, E. A star on a prograde circular orbit in the disk plane has  = 1; stars with  > 0.65 are considered a part of the disk components. Particles that do not satisfy this requirement are taken to be part of the spheroidal components. Motivated by observations of the MW spheroid that exhibit differences in stellar kinematics and chemical abundances as one moves outwards (Carollo et al. 2007, 2010; Zoccali et al. 2008), we separate our spheroids into two stellar populations according to their binding energy. The central spheroid (bulge) is defined by stars more bound than those with the minimum energy (E_cen) at r ≥ 0.5 × r_opt (r_opt is defined as the radius that encloses ∼80% of the stellar mass identified by the SUBFIND algorithm (Springel et al. 2001). Stars more weakly bound than E_cen are taken as part of the stellar halo. In this Letter, and to confront with observations, the stellar halos are not separated into inner and outer components as done in previous works, but rather are taken as an ensemble system. These criteria are chosen so that the definition of the stellar halos adapts to the overall size of each individual galaxy, and is the same definition used in the aforementioned series of Aquarius Project stellar-halo papers. For the analysis that follows, we have removed satellite galaxies identified by the SUBFIND as individual systems, but disrupted satellites and/or any residual stellar streams remain.

Tissera et al. (2013, 2014) analyzed the spatial distributions, chemical abundances, and formation histories of the stellar populations of the Aquarius halos and found evidence that they are mainly built by stars formed in satellite galaxies, and later accreted onto the main halo. However, and as in agreement with previous numerical results, an important contribution of in situ stars are detected in the central regions (Section 1). In addition, Tissera et al. (2014) showed that stellar halos dominated by accreted massive satellite galaxies (M > 10^9.5 M_⊙) exhibit steeper metallicity gradients (see also Cooper et al. 2010). They also reported that low-metallicity stars are mainly contributed by low-mass satellites, and are more frequent in the outskirts of halos. The central regions of such systems (within ∼10 kpc, including the bulge), where a significant contribution of old stars was found, have been analyzed in detail by Tissera et al. (2018). Moreover, one of the simulated halos (Aq-C-5) showed a good match with the spatial, kinematic, and metallicity properties observed in the MW (Zoccali & Valenti 2016). The analysis of the assembly history of the central regions shows that this halo did not accrete satellite galaxies more massive than ∼10¹⁰ M_⊙ during its assembly. This characteristic is also relevant in the analysis of the age gradients in our simulated halos.

For the five simulated halos, we identify those stars that formed in situ or in accreted satellite galaxies by adopting the following criteria: in situ stars were assumed to form within the virial radius, while accreted stars were assumed to form in separate galaxies prior to accretion (i.e., before entering the virial radius of the progenitor galaxy). Tissera et al. (2013) defined three different sub-populations of in situ stars: (1) those born from gas accreted in the first stages of assembly, (2) disk-heated stars formed in the disk structure of the main progenitor galaxy, then heated kinematically, and (3) those formed from gas carried in by gas-rich satellite galaxies (endo-debris). As described in previous works, disk-heated and endo-debris stars exhibit distinct chemical properties that can help to link observations to the galaxy-formation models (Brook et al. 2004b; Tissera et al. 2013). Nevertheless, in this Letter, for the sake of clarity and simplicity, all stars born inside the virial radius are grouped and analyzed as in situ stars, as they dominate the inner region of the stellar halos, while the outer regions are mainly populated by accreted stars.¹⁰

For illustration purposes, Figure 1 shows the smoothed stellar age-map distribution projected onto the (x, z) plane (z is the direction of rotation of the galaxy) for three of the analyzed halos.¹¹ The maps have been built by removing the contribution of the bulge according to the criteria given in Section 2. The upper panels represent the entire population of stars assigned to the halos, while the lower panels show the maps for the accreted stellar population only; both maps extend to a radius of 100 kpc. As can be seen, the halos exhibit clumpy age structure at large distances from the center, reflecting the mixture of stars with different ages. The presence of younger structures and their increasing number with galactocentric distance are globally consistent with the observational results of the MW reported in Carollo et al. (2016). By comparing the age distributions of the entire stellar halo population (upper panels) with those of the accreted stars only (lower panels), different features emerge: when the in situ stars are removed from the age map, an age distribution qualitatively more similar to that observed in the MW's halo is visible, consisting of younger structures at large distances from the galactic center and a concentration of older stars in the central region.

Zoom In Zoom Out Reset image size

Examination of the panels in Figure 1 reveals that the accreted central region of Aq-C and Aq-A is dominated by old stellar populations (dark blue area), in agreement with the Ancient Chronographic Sphere (ACS) observed in the MW by Santucci et al. (2015) and Carollo et al. (2016). We note that the in situ components also contribute with old stars as shown in Tissera et al. (2018). Apart from this common feature, each halo has its own peculiarities. In the case of Aq-C, when the entire stellar-halo population is considered, the central region of the galaxy exhibits a disk-shaped feature, due to the contribution of younger, disk-heated stars. In fact, the mass fraction of halo stars that originated in the disk is ∼40%, consistent with the fractions found by Brook et al. (2004b). Interestingly, Aq-D has a different age distribution, with slightly younger accreted stars, compared to the other two halos.

To quantify the age variations, we estimate the age profiles for each stellar component as a function of galactocentric radius. The age profiles are derived by taking the median values in concentric shells of 2 kpc radial extension. Figure 2 shows the age trends for the in situ (left panel), accreted (middle panel), and combined (right panel) stellar halo populations (all represented by solid lines). The age gradient found in observational studies of the MW is denoted by the green dotted–dashed line. As can be seen, on average the in situ sub-populations possess flat age gradients, and in some cases (Aq-A, Aq-C, Aq-D) a significant contribution of younger stars in the very central regions (10–20 kpc). Such contributions are caused by gas brought in by gas-rich satellite galaxies that enter the virial radius (these stars are classified as the endo-debris population and included in the in situ component) or by disk-heated stars. Each age profile reflects its particular history of assembly.

Zoom In Zoom Out Reset image size

The accreted stellar population exhibits a steeper age gradient in most of the analyzed halos, as can be seen from the linear regressions applied to the age profiles of these components (Figure 2, middle panels). In fact, the inspection of Figure 2 reveals that the accreted components are younger than the in situ stellar populations for Aq-B and Aq-G by 1–3 Gyr, while for Aq-A and Aq-C the accreted component is younger than the in situ one only in the outer halo region (r > 25 kpc) by less than 1 Gyr. The accreted population in Aq-B exhibits a quadratic-shape age gradient, but it does not appear statistically significant, while model Aq-C shows a significant linear age gradient in the range from ∼−11 to ∼−12 Gyr (not labeled in the panel).

In the right column of Figure 2, the total age profiles for the entire stellar halo are shown together with the linear regression fits (solid black and red lines, respectively). Negative age gradients are found for Aq-A, Aq-B, and Aq-C, with values of −8.0 ± 2.0 Myr kpc⁻¹, −12.3 ± 3.5 Myr kpc⁻¹, and −11.8 ±1.0 Myr kpc⁻¹ (not labeled in the panels), respectively, while Aq-D and Aq-G show very weak age gradients. The gradients derived for the entire stellar halos of the Aquarius set are shallower than those determined in observational studies of the MW's stellar halo¹² : −25 Myr kpc⁻¹ (Carollo et al. 2016) and −30 Myr kpc⁻¹ (Das et al. 2016).¹³

It is worth noting that observations of halo stars, both within the MW and in external galaxies, are optimized to reduce disk-component contamination. In particular, observations carried out in external galaxies are generally performed along the major axis of disk rotation (Monachesi et al. 2016a; Harmsen et al. 2017). To better match the observational conditions, we recalculate the age profiles and fractions by excluding all stars within 5 kpc of the mid-plane. By doing this, the influence of stars that might still belong to an extended vertical disk is minimized (Harmsen et al. 2017). The age profiles for these sub-samples (hereafter referred to as the fiducial stellar halos) are represented in Figure 2 with dotted–dashed lines. As can be seen from the figure, most of the discarded stars belong to the in situ stellar populations (larger discrepancies in the in situ age profiles with respect to the original sample). The recalculated age profiles are, in most cases, steeper than those derived by adopting the entire halo sample and are labeled in the panels. The age gradients for the fiducial stellar halo populations (in situ and accreted stars combined) are −6.8 ± 2.8 Myr kpc⁻¹, −19.5 ± 4.0 Myr kpc⁻¹, −19.2 ± 1.3 Myr kpc⁻¹, +0.2 ± 1.0 Myr kpc⁻¹, and −9.6 ±1.3 Myr kpc⁻¹, for Aq-A, Aq-B, Aq-C, Aq-D, and Aq-G, respectively. These slopes are shaped significantly by the particular assembly history of each halo.

The global age trend in a given halo is affected not only by the median age of the in situ and accreted stars, but also by their relative contribution as a function of the radius. Figure 3 shows the stellar mass fraction versus galactocentric radius for the in situ (red), accreted (green), and the total (black) stellar populations in the simulated halos. The stellar mass fraction is calculated within the virial radius for both the entire (solid line) and the fiducial (dotted–dashed line) stellar halos. It is also important to mention that Aquarius stellar halos have a contribution of stars younger than 10 Gyr representing ∼1%–8% of the total stellar halo mass, except for Aq-B (29%), which has experienced a more recent massive accretion. These stars are mostly associated to the in situ component and to the accretion of more massive satellites. Such younger stellar population is not represented in current observations that make use of BHB stars to derive the age structure of the halo system. Nevertheless, in the Aquarius simulations the presence of stars younger than 10 Gyr do not strongly affect the overall trends of the age gradients, except in the central regions.

Zoom In Zoom Out Reset image size

As can be seen in Figure 3, the stellar mass fraction of the in situ population is dominant out to ∼20 kpc for all the halos except Aq-B. Beyond ∼20 kpc, the accreted component dominates over the in situ one at all galactocentric distances, within the virial radius, and for all the simulated halos. This is consistent with the findings of previous works (e.g., Tissera et al. 2014; Cooper et al. 2015). The negative age trends are determined principally by the accreted stars, assembled as the satellite galaxies fall into the potential well of the main galaxies, and then are disrupted. Different mechanisms take place during the assembly of the stellar halos and the mass of the accreted satellite galaxies, as well as their time of accretion, play a major role. In lower-mass satellites, the star formation is truncated earlier due to the gas exhaustion, gas outflows driven by supernovae, tidal stripping, and/or reionization. Such quenching likely occurs before these clumps merged with the main galaxy, thus these satellites possess mainly very old stars. On the contrary, more massive satellites experience a more prolonged star formation activity due to their efficiency in retaining gas in the deeper potential wells. These massive satellites have both young and old populations (Tissera et al. 2014). A combination of accretion time and the mass of the satellites will set the age profile. A negative age profile could arise when low-mass satellites are accreted very early on. These will contribute mostly to the inner regions with their old stars (Tissera et al. 2018). Intermediate- and high-mass satellites accreted later on will have younger stellar populations (because they continued forming stars for longer periods) and their stars will populate both the inner and outer regions. The presence of a larger fraction of the oldest early-on accreted stars in the center of the galaxy (small radii) will set the negative age gradient, but the strength of such profiles depends upon the particular assembly history of each galaxy. The slope is also affected by the generally flat age profile of the in situ stars and their dominance in the inner-halo region. The in situ component comprises a combination of well-mixed stellar populations primarily located in the inner region of the stellar halos. Disk-heated stars populate the inner regions of halos, increasing the fraction of younger stellar populations and contributing to a flattening of the age profile. In some cases, like in Aq-D, a flat age gradient can be generated by the opposite age dependence of the in situ and accreted stellar populations.

The Aquarius stellar halos have greater masses than expected from observations (Harmsen et al. 2017). As the stellar-halo profiles are consistent with an Einasto profile, most of the mass is concentrated in the central regions (Tissera et al. 2014). Part of the mass excess could be due to the inefficient regulation of the star formation activity or an overproduction of disk-heated stars. In the inner halo, the fraction of disk-heated stars can differ from one halo to another, as reported in Table 1 of Tissera et al. (2013; in percentages: 31 (Aq-A), 3 (Aq-B), 24 (Aq-C), 26 (Aq-D), and 35 (Aq-G)). These stars could also originate through a misclassification of thick-disk stars or by the presence of endo-debris stars that contribute significantly (from 20% to 40%), as given in Table 1 of Tissera et al. (2013). The contribution of endo-debris stars can be diminished by improving the efficiency of the signal-to-noise (S/N) feedback. If the in situ contributions were removed, or diminished significantly, then the stellar mass of these halos would be more in agreement with current MW observations, and the stellar age profiles of some of them would likely be closer to the reported values (see Monachesi et al. 2018 for similar conclusions using the Auriga simulations.)

In order to close the interpretation of the age profiles, we make use of the analysis presented in Tissera et al. (2014), where the mass contributions of satellites with different dynamical masses is investigated. From Figure 6 in that paper, it is clear that Aq-A and Aq-C formed their stellar halos with important contributions from less-massive satellites, while the remaining halos accreted stars from satellites more massive than 10^9.5 M_⊙.¹⁴ This is particularly relevant for the properties of the central regions. Indeed, Aq-A and Aq-C show central age distributions that resemble the ACS observed in the MW. Hence, the analysis of these simulations suggests that the MW should have formed its central regions by the accretion of less-massive satellites ($\lt \sim {10}^{9.5}$ M_⊙), and did not have a significant major-merger contribution. This is consistent with previous works where other methods and models have been adopted (Deason et al. 2017; D'Souza & Bell 2018; Monachesi et al. 2018).

In this Letter we focused on an analysis of the age gradients in the stellar halos of a subset of MW-mass galaxies drawn from the Aquarius Project. We found that these stellar halos exhibit a diversity of age profiles, reflecting their different histories of formation and assembly. Our results suggest that negative age gradients are determined principally by the accreted component, and that the in situ stars affect the slopes in the inner regions, depending upon their relative importance. The flatter age profile of the in situ component is caused by the presence of different well-mixed stellar populations, including the contribution of disk-heated stars and those formed by the gas transported inward by more massive satellites. The negative age gradient set by the accreted stellar component reflects the inside-out assembly of the halo with the contribution of the latest merger events to the outskirts. The characteristics of the accreted satellites such as mass, gas fraction, and their accretion time contribute to modulate the slope, making it less negative if they are able to reach the inner regions carrying younger stars and gas to fuel star formation activity. In general, the in situ component flattens the gradients of the global profiles. The two halos that show the steepest age gradients are those that formed primarily from the contributions of small- and intermediate-mass satellite galaxies. Halos assembled with significant contributions from more massive satellites tend to have shallower age gradients, because these systems carried in younger stars and gas to feed new star formation activity. Our analysis shows that similar slopes of the age gradient reported in the MW's halo can be reproduced by considering only the contribution of the accreted stars. This suggests that the simulated galaxies might be producing an excess of in situ stars, and that the subgrid physics should be improved to reduce such an overproduction of stars. The strong negative age gradient observed in the MW's halo is found in the simulated halos with important contributions from less-massive satellites, suggesting that the MW's halo was assembled from satellites with total dynamical mass lower than ${10}^{9.5}\,{M}_{\odot }$ and a minimal contribution from in situ stars.

We would like to thank the anonymous referee for giving constructive comments that helped to improve this work. P.B.T. acknowledges partial support from Fondecyt Regular 1150334 (CONICYT) and UNAB Project 667/2015. T.C.B. and D.G. acknowledges partial support for this work 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. B.K.G. acknowledges the support of STFC through the University of Hull Consolidated Grant ST/R000840/1, and access to viper, the University of Hull High Performance Computing Facility.