A Tilt in the Dark Matter Halo of the Galaxy

Modern stellar surveys (e.g., Two Micron All Sky Survey, Sloan Digital Sky Survey, Gaia) have carried out a remarkably successful archeological excavation of our Galaxy's fossil graveyard: the stellar halo. Since the foundational works of Eggen et al. (1962) and Searle & Zinn (1978), our picture of the Galactic halo formation history has evolved to encompass both in situ and accretion processes, some of which are continuing to the present day as evident in the Sagittarius stream (Ibata et al. 1994; Majewski et al. 2003) and the infalling Magellanic Clouds (Mathewson et al. 1974; Besla et al. 2007; Zaritsky et al. 2020). Further back in the Galaxy's timeline, a striking discovery from the Gaia mission suggests that the bulk of the stellar halo formed from a single merger event 8–10 billion years ago. Referred to here as Gaia–Sausage–Enceladus (GSE; Belokurov et al. 2018; Helmi et al. 2018), the debris of this merger comprises the majority of the inner stellar halo (Galactocentric radius r_gal ≲ 30 kpc; Iorio & Belokurov 2019; Naidu et al. 2020), providing compelling evidence for hierarchical formation of the Galaxy as predicted by the Λ-CDM standard cosmological model (Press & Schechter 1974; Blumenthal et al. 1984; Bullock & Johnston 2005).

In light of this discovery, two major stellar overdensities in the halo, the Hercules–Aquila Cloud (HAC; Belokurov et al. 2007) and the Virgo Overdensity (VOD; Vivas et al. 2001), have been interpreted as apocentric pileups (Deason et al. 2018; Simion et al. 2019; Naidu et al. 2021) of the GSE merger debris. Supporting evidence of this association is that they occupy diagonally opposing octants in the Galaxy (Iorio & Belokurov 2019), overlap in orbital properties (Simion et al. 2019), and exhibit metallicity distributions indistinguishable from that of GSE (Naidu et al. 2021, hereafter N21). In addition, Han et al. (J. Han et al, in preparation) measure the global shape and orientation of the stellar halo using a sample of giants from the H3 survey (Conroy et al. 2019) equipped with 6D phase space, metallicity, and [α/Fe] abundance measurements beyond distances d ∼ 50 kpc. Using a chemodynamical selection to identify GSE giants, they find that the diffuse halo substructure is best fit by a prolate spheroid, tilted off the Galactic plane pointing toward VOD in the north, and HAC in the south. This is consistent with the distribution of Gaia DR2 RR Lyrae in the halo reported by Iorio & Belokurov (2019).

Parallel to these efforts constraining the spatial extent of GSE, the timing of the merger has also been investigated. Stellar ages estimated from various techniques have been used to infer a timeline of the merger (Gallart et al. 2019; Belokurov et al. 2020; Bonaca et al. 2020; Chaplin et al. 2020; Borre et al. 2021; Grunblatt et al. 2021; Montalbán et al. 2021), building a consensus that the merger completed 8–10 Gyr ago. This merger timeline is also supported by the kinematics of in situ halo stars, which are consistent with being heated from the thick (high α) disk due to GSE (Zolotov et al. 2009; Purcell et al. 2010; Bonaca et al. 2017). In simulations, Fattahi et al. (2019) identify halos in the Auriga project (Grand et al. 2017) that host an accreted, highly radial stellar halo to find that these halos are preferentially created from a single merger with a dwarf galaxy 6–10 Gyr ago. Furthermore, Grand et al. (2018) and Mackereth et al. (2018) demonstrate using Auriga and EAGLE simulations that the bimodality of the [α/Fe]–[Fe/H] sequence in the Galactic disk could be the consequence of an early (z > 1) star formation event triggered by a gas-rich merger. This is once again consistent with the observationally inferred timing of the GSE merger.

While cosmological simulations have been helpful to constrain general properties of the merger, a tailored simulation is required to study the initial conditions of the merger and its subsequent evolution in detail. To this end, N21 present a grid of ∼500 N-body simulations from which they identify a fiducial simulation that closely resembles GSE based on H3 data (Conroy et al. 2019). This simulation is a 2.5:1 total mass ratio merger of the Milky Way (MW) at z = 2 with a dwarf galaxy on an inclined, retrograde orbit. The resulting z-component of the angular momentum L_z and Galactocentric distance r_gal distributions of the merger match those of the GSE sample from the H3 survey (Naidu et al. 2020). Remarkably, the simulation also reproduces HAC/VOD-like stellar overdensities in the correct Galactocentric locations, which was not a condition used to select the GSE analog. Furthermore, these HAC/VOD analogs persist over several gigayears, which is consistent with the observationally inferred timeline of the GSE merger. Lastly, the dark matter contribution from GSE is ∼20% within r_gal < 30 kpc in this simulation, which is consistent with hydrodynamical simulations of GSE-like major mergers in MW-like galaxies (Fattahi et al. 2019; Dillamore et al. 2021).

Thus, a coherent picture emerges: the stellar halo out to 30 kpc is dominated by debris from an ancient (∼8–10 Gyr) radial merger and is globally asymmetric, bookended by the two diagonal apocenters at VOD (above the disk) and HAC (below the disk).

Alternatives to this picture have been proposed. For example, Balbinot & Helmi (2021) integrate the orbits of a local (heliocentric distance d < 2.5 kpc) Gaia EDR3 6D sample to predict the global distribution of halo substructures. After integrating GSE stars over 8 Gyr, they report a symmetric distribution of orbits on the sky. From this, they argue that the current distribution of GSE cannot be globally asymmetric, and thus the observed stellar halo asymmetries in Iorio & Belokurov (2019) could be an artifact of observational biases arising from the Gaia scanning pattern. Also based on orbital integration, Donlon et al. (2019, 2020, 2021) conclude that the radial debris observed in the stellar halo must have been deposited recently (∼ 2 Gyr) in order to explain the spatial coherence of VOD. According to these authors, a large-scale asymmetry in the stellar halo phase mixes too quickly to be compatible with the proposed merger epoch of GSE.

Both of these studies assume a spherical dark matter (DM) halo to integrate orbits of halo stars. Indeed, that is common practice: although observational constraints on the shape of the DM halo are highly uncertain (e.g., Read 2014; Bland-Hawthorn & Gerhard 2016; de Salas & Widmark 2021), a spherical DM halo serves as the default model for the purpose of orbital integration. More sophisticated models still consider oblate (Vera-Ciro & Helmi 2013; Hattori et al. 2021), prolate (Bowden et al. 2016), or triaxial halos (Law & Majewski 2010) that are aligned with the Galactic disk. However, there are well-grounded motivations to consider a DM halo that is tilted with respect to the disk. First, observations of the stellar halo suggest a large-scale tilt on the relevant scales of r_gal ∼ 30 kpc, hinting that the DM halo may embody a similar tilt. Second, cosmological simulations show that it is common for DM halos in MW-like galaxies to exhibit nonspherical shapes that are misaligned with the disk (e.g., Prada et al. 2019; Emami et al. 2021; Dillamore et al. 2021) for a variety of reasons linked to their accretion history. In the Galaxy, such a disk–halo misalignment has been proposed to explain the stability of the Sagittarius stream (Debattista et al. 2013) or the coplanar orbits of Galaxy's satellites (Shao et al. 2021).

Inspired by these clues, in this paper we study the evolution of radial stellar debris in a spherical halo and a tilted halo, and thus the plausibility of unmixed merger debris such as HAC/VOD persisting over several gigayears in the Galactic halo. The rest of the paper is organized as follows. In Section 2 we define the Galactic potentials of a spherical and tilted DM halo. In Section 3 we present the results of integrating stellar halo orbits in these potentials. Finally, in Section 4 we interpret these results in the context of the Galaxy.

In this section, we construct the Galactic potentials used to study the evolution of radial stellar debris in the halo.

We begin by defining the spherical and tilted DM halo models. The former is the canonical picture of the Galactic halo and has been explored to a great extent. For this study, we choose the DM halo described by Bovy (2015), following an Navarro–Frenk–White (NFW) profile (Navarro et al. 1997) with scale radius r_s = 16 kpc and virial mass M_vir = 8 × 10¹¹ M_⊙. The purpose of a tilted halo model in this study is to explore the evolution of large-scale stellar halo asymmetries like HAC/VOD, which N21 successfully reproduces. Thus, we adopt the N21 DM distribution as the fiducial tilted halo model.

Figure 1 shows the DM distribution in N21. Each panel is a 2D histogram of the positions of DM particles, adopting a greyscale that is normalized to include 98% of the sample. Overplotted in solid lines are equipotential contours on a logarithmic scale, and dotted lines mark corresponding contours of the spherical halo from Bovy (2015) to guide the eye. We show the simulated GSE DM distribution at t_lb = 5 Gyr in the top panels, and at present day in the middle panels. These two rows demonstrate that the GSE DM distribution is tilted and stable over the last 5 Gyr. The bottom panels show the combined DM distribution of the GSE and MW analogs at present day. By fitting a triaxial ellipsoid to these distributions, we find that the GSE DM is triaxial (axis ratio 1:0.9:0.57) and the MW DM is oblate (axis ratio 1:1:0.9), resulting in an overall oblate (axis ratio 1:0.96:0.8) DM distribution that is tilted ∼30° above the Galactic disk toward the VOD in the north and HAC in the south.

Zoom In Zoom Out Reset image size

To study the evolution of stellar orbits in these halos, we now calculate the Galactic potential by combining models of the halo, disk, and bulge.

The potential generated by the spherical halo can be analytically derived from the NFW profile (e.g., Łokas & Mamon 2001). Meanwhile, the tilted DM halo does not follow any particular analytical density function, so we employ the self-consistent field method (SCF; Hernquist & Ostriker 1992; Lowing et al. 2011) implemented in the python package gala (Price-Whelan 2017; Price-Whelan et al. 2020) to approximate the potential. The advantage of this framework is its flexibility to approximate nonspherical distributions to the desired level of accuracy. After analyzing the level of variance in the terms of the SCF expansion, we include expansion terms up to n = 20 and l = 4 to sufficiently reproduce the potential while keeping computational costs manageable. The amplitude of the SCF expansion is determined by matching the energies of orbits near the solar neighborhood to that derived from the spherical potential.

The Galactic disk is another important component of the Galactic potential. While N21 self-consistently describes the dynamics of the merging Galaxy at z = 2, its disk model mimics the thick disk of the Galaxy and stays at a fixed mass. In the real Galaxy, the thin (low α) disk is ∼6 times more massive than the thick (high α) disk at the present day (Bland-Hawthorn & Gerhard 2016). Observations suggest that the bulk of the thin disk grew steadily after z ∼ 1–2 (e.g., Haywood et al. 2013; Fantin et al. 2019; Bonaca et al. 2020), following a rapid star formation and quenching history of the thick disk at z ∼ 2. Over 5 Gyr, an average GSE star can cross the disk up to ∼10 times, so the significant growth of the thin disk would clearly impact the dynamics of halo stars. Thus, to examine the effect of the growing thin disk, we consider two models. We first consider a static thick disk with fixed mass 6 × 10⁹ M_⊙, scale length 2 kpc, and scale height 900 pc. This model is analogous to the disk in N21 and the present day thick disk of the Galaxy. We then consider the thick disk coupled with a linearly growing thin disk with scale length 2.6 kpc and scale height 300 pc as described in Bland-Hawthorn & Gerhard (2016). We increase the mass of the thin disk at every megayear, starting from zero and ending at 3.5 × 10¹⁰ M_⊙. This model aims to capture the growth of the thin disk in the Galaxy. Throughout this paper, we refer to the static thick disk model as the "light disk" and the thick and thin disk model as the "growing disk."

Finally, for the bulge, we adopt a Hernquist profile (Hernquist 1990) with mass of 1.4 × 10¹⁰ M_⊙ and scale length of 1.5 kpc (Bland-Hawthorn & Gerhard 2016).

We thus consider four Galactic potentials in this study: a tilted/spherical halo with a light disk, and a tilted/spherical halo with a growing disk. In these potentials, we integrate the orbits of stars from the t_lb = 5 Gyr snapshot in N21 to the present day. Figure 2 shows the N-body snapshots of the stellar halo at these two epochs. We integrate orbits using the Dormand & Prince (1978) integrator at 1 Myr time step, implemented in the python package gala (Price-Whelan 2017; Price-Whelan et al. 2020).

Zoom In Zoom Out Reset image size

Here we present the results from integrating stellar halo orbits in the four Galactic potentials. We outline results separately for the global stellar halo and the local (solar neighborhood) halo, the latter being relevant to studies that utilize local 6D data (e.g., the commonly used Gaia DR2 and EDR3 Radial Velocity Spectrometer (RVS) samples).

3.1. Global Halo

We show the evolution of the global distribution of halo stars in Figure 3. The leftmost panel shows the t_lb = 5 Gyr snapshot in N21, followed by two columns showing the t_lb = 4 Gyr and present day distributions. In the top row we show the N-body simulations. In the middle and bottom rows, we show the positions of stars derived from orbital integrations at their respective t_lb and the preceding 100 Myr. We choose to display the Galactic YZ plane, and note that the XZ and XY planes show the same trend.

Zoom In Zoom Out Reset image size

Figure 3 showcases how the shape of the halo can greatly influence the evolution of stellar orbits, and we highlight two important features. First, the N-body snapshots and tilted halo orbits closely resemble each other, demonstrating the sufficiency of the SCF approximation of the tilted halo potential. Second, the spherical halo erases the tilt in the stellar halo within the first Gyr, whereas the tilted halo preserves the tilt over 5 Gyr. This is somewhat to be expected, since spherical potentials host rosette orbits that can fill the entire available space, while triaxial potentials yield more complex orbits.

Next, Figure 4 shows the effect of a growing disk on stellar halo orbits. We show the present day stellar distribution from orbital integration in a tilted halo and light disk in the top panels, and tilted halo and growing disk in the bottom panels. Most notably, the tilt in the stellar halo with respect to the disk is preserved in both potentials. As expected, orbits in the XZ and YZ planes are somewhat concentrated onto the disk, but the spatial asymmetry is still visibly clear. This indicates that the growth of the thin disk in the Galaxy would not erase HAC/VOD-like overdensities in the halo.

Zoom In Zoom Out Reset image size

Figure 5 further quantifies the degree to which the stellar halo is spatially mixed. We compute the "diagonal excess" of the stellar halo at each time step, defined as

$$\begin{eqnarray}&&\mathrm{Diagonal}\ \mathrm{Excess}\equiv \frac{{f}_{+{\rm{X}}+{\rm{Y}}-{\rm{Z}}}+{f}_{-{\rm{X}}-{\rm{Y}}+{\rm{Z}}}}{2}-\frac{{f}_{\mathrm{other}}}{6},\end{eqnarray} \tag{ 1 }$$

where f_+X+Y−Z and f_−X−Y+Z are the fraction of stars contained in two diagonal octants (respectively containing the bulk of HAC/VOD in the Galaxy) compared to the fraction of stars in the remaining six octants (f_other). For example, a stellar halo that is tilted and aligned with HAC/VOD will show a positive diagonal excess, while a spherical stellar halo will show a diagonal excess of zero. Thus, if the GSE debris becomes fully mixed, the diagonal excess will converge to zero, while if a large-scale asymmetry persists, the diagonal excess will be positive. Figure 5 clearly demonstrates this trend. While the N-body simulation (black line) and tilted halo models (red lines) show mixing timescales greater than 5 Gyr, the spherical halo (blue lines) fully mixes its orbits by the first Gyr. The growing disk (dashed lines) acts to catalyze the spatial mixing of orbits, but not strongly enough to erase global stellar halo asymmetries evolving in a tilted DM halo. We note that when we integrate orbits in a spherical halo only (excluding the disk and the bulge), there is a remaining diagonal excess of ∼0.05 by 5 Gyr. This demonstrates that while the sphericity of the potential is the driving factor in preserving or destroying the large-scale asymmetry of GSE, the disk also plays a role in breaking the spherical symmetry to catalyze the mixing of orbits.

Zoom In Zoom Out Reset image size

Lastly, Figure 6 provides insight into why GSE could be a sensitive probe of the underlying potential. We plot the orbital eccentricity of each star in N21 against ∣Φ_spherical − Φ_tilted∣, where the latter is defined as

$$\begin{eqnarray}&&| {{\rm{\Phi }}}_{\mathrm{spherical}}-{{\rm{\Phi }}}_{\mathrm{tilted}}| \equiv \left|\displaystyle \sum _{i=0\mathrm{Myr}}^{600\mathrm{Myr}}{\phi }_{\mathrm{spherical},i}-{\phi }_{\mathrm{tilted},i}\right|,\end{eqnarray} \tag{ 2 }$$

where ϕ_spherical,i and ϕ_tilted,i are the spherical/tilted potentials evaluated at a star's location at i Myr. The 600 Myr summation timescale is chosen such that most stars would have had enough time for at least one pericenter passage. The eccentricity is calculated in the tilted potential, and we make note that the result does not change when eccentricity is calculated in a spherical potential. There is a strong positive correlation in eccentricity and ∣Φ_spherical − Φ_tilted∣, and the stars on disky orbits show ∣Φ_spherical − Φ_tilted∣ converging to zero. This correlation indicates that eccentric orbits that probe the full range of radii in the halo are more sensitive to the underlying potential than disky orbits. Since GSE stars are predominantly on highly radial/eccentric orbits, it is a sensitive tracer of the shape of the DM halo that its progenitor deposited in the Galaxy.

Zoom In Zoom Out Reset image size

3.2. Local Halo

We define the local (solar neighborhood) halo to be a sphere of d ≲ 3 kpc centered on the Sun. The local halo sample has been historically important as it contains the most complete 6D information of stars. Numerous studies have extrapolated the motions of local halo stars in order to understand the global halo (e.g., Carollo et al. 2010), and in this section we assess how well the local halo sample represents properties of the global halo.

In Figure 7 we identify local halo stars in N21 that satisfy d < 3 kpc at the present day (top row, marked in red), and track where they were 500 Myr ago (middle bottom rows). If the local halo is an unbiased sample of the global halo, then we expect these two distributions to look similar. However, in the XY and XZ projections, the t_lb = 500 Myr distribution of the global and local halo are notably different. This suggests that the orbits of local halo stars do not reflect the spatial distribution of the global halo, but are rather on more spatially mixed orbits. This is not surprising, since the solar location is off-axis from the tilt of the global halo and off-centered from the Galactic center.

Zoom In Zoom Out Reset image size

Examining the YZ projection in Figure 7, we see that the bottom two panels are similar to each other. This suggests that the local halo stars retain memory of the global halo in the YZ plane, since the solar location on the YZ plane is not off-centered with respect to the Galactic center. The exact relationship between the local and global halo is a complicated product of the spatial/dynamical distribution of the global halo and the solar position, and cannot be decoupled from one's assumptions about the global halo.

In this study we have investigated the evolution of stellar debris from an inclined merger remnant in the Galactic halo. Starting with a simulated stellar halo from Naidu et al. (2021), we integrated stellar orbits in tilted and spherical DM halos with both light and growing disks to present day. After 5 Gyr, orbits in a spherical halo are fully phase mixed, while orbits in a tilted halo preserve a large-scale tilt even in the influence of a growing (massive) disk. We further quantified stellar halo asymmetries by calculating the excess fraction of stars in the two octants containing HAC/VOD, from which we evaluated the timescales of spatial mixing in the Galactic halo. Finally, we examined the orbits of local halo stars, and showed that these orbits do not trace the spatial distribution of the global halo.

We now interpret these results in the context of the Galactic halo. The interpretation of HAC/VOD as apocenter pileups from GSE is supported by compelling evidence from spatial, dynamical, and chemical considerations. A key issue is whether one would expect HAC and VOD to be long-lived structures, in light of the evidence that the GSE merger occurred ∼8–10 Gyr ago. If the DM halo of the Galaxy is spherical (or in general is aligned with the disk), then any large-scale asymmetries traced by radial merger debris will be quickly erased. In this case it would be difficult to understand how HAC and VOD could be associated with an ancient merger. However, a tilt in the DM halo, as we have explored in this work, enables these overdensities to persist for many gigayears. The evidence supporting an association between HAC/VOD and GSE, combined with the estimated age of the merger, leads us to conclude that the DM halo of the Galaxy is likely tilted with respect to the disk in the direction of the HAC and VOD overdensities.

In a globally tilted halo, the local halo occupies a spatially and dynamically biased location due to its offset from the Galactic center. In particular, any sample of stars on radial orbits in the local halo will preferentially contain orbits that have been phase mixed in the Galaxy. In certain Galactocentric projections (e.g., YZ) where the solar neighborhood is not offset from the Galactic center, these orbits can still retain information about the global halo. In other projections, most of this information is lost. Inferring properties of the global halo by extrapolating the local halo must therefore be carried out with caution.

The nontrivial shape of DM halos has been the topic of many recent studies using cosmological simulations. For example, Emami et al. (2021) report that only 32% of Milky Way-like galaxies in the Illustris TNG50 simulation (Nelson et al. 2019; Pillepich et al. 2019) show simple DM halos, and the rest of the sample exhibit either twisted (gradually rotating) or stretched (abruptly rotating) halos. Prada et al. (2019) find that 20% of Milky Way-like galaxies in the Auriga project (Grand et al. 2017; Monachesi et al. 2019) have twisted halos, while Dillamore et al. (2021) find that most galaxies in the ARTEMIS simulations (Font et al. 2020) that show GSE-like features exhibit a global change in orientation of their DM halo. In EAGLE simulations (Crain et al. 2015; Schaye et al. 2015), Shao et al. (2021) find that the DM halo of galaxies are well aligned with the common orbital plane of their satellite galaxies. They notice that the common orbital plane of the Galaxy's satellites is misaligned with the disk, from which they argue that the DM halo of the Galaxy must be twisted. In isolated simulations, Debattista et al. (2013) predicted a halo–disk misalignment in the Galaxy based on models of the Sagittarius stream. Thus, the halo–disk misalignment that we find in this study is consistent with theoretical expectations.

An important caveat of this study is that we do not consider the mutual interaction between the halo and its environment: the disk and the satellites of the Galaxy. For example, the disk will pull the tilt of the halo closer to the plane, while the halo may induce a warp in the disk at large radii. This could mean that the initial tilt of the halo has to be larger than the current day configuration to explain the spatial distribution of HAC/VOD. We note that GSE contributes only ∼10% of the dark matter within r_gal < 10 kpc in the N21 simulation, so where the disk gravity is dominant the contribution from GSE is smaller than in the outer halo. Among the satellites of the Galaxy, the Large Magellanic Cloud (LMC) has a significant influence on the global shape of the halo (e.g., Garavito-Camargo et al. 2019; Erkal et al. 2020; Conroy et al. 2021; Petersen & Peñarrubia 2021). However, this effect is most prominent at larger distances (r_gal > 50 kpc) than what we consider here, and within r_gal < 30 kpc the influence of the LMC is likely to be much smaller compared to GSE (though see Lucchini et al. 2021 for an alternative scenario). Smaller satellites also contribute to the steady growth in the mass of the halo, but this accretion again occurs mainly in the outer halo while the inner halo is expected to remain quiescent after t_lb ∼ 8 Gyr (e.g., Bullock & Johnston 2005; Font et al. 2006; De Lucia & Helmi 2008; Cooper et al. 2010; Pillepich et al. 2014; Bland-Hawthorn & Gerhard 2016).

Future work will address these issues. Additional observations are necessary to constrain the precise locations and extent of HAC and VOD and the overall tilt of the stellar halo, which will inform the configuration of the DM halo. In light of ongoing and upcoming spectroscopic surveys such as H3, DESI, SDSS-V, WEAVE, and 4MOST, these are timely endeavors.