The Stellar Halo of the Galaxy is Tilted and Doubly Broken

We model each star as drawn from an inhomogeneous Poisson point process. This likelihood-based method to constrain density profiles was pioneered by studies such as Bovy et al. (2012), Rix & Bovy (2013), and Xue et al. (2015). We define the fundamental variables that determine the Poisson state of a star to be the 3D location r , initial mass m, and metallicity [Fe/H]. We express r in either galactocentric coordinates (X_gal, Y_gal, Z_gal) or heliocentric coordinates (l, b, d). For the galactocentric coordinates, we adopt a right-hand coordinate system with the solar location at (−8.122 kpc, 0 kpc, 0.02 kpc) and the Galactic center at the origin. While stellar age is another fundamental variable, most GSE stars are within a stellar age range of 8–10 Gyr (e.g., Gallart et al. 2019; Bonaca et al. 2020). Hence, stellar ages do not significantly affect the observable parameters, and we assume a uniform stellar age of 10 Gyr for this study. The intrinsic distribution of the state variables, the rate function λ, is assumed to be separable with the initial mass function (IMF) and the metallicity distribution function (MDF) as follows:

$$\begin{eqnarray}&&\lambda =\rho ({\boldsymbol{r}})\cdot \mathrm{IMF}(m)\cdot \mathrm{MDF}([\mathrm{Fe}/{\rm{H}}]).\end{eqnarray} \tag{ 1 }$$

The assumption of separability is equivalent to assuming that the MDF and IMF do not depend on Galactic location. Because we select a chemically homogeneous population, the MDF to first order does not depend on location. In detail, an imprint of the intrinsic metallicity gradient of the GSE progenitor could be left behind from different stripping episodes (i.e., the outermost stars that are likely more metal-poor are stripped first), but Conroy et al. (2019a) and N21 show that this gradient over the range 6 < r_Gal < 60 kpc is essentially zero. Thus, the assumption of a separable MDF is a good approximation. We adopt a global IMF (Kroupa 2001) and an empirical MDF derived from a fitted accretion model (see Conroy et al. 2019a for details). Folding in the H3 selection function, the observed rate can be written as

$$\begin{eqnarray}&&{\lambda }_{{obs}}=\lambda ({\boldsymbol{r}},m,\mathrm{Fe}/{\rm{H}})\cdot S(\mathrm{tile}),\end{eqnarray} \tag{ 2 }$$

where S is a tile-dependent selection function on the observed stellar photometry and Gaia parallax π. To transform the fundamental variables m and [Fe/H] into photometry, we use MIST v1.2 isochrones (Dotter 2016; Choi et al. 2016) at 10 Gyr and place them at a heliocentric distance d derived from r . As discussed earlier, the stars in the GSE sample show stellar ages between 8−10 Gyr: at those ages the changes in stellar isochrones are minute, and the specific choice of isochrone age between 8−10 Gyr does not introduce significant systematics in either the selection function or the rate function. The H3 target selection involves a cut on r-band magnitude, and the limiting magnitude of our giant selection for a given tile further depends on factors such as sky brightness and weather. A histogram of limiting magnitude for S/N = 5 is presented in Figure 14 in the Appendix, and we take this into account in our selection function. While S should intuitively be a binary function (a star is either observable or not observable), it is actually a continuous function from 0 to 1. First, Gaia parallaxes have a statistical uncertainty that causes stars at the edge of detectability to have nonnegligible scatter. We account for this effect by marginalizing the rate function over the Gaia parallax uncertainty, which is a function of the ecliptic angle, G-band magnitude, and V − I color. Second, we can only place a maximum of 200 fibers within any given field, sometimes even fewer than 200 to avoid fiber collisions among the robotic positioners. If there are more targets than fibers available, we must include the sampling fraction (defined as N_{received fiber}/N_target) in the observed rate function. This sampling fraction is calculated based on the parent Pan-STARRS catalog, and we include it in the rate function as a constant multiplicative factor for each tile. Finally, H3 implements a ranked targeting scheme, where rank 1 stars (BHB, RRL, photometrically selected luminous K giants) are given priority over rank 2 stars (selected on parallax and r magnitude). Thus, the sampling fraction is also a function of the rank of the target. We note that rank 1 stars only comprise 6% of the sample. Accounting for the ranks, the observed rate function can be written as

$$\begin{eqnarray}&&{\lambda }_{{obs}}=\lambda ({\boldsymbol{r}},m,[\mathrm{Fe}/{\rm{H}}])\cdot \left[{ \mathcal B }\cdot {f}_{{r}_{1}}+{ \mathcal B }\cdot {f}_{{r}_{2}}\right],\end{eqnarray} \tag{ 3 }$$

where ${ \mathcal B }$ is a Boolean function on the state variables that determine if the star is rank 1, rank 2, or nonobservable; ${f}_{{r}_{1}}$ and ${f}_{{r}_{2}}$ are sampling fractions of the respective target ranks. The definition of rank has evolved over the duration of the survey, and we report these changes in the Appendix. We note that there also exist rank 3 (faint filler targets) and rank 4 (nearby filler targets), but these targets comprise less than 1% of the giants with S/N > 5 and we exclude them from our analysis.

Once we have specified the observed rate function, we can evaluate the likelihood function for an inhomogeneous Poisson process:

$$\begin{eqnarray}&&\mathrm{log}{ \mathcal L }=-{N}_{\lambda }+\sum _{i}\mathrm{log}{\lambda }_{{obs}}({{ \mathcal D }}_{i}),\end{eqnarray} \tag{ 4 }$$

where ${{ \mathcal D }}_{i}\equiv ({{\boldsymbol{r}}}_{i},{m}_{i},{[\mathrm{Fe}/{\rm{H}}]}_{i})$, and

$$\begin{eqnarray}&&{N}_{\lambda }=\int {\lambda }_{{obs}}\ {\rm{d}}{\boldsymbol{r}}\ {\rm{d}}m\ {\rm{d}}[\mathrm{Fe}/{\rm{H}}].\end{eqnarray} \tag{ 5 }$$

Intuitively, N_λ represents the expected number of events (stars) observed through the selection function given a set of model parameters. While this integral would be computationally formidable to perform over the full 5D space (three spatial dimensions, m, and metallicity), it is tractable if we approximate all of the stars in one tile to lie along one line of sight at the center of the tile, and thus collapse three spatial dimensions to one. Since the size of each tile (05 radius) is much smaller than the total footprint (∼15,000 square degrees), we confidently make this approximation.

We now discuss the functional form of ρ( r ). Assuming that the halo can be approximated by an ellipsoid (which includes spherical and flattened spherical shapes), there are two parts to consider in parameterizing ρ( r ). The first is the radial profile, representing the density along the principal axes of the ellipsoid. Previous studies have fit a singly broken power law or a Sérsic profile to stellar halo tracers. Here, we consider a multiply broken power law that encompasses a no-break to doubly broken power law.

The second component of ρ( r ) is the angular profile, representing the overall orientation of the ellipsoid. This profile has historically not been explored as in depth as the radial profile. However, there is no a priori reason to assume that a merger remnant would occupy an azimuthally symmetric distribution in the Galaxy; indeed, there is increasing evidence to consider a halo that is allowed to be rotated with respect to the Galactic plane (Iorio & Belokurov 2019; Naidu et al. 2021b; Han et al. 2022). In this study, we consider the simplest form of the angular profile, which is to assume the ellipsoid as a solid body that has been rotated. More flexible models could consider the rotation angles changing as a function of radius, and we leave this to future work. Combining the radial and angular profiles, we write the 3D density profile as follows:

$$\begin{eqnarray}&&\rho ({\boldsymbol{r}})=f(M,{\alpha }_{1},{\alpha }_{2},{\alpha }_{3},{r}_{b,1},{r}_{b,2},p,q,\mathrm{pitch},\mathrm{yaw}).\end{eqnarray} \tag{ 6 }$$

Here, f is a multiply broken power law that is described by 10 parameters: total stellar mass enclosed in Galactic radius 1–100 kpc, M, first power-law slope α₁, second power-law slope α₂, third power-law slope α₃, first breaking radius r_b1, second breaking radius r_b2, intermediate axis flattening p, minor axis flattening q, pitch angle, and yaw angle (all angles measured clockwise). The power law is a function of the flattened radius ${r}_{q}\equiv \sqrt{{X}^{2}+{(Y/p)}^{2}+{(Z/q)}^{2}}$, where p and q are values between 0 and 1. The X, Y, and Z axes are, respectively, defined to the major, intermediate, and minor axes. The power-law slopes and breaking radii are allowed to overlap such that f can describe a single power law to a triple power law (i.e., doubly broken). The pitch angle represents the angle of the major axis of the ellipsoid with respect to the Galactic plane (colloquially referred to as tilt), and the yaw angle represents the azimuthal angle of the major axis with respect to the Galactic Z-axis. As we discuss in Section 4, the best-fit ellipsoid model is near-prolate, and thus has approximate internal azimuthal symmetry. Hence, we do not report the rotation of the ellipsoid about the long axis, i.e., roll angle. When we apply a fit of all three Euler angles, the shape and flattening parameters show an insignificant change, and the pitch and yaw angles change within the statistical uncertainties of the model.

Finally, to write the density function in terms of observed coordinates (l, b, d), we must account for the change in the volume element as a function of observed coordinates. This is the absolute value of the Jacobian determinant, and can be derived as $| | J| | (l,b,d| {\boldsymbol{r}})={d}^{2}\cos (b)$.

The task at hand of computing the likelihood function can naturally be parallelized over tiles. We use OpenMPI (Gabriel et al. 2004) to distribute the likelihood function over 1000 CPUs, enabling on-the-fly calculations of the likelihood instead of having to rely on precomputed grids. We assume a uniform prior for all parameters, and use emcee (Foreman-Mackey et al. 2013) to sample from the posterior. We employ 200 walkers and a nominal burn-in phase of 200 steps. The number of final steps vary across models, but are typically on the order of a few hundred.

The forward model framework described above is also a generative model with which we can create mock data sets. Given true model parameters, we can sample stars from the spatial density, IMF, and MDF. We can then pass these stars through the full H3 selection function to create a mock data set. Such mock data sets are useful for several purposes. First, we can fit the mock data set through our full fitting pipeline to validate our inference framework. An example of such a fit is displayed in Figure 4. The diagonal panels show the posterior of each parameter as sampled with emcee, while the other panels show the covariance among different parameters. The blue lines indicate the sample mode (also written as text on top of the diagonal panels), and the red lines indicate the true parameters that were inputs to the mock data. All of the true parameters are well within the statistical uncertainties of the inferred parameters, demonstrating the robustness of our inference framework.

Zoom In Zoom Out Reset image size

Second, we use the mock catalog to understand how model assumptions can affect our interpretation of the data. For example, we create a tilted, triaxial, doubly broken mock data set (from the same parameters as Figure 4), and then fit it with a nontilted, spheroidal (one flattening parameter q), singly broken power law. This finding is analogous to previous measurements of the stellar halo. Figure 5 shows the posterior distribution of such a fit. We see that the model has been strongly constrained to an oblate spheroid (q = 0.84) with a single breaking radius in between the two true breaking radii. This result is a cautionary tale: despite the posterior being strongly constrained to be an oblate spheroid with a single break, the underlying true density can be triaxial and tilted. Thus, there is strong motivation to project the best-fit models back onto data space, where clear systematics (such as misinterpreting a triaxial tilted ellipsoid as an oblate spheroid) would be visible as residuals between the data and the model. This lesson forms the basis of how we examine and present our results in the following section.

Zoom In Zoom Out Reset image size

Our fiducial model for the stellar halo is a triaxial, doubly broken power law with an orientation not constrained to lie in the Galactic plane. There are three key features to our best-fit model: (1) we find two breaking radii at 12 kpc and 28 kpc; (2) the principal axes ratios are 1:p:q = 10:8:7, which is nearly prolate; (3) we find that the ellipsoid is tilted by 25° off of the Galactic plane and 24° away from the Sun–Galactic center axis. Figure 6 shows the full posterior of each parameter of the fiducial model. That the covariance among parameters is mostly small indicates that our model parameters are not degenerate with one another. This is consistent with our expectation that each parameter should affect distinct spatial features of the ellipsoid. In Figure 7 we plot the isodensity surface of the fiducial model in galactocentric coordinates, evaluated at r_q = 20 kpc. The solar location is marked as a red star, and an artist's impression of the Galactic disk (NASA/JPL-Caltech/R. Hurt 2008) is overplotted to orient the reader.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 8 we show the on-sky stellar density of the data and of the models projected onto data space. The top panel shows the density of GSE, and the middle left (right) shows the density of a nontilted (fiducial) model. Densities are normalized by the mean density in each panel. The bottom-left (right) panel shows the residuals of the data subtracted by the nontilted (fiducial) model, normalized by the mean density of the data. The fiducial model is a tilted, triaxial ellipsoid, and the nontilted model is created by fitting a model with rotation angles (both pitch and yaw) fixed to zero. All panels are smoothed with a 20° Gaussian kernel. The faint white dots mark H3 tiles. The top and middle panels reflect the H3 selection function to first order. However, since the data and the models are subject to the same selection function, the difference in bottom-panel model residuals are insensitive to the selection function. In this context, we observe a striking dipolar over-underdensity pair in (−l, +b) and (+l, −b) for the nontilted model residual plot. This dipole is bookended by the Virgo Overdensity in the north and the Hercules-Aquila Cloud in the south. The dipole residual is greatly reduced in the tilted (fiducial) model, demonstrating that an overall rotation to the ellipsoid can explain the observed dipole in the nontilted residual. However, the residuals of the fiducial model is still imperfect, which we discuss in Section 5.

Zoom In Zoom Out Reset image size

In Figure 9 (top panel) we plot a log-bin histogram of observed flattened radii of GSE (black dots) and compare to the fiducial model (red line) as well as best-fit singly broken (green line) and no-break (blue line) power-law models. In the bottom panel we plot the residuals of each model to the data, and indicate the Poisson uncertainty of the data in red shading, evaluated from 100 mock catalogs. While the single-break power law fits the r_q > 15 kpc halo well, it overpredicts the stellar density at inner radii. To model both the inner and outer halo, the data requires a minimum of two breaking radii. Using the Akaike information criterion (AIC; Akaike 1974), we can use the likelihood values to determine which model better describes the data. AIC is defined as $2k-2\mathrm{ln}{ \mathcal L }$, where k is the number of free parameters and ${ \mathcal L }$ is the maximum likelihood value for the model. Based on this criterion, the probability that the singly broken power law minimizes the information loss compared to the doubly broken power law is e⁻²³⁷²⁴, essentially zero. We thus confidently rule out the singly broken power law as the preferred model.

Zoom In Zoom Out Reset image size

It is worth nothing that the r < 10 kpc halo strongly overlaps with in situ stars (both the disk and in situ halo), and many prior studies of the stellar halo therefore exclude this region from their sample. Indeed, the inner breaking radius would be difficult to observe if we were not able to separate GSE from the in situ stars. For example, the RRL sample used by Iorio et al. (2018) that spans the inner halo down to r = 1.5 kpc would include many in situ halo and disk stars, which could be a reason why they did not measure a significant tilt in the RRL data. Owing to the clean chemistry and kinematic selection enabled by H3, we are able to probe the density profile to relatively small galactocentric radii. Furthermore, if our sample were significantly contaminated by disk stars, the innermost power-law slope would likely be steeper and the difference between α₁ and α₂ would be smaller than what we measure. Hence, our measurement of the innermost breaking radius is likely not a product of disk contamination. Lastly, α₁ and α₂ are ∼8σ apart from each other and α₂ and α₃ are ∼15σ apart, demonstrating the statistical significance of the doubly broken power law.

In this paper we have explored a flexible model (triaxial ellipsoid that is allowed to rotate) and new data (a high-purity GSE sample obtained from chemo-dynamical selection) to measure the density profile of the most dominant contributor to the stellar halo. It is important to understand how our results depend on each component.

From Figure 5, we have learned that an intrinsically triaxial and tilted ellipsoid can be incorrectly constrained to be an oblate spheroid if the model is restricted to be aligned with the Galactic plane. Thus, if we fit our fiducial GSE sample with a model that is restricted to be aligned with the plane, we would obtain an oblate spheroid as our best-fit model. However, we know from Figures 8 and 9 that this oblate spheroid model does not reproduce either the radial or angular distribution of the data. Thus, this exercise once again reveals the limitations of fitting a restricted model to the data, and the value of a generative framework that allows one to evaluate the model in data space.

To connect our data to previous studies, we fit a restricted model to a sample (an oblate spheroid with a singly broken power law) to a halo sample selected only on ∣Z∣ > 4 kpc, with no chemical or kinematic selection. Such a geometric selection of halo stars being above some Galactic height is characteristic of pre-Gaia halo samples (e.g., Xue et al. 2015). Our sample is, of course, still specific to the H3 footprint, so it is not identical to prior studies. Figure 10 shows the posterior of this restricted fit. Similar to previous work, we find that the data can be strongly constrained to an oblate spheroid (q ∼ 0.9) with one breaking radius at ∼20 kpc. The inner power-law slope, α₁ = 3.3, is on the steeper end of literature values, and the outer power-law slope, α₂ = 4.6, is near the average of literature values. We note that when we fit the fiducial model to the ∣Z∣ > 4 kpc sample, the pitch (186 ±,182) and yaw (−125 ± 154) angles are consistent with zero and the triaxiality is reduced (p =0.84 ± 0.7 and q = 0.88 ± 0.04). We further elaborate on placing our results in the context of the literature in Section 5.

Zoom In Zoom Out Reset image size

We now assess the specific choice of eccentricity used to define GSE. To this end, we return to Figure 3, where we show the distribution of eccentricity against metallicity of the parent sample. While GSE dominates the high-e sample, it does extend down to lower eccentricities. It is thus worthwhile to explore how our result changes if we relax the eccentricity selection of GSE to be e > 0.5. This increases our GSE sample size to 7886. The result is summarized in Table 1. Other than the total stellar mass, which increases proportionally to the sample size as expected, none of the density parameters change in any appreciable way. This also demonstrates that our kinematic cut does not induce large selection effects in the radial range that we explore in this study.

Table 1. Summary of Best-fit Parameters

Sample	$\mathrm{log}(M/{M}_{\odot })$	α1	α2	α3	rb,1 (kpc)	rb,2 (kpc)	p	q	Pitch (°)	Yaw (°)
Fiducial Model
Fiducial GSE	${8.76}_{-0.02}^{+0.02}$	${1.70}_{-0.24}^{+0.16}$	${3.09}_{-0.11}^{+0.10}$	${4.58}_{-0.11}^{+0.10}$	${11.85}_{-0.71}^{+0.92}$	${28.33}_{-1.55}^{+1.05}$	${0.81}_{-0.03}^{+0.03}$	${0.73}_{-0.02}^{+0.02}$	$-{25.39}_{-3.20}^{+3.11}$	$-{24.33}_{-5.51}^{+4.94}$
e > 0.5 GSE	${8.88}_{-0.02}^{+0.02}$	${1.70}_{-0.06}^{+0.05}$	${3.10}_{-0.06}^{+0.06}$	${4.65}_{-0.04}^{+0.03}$	${10.61}_{-0.46}^{+0.43}$	${29.26}_{-1.07}^{+0.90}$	${0.82}_{-0.02}^{+0.02}$	${0.74}_{-0.02}^{+0.02}$	$-{25.13}_{-1.97}^{+2.02}$	$-{23.88}_{-3.36}^{+3.36}$
Parent	${9.13}_{-0.03}^{+0.04}$	${2.08}_{-0.20}^{+0.24}$	${3.45}_{-0.09}^{+0.10}$	${4.45}_{-0.16}^{+0.15}$	${9.48}_{-0.91}^{+0.86}$	${27.68}_{-1.67}^{+1.92}$	${0.88}_{-0.02}^{+0.02}$	${0.78}_{-0.01}^{+0.01}$	$-{10.40}_{2.09}^{+2.08}$	$-{11.06}_{-4.15}^{+4.53}$
Restricted Model: single break, p = 1, pitch = 0, yaw = 0
∣Z∣ > 4 kpc	${9.10}_{-0.05}^{+0.06}$	${3.30}_{-0.09}^{+0.10}$	${4.59}_{-0.15}^{+0.26}$	⋯	${21.51}_{-1.25}^{+2.00}$	⋯	1 (fixed)	${0.87}_{-0.03}^{+0.02}$	0 (fixed)	0 (fixed)

Download table as:  ASCIITypeset image

Finally, we explore how our result changes in the complete absence of a chemo-dynamical selection by fitting the whole parent sample with the fiducial model. This sample includes the thick disk and in situ halo, and reflects the total sum of distant (d > 1 kpc) giants in H3 that are not part of a kinematically cold structure (i.e., Sagittarius stream, Sextans, Ursa Minor, Palomar 13, Draco, M31, M33, M92, and GD-1). The result is summarized in Table 1. The tilt angle, while still nonzero, is notably smaller (∼10°) and the shape is closer to oblate (10:9:8), which is what we expect since the thick disk and in situ halo occupy a distribution that is different from GSE and is aligned with the Galactic plane. Additionally, the inner power-law slope is steeper (α₁ ∼ 2.1) than the fiducial model, again likely due to the thick disk and in situ halo as their distributions are steeper than GSE and occupy the inner 10 kpc.

In this study, we have explored a triaxial ellipsoid model to measure the shape and density profile of the most dominant component of the accreted stellar halo of the Galaxy. This model reproduces key features of the data in both radial and angular projections, corresponding to the overall tilt of the halo and the double-break in the power law. We thus believe the fiducial model to be representative of the underlying data. However, the true stellar halo is surely not a perfect ellipsoid as we have drawn in Figure 7. While the on-sky density residuals of the fiducial model is clearly better than that of the nontilted model in Figure 8, we still observe a spatially correlated offset from the data. This is likely due to a more complex shape of the halo, such as its overall orientation changing as a function of radius (e.g., Prada et al. 2019; Emami et al. 2021; Shao et al. 2021 show this in simulations) and/or having different radial density distributions on either side of the Galactic center due to dynamical friction as the progenitor sinks into the Galaxy (Naidu et al. 2021b). While our work allows for a more general shape than in previous work, there are clear motivations to consider even more flexible models. Factoring in the in situ disk and halo, the combined halo is likely a superposition of at least two distinct shapes: a triaxial ellipsoid accreted halo, and an oblate in situ halo. Furthermore, some studies have measured a radially varying flattening parameter (Xue et al. 2015; Das & Binney 2016; Iorio et al. 2018; Iorio & Belokurov 2019), which could partially be due to contamination from the in situ halo and other substructures. Future studies that can employ more data points should be able to simultaneously constrain radially varying rotation angles and shape parameters.

Another limitation of this study is that H3 does not contain data below δ < −20°, so we cannot measure features to the halo profile contained in those regions. Future spectroscopic surveys in the southern hemisphere such as 4MOST and the Sloan Digital Sky Survey (SDSS)-V will provide important constraints on the behavior of the stellar halo. Furthermore, the distances estimated by MINESWeeper have ∼10% statistical uncertainties that we do not directly factor into the inference framework. Instead, we verify that these uncertainties do not significantly affect our results by fitting mock catalogs that include distance uncertainties. All parameters are inferred to be well within statistical uncertainty from their true values, as we show in Figure 4. We have also removed very distant stars beyond r > 60 kpc that are most affected by distance uncertainties and have poor astrometric quality. To properly model the outer halo, where tracers are rare and distance uncertainties can sway the results, one should marginalize the likelihood function over the distance posterior. The constituents of the distant r > 60 kpc halo are largely uncharted, and will be an important direction for future work. Our work also does not attempt to model the spatial-dependent kinematics of GSE, which could introduce systematic errors (e.g., orbits becoming more circular in the outer halo, where unrelaxed tidal debris may still exist). Lastly, while the fiducial e > 0.7 sample should largely be comprised of GSE, the e > 0.5 sample may contain other mildly eccentric, accreted substructures such as Sequoia (Matsuno et al. 2019; Myeong et al. 2019), Arjuna, and I'itoi (N20). However, these structures are subdominant (Sequoia, for example, is estimated to be ∼40 times less massive than GSE). Thus, their contribution to the density profile is very small compared to GSE. Nonetheless, understanding the spatial distribution of these structures is important to completing our knowledge of the accreted stellar halo.

Due to our novel sample comprising high-probability GSE stars, it is not a straightforward task to directly compare our results to previous studies. Instead, here we place our fiducial density model in the context of previous studies of the stellar halo.

We first highlight the double-break feature of our density profile. Figure 11 shows the distribution of breaking radii r_b in the literature compared to this work. Previous results report a single breaking radius, and we visualize these radii as Gaussian distributions based on the reported uncertainties. The breaking radii are in flattened units, so the physical radii can change by a factor determined by the flattening parameter. However, this change is not larger than a few kiloparsecs for the most flattened models, and does not affect the evident dichotomy in the literature values of r_b . Specifically, there exists a cluster of r_b ∼ 15–20 kpc and a cluster of r_b ∼ 25–30 kpc. If we interpret the double-break power law that we measure to be the outer and inner apocentric pile-ups of GSE, we can understand the dichotomy as an artifact of the survey footprint and radial range probed by previous studies. Depending on the specific sight line, the outer (inner) apocentric pile-up can be more prominent than the inner (outer) pile-up due to the tilt and orientation of GSE. It is therefore likely not coincidence that the two breaking radii we measure coincides with the two clusters of results in the literature.

Zoom In Zoom Out Reset image size

In Figure 12 we plot the density profiles from the literature over their respective radial range of the study. In Figure 13 we make a similar comparison for the power-law index, where the transition regions between inner and outer power-law indices correspond to reported uncertainties of the breaking radius. We see that the density profile presented in this study approximately traces the average literature profile in the inner (r < 15 kpc), mid (15 kpc < r < 25 kpc), and outer (r > 25 kpc) halo.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A noteworthy feature is that our inner breaking radius is smaller than literature estimates. Several factors could lead to this result. First, our halo sample includes stars at smaller galactocentric radius than most studies owing to our ability to effectively remove in situ stars. Second, from the mock test shown in Figure 5 we have seen that by approximating a singly broken power law to an intrinsically doubly broken power law, we estimate a breaking radius that is greater than the inner breaking radius. From this we can infer that even in a study that surveys sight lines dominated by the inner breaking radius, stars from the outer breaking radius will shift the measured breaking radius outwards. On the contrary, previous samples dominated by stars at greater distances are not sensitive to the inner breaking radius.

We now discuss the shape and orientation of our model. The vast majority of previous studies report an oblate spheroid that is aligned with the disk, which has principal axes ratios of 1:1:0.6–0.8 (Bland-Hawthorn & Gerhard 2016, and references therein). Among the few studies that have fit a triaxial ellipsoid, Iorio et al. (2018) report axes ratios of 1:0.8:0.5−0.6 from Gaia DR1 RRL cross-matched with 2MASS out to r_Gal < 28 kpc. This axis ratio is comparable to what we find, which is 1:0.8:0.7, as shown in Figure 7. In terms of orientation, most studies assume a fixed orientation of the halo that is aligned with the disk. However, Iorio & Belokurov (2019) report a misalignment (yaw) angle of the major axis with the Galactic X-axis by 70° based on DR2 RRL. The yaw angle reported by Iorio & Belokurov (2019) and our study both point toward the same quadrants in the Galactic X–Y plane, although the precise value of the angle is different (24° from the X-axis in our study). Given that the sample used in this study is significantly different from the Gaia RRL due to our chemo-dynamical selection, this difference in yaw angle may not be significant. Future work should clarify this issue. Furthermore, Iorio & Belokurov (2019) speculate the existence of a tilt of the halo with respect to the Galactic plane based on the X–Z projection of RRL, and the direction of the tilt is consistent with what we report in this paper. In both studies, the ellipsoid is tilted above the plane in the direction of the Sun.

Finally, we discuss our estimate of the total stellar mass of GSE. As summarized in Table 1, we infer a total GSE stellar mass of 5.8−7.6 × 10⁸ M_⊙, depending on if we define GSE to be e > 0.7 or e > 0.5. We note again that we calculate the total mass to be enclosed in 1–100 kpc. As we do not have data within r_Gal < 5 kpc, the mass enclosed in this region is extrapolated from the inferred density profile, accounting for ∼10% of the total mass. The inferred mass range of GSE comfortably brackets the literature values of the GSE mass such as Helmi et al.'s (2018) 6 × 10⁸ M_⊙, Vincenzo et al.'s (2019) 5 × 10⁸ M_⊙, Feuillet et al.'s (2020) 7 × 10⁸ M_⊙, and N20's 4−7 × 10⁸ M_⊙. There exist exceptions, such as Mackereth et al. (2018), who find a lower accreted stellar halo mass of 3 × 10⁸ M_⊙ at e > 0.7. We note that the restricted fit we present in Figure 10 yields a significantly higher mass estimate of 1.3 × 10⁹ M_⊙ than the fiducial fit, which shows that the triaxiality and the existence of an inner break (thus allowing for a shallower innermost power-law slope) can significantly affect the total stellar mass estimate. This higher mass estimate is remarkably similar to what Deason et al. (2019) find for the total stellar halo mass based on RGB stellar counts, 1.4 × 10⁹ M_⊙.

A tilted, doubly broken density profile of the stellar halo has a number of implications for the galaxy. First, the overall tilt of the stellar halo with respect to the Galactic disk serves as independent evidence for the accreted origin of the halo, as it is difficult to achieve such a configuration through purely in situ processes. The two density breaks corroborate the scenario in which GSE goes through two apocentric passages before merging with the galaxy, as predicted by Naidu et al. (2021b). The second apocenter occurs at smaller galactocentric radius than the first due to dynamical friction, and the distance between them is an observable quantity, which we have measured in this work, that can be used to constrain properties of the merger. For example, the total mass ratio and the concentration of the progenitor should strongly affect how many apocentric passages the progenitor can go through before fully merging, and how deeply the progenitor plunges at each passage. Future work will address such questions.

Second, the association of the tilted stellar halo with an ancient (t_lb ∼ 10 Gyr) merger holds a clue to the shape of the dark matter halo. Han et al. (2022) show that a present-day spatial asymmetry of GSE, coupled with its estimated old age, implies that at least some fraction of the underlying dark matter halo is also tilted (hence necessarily aspherical) and oriented in a similar direction as GSE. Understanding the dark matter distribution in the galaxy is important for both Galactic archeology and local direct detection experiments, as the nature of the dark matter halo affects the local DM density and velocity distributions (e.g., Evans et al. 2019; Necib et al. 2019). Furthermore, while the evolution of the stellar disk induced by the dark matter halo has been studied before (e.g., Debattista et al. 2013; Yurin & Springel 2015; Dillamore et al. 2022), the stability and growth of the disk specifically in a tilted, triaxial dark matter halo is an open question, and will be investigated in future studies.

Lastly, it has been recognized that the Large Magellanic Cloud (LMC) likely 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). While this effect is most prominent at larger distances (r_Gal > 50 kpc) than what we consider here (see Lucchini et al. 2021 for an alternative scenario), the direction of overdensities due to the LMC and GSE have similar on-sky orientations. Whether this curious alignment is due to coincidence or a manifestation of the large-scale environment that our Galaxy resides in will be addressed by future studies.