The Global Structure of the Milky Way's Stellar Halo Based on the Orbits of Local Metal-poor Stars

In this work, for the purpose of calculating each orbital density analytically, as described below, we adopt the Stäckel form for the MW's gravitational potential. The Stäckel potential, ψ, is generally expressed as

$$\begin{eqnarray}&&\psi (\lambda ,\nu )=-\displaystyle \frac{(\lambda +\gamma )G(\lambda )-(\nu +\gamma )G(\nu )}{\lambda -\nu },\end{eqnarray} \tag{ 2 }$$

where G(τ) is an arbitrary function (de Zeeuw 1985; Dejonghe & de Zeeuw 1988). λ and ν are the spheroidal coordinate system (λ, ν), defined as the solution to the following equation for τ,

$$\begin{eqnarray}&&\displaystyle \frac{{R}^{2}}{\tau +\alpha }+\displaystyle \frac{{z}^{2}}{\tau +\gamma }=1,\end{eqnarray} \tag{ 3 }$$

where α and γ are constants and (R, z) are the Galactocentric cylindrical coordinates, assuming that the position of the Sun is (R_⊙, z_⊙) ≡ (8.3 kpc, 0) (Gillessen et al. 2009). The contours of λ, ν = constant describe ellipses and hyperbolas, respectively.

For the arbitrary function G(τ) in Equation (2), we consider the contribution of gravitational force from the stellar disk and dark halo components. Here we adopt a Kuzmin–Kutuzov potential (Dejonghe & de Zeeuw 1988; Batsleer & Dejonghe 1994; Chiba & Beers 2001) for both components given in the first and second terms below, respectively,

$$\begin{eqnarray}&&G(\tau )=\displaystyle \frac{{G}_{\mathrm{grav}}{kM}}{\sqrt{\tau }+\sqrt{-\gamma }}+\displaystyle \frac{{G}_{\mathrm{grav}}(1-k)M}{\sqrt{\tau +b}+\sqrt{-\gamma +b}},\end{eqnarray} \tag{ 4 }$$

where G_grav is the gravitational constant, M is the total mass of the MW, k is the ratio of the disk mass to the total mass, and b is the parameter for making this two-component model in the Stäckel form.

In this work, we adopt M = 10¹² M_⊙, which is nearly in agreement with several recent measurements, e.g., 0.7 ×10¹² M_⊙ (Eadie & Jurić 2019), 1.3 × 10¹² M_⊙ (Grand et al. 2019; Posti & Helmi 2019), 1.0 × 10¹² M_⊙ (Deason et al. 2019), and 1.2 × 10¹² M_⊙ (Deason et al. 2021). For the other parameters, k, α, γ, and b, we take the same values as those in Chiba & Beers (2001) so as to make the current model resemble the real galaxy: the circular velocity at R > 4 kpc is constant and ≃ 220 km s⁻¹, the local mass density at R = R_⊙ is 0.1 ≃ 0.2 M_⊙ pc⁻³, and the surface mass density at R = R_⊙ is 70 M_⊙ pc⁻² within z ≲ 1.1 kpc. We adopt k = 0.09 for the disk fraction, and α and γ are associated with the disk's scalelength in R, a_D , and scaleheight in z, c_D , respectively, as $\alpha =-{a}_{D}^{2}$ and $\gamma =-{c}_{D}^{2}$. The corresponding scales for the dark halo, a_H and c_H , are set as ${a}_{H}^{2}={a}_{D}^{2}+b$ and ${c}_{H}^{2}={c}_{D}^{2}+b$, respectively. These values of scales are given as c_D = 0.052 kpc, c_H = 17.5 kpc, c_D /a_D = 0.02, and c_H /a_H = 0.99.

The rotation curve of this gravitational potential is shown in Figure 2. It captures the features that the rotational speed near the Sun, R ∼ 8 kpc, is about 220 km s⁻¹, and the curve farther than 8 kpc is almost flat. They are consistent with the observational properties (e.g., Mróz et al. 2019; Sofue 2020), so we believe this model is appropriate for the following analysis. It is true that the current Stäckel potential does not perfectly represent the actual Galactic potential, especially near the disk plane, in contrast to other models, such as MWPotential2014 (Bovy 2015). However, since we focus on the distribution of the stellar halo, which is mostly spread out widely in the high-z region, we consider that the difference of the potentials near the disk does not affect the result critically.

Zoom In Zoom Out Reset image size

In this work, we define the heliocentric distance as the inverse of parallax. Bailer-Jones et al. (2021) showed a more accurate distance, considering the effects of extinction, magnitude, and color. We find that the difference between ϖ⁻¹ and the distance of Bailer-Jones et al. (2021) lies within 1σ for 75% of our sample and within 2σ for 90%, under the condition of σ_ϖ /ϖ < 0.2. Thus, we regard the inverse of parallax as the distance in this work.

The Stäckel potential has a unique feature that all three integrals of motion can be written analytically (de Zeeuw & Lynden-Bell 1985). This axisymmetric system allows three integrals of motion, (E, I₂, I₃). E and I₂ are defined as

$$\begin{eqnarray}&&E=\displaystyle \frac{1}{2}| {\boldsymbol{v}}{| }^{2}+\psi ({\boldsymbol{x}})\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{I}_{2}=\displaystyle \frac{1}{2}{L}_{z}^{2}\,\mathrm{with}\,{L}_{z}={{Rv}}_{\phi },\end{eqnarray} \tag{ 6 }$$

where v_ϕ is the azimuthal component of the velocity vector. The remaining third integral cannot generally be described analytically in any axisymmetric system, except for the case of a Stäckel potential. In the case of a Stäckel model, the third integral I₃ is expressed as

$$\begin{eqnarray}&&{I}_{3}=\displaystyle \frac{1}{2}({L}_{x}^{2}+{L}_{y}^{2})+(\gamma -\alpha )\left[\displaystyle \frac{1}{2}{v}_{z}^{2}-{z}^{2}\displaystyle \frac{G(\lambda )-G(\nu )}{\lambda -\nu }\right],\end{eqnarray} \tag{ 7 }$$

where (x, y, z) are the Cartesian coordinates with $x=R\cos \phi$ and $y=R\sin \phi$.

The distribution of our stellar data in phase space (E, I₂, I₃) is shown in Figures 3 and 4. The rectangle in Figure 3 means that the stars inside it belong to the GSE-like subsample.

Zoom In Zoom Out Reset image size

The physical meaning of I₃ is the generalized horizontal angular momentum L_∥: in a spherically symmetric case, two constants α and γ are equal, so ${I}_{3}=(1/2)({L}_{x}^{2}+{L}_{y}^{2})=(1/2){L}_{\parallel }^{2}\,=(1/2)({L}^{2}-{L}_{z}^{2})$. The total angular momentum L is an integral of motion in a spherical model, so L_∥ is also an integral of motion. Therefore, $\sqrt{2{I}_{3}}$ can be regarded as the generality of L_∥ (this is why the horizontal axis in Figure 4 describes $\sqrt{2{I}_{3}}$, not simply I₃).

Zoom In Zoom Out Reset image size

It it known that I₃ has a correlation with the maximum inclination angle, ${\theta }_{\max }$, of an orbit, where ${\theta }_{\max }$ is defined as the angle from the Galactic plane. Our ${\theta }_{\max }$ versus I₃ diagram is shown in Figure 5. We find that the correlation appears clearly. Additionally, Figure 6 shows the number distribution of ${\theta }_{\max }$ for the current sample stars. It suggests that the distribution is peaked at ${\theta }_{\max }\sim 15^\circ$ and smoothly declines at larger ${\theta }_{\max }$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since the integrals of motion contain complete information about the orbit of a star, the fact that all integrals of motion are described analytically means that we can reproduce the stellar orbit accurately and completely. Under the Stäckel potential, the orbital density of the i-th star at position x is expressed as (de Zeeuw 1985; Dejonghe & de Zeeuw 1988; Sommer-Larsen & Zhen 1990; Chiba & Beers 2001)

$$\begin{eqnarray}&&{\rho }_{\mathrm{orb},i}({\boldsymbol{x}})=\displaystyle \frac{2\sqrt{2}}{R}\displaystyle \frac{1}{\sqrt{{N}_{i}(\lambda )(\lambda +\gamma )}\sqrt{-{N}_{i}(\nu )(\nu +\gamma )}\sqrt{{I}_{2,i}}}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\mathrm{where}\,{N}_{i}(\tau )=G(\tau )-\displaystyle \frac{{I}_{2,i}}{\tau +\alpha }-\displaystyle \frac{{I}_{3,i}}{\tau +\gamma }+{E}_{i}.\end{eqnarray} \tag{ 9 }$$

We reconstruct the density distribution of the sample halo stars, ρ( x ), by overlaying the orbital densities with proper weights, c_i .

$$\begin{eqnarray}&&\rho ({\boldsymbol{x}})=\sum _{i=1}^{N}{c}_{i}{\rho }_{\mathrm{orb},i}({\boldsymbol{x}}).\end{eqnarray} \tag{ 10 }$$

Here, the orbit-weighting factors c_i are determined by a maximum likelihood method, as follows. When the i-th star has integrals of motion (E_i , I_2,i , I_3,i ), the probability P_ij that the star is observed at the position x _j is expressed as

$$\begin{eqnarray}&&{P}_{{ij}}=\displaystyle \frac{{c}_{i}{\rho }_{\mathrm{orb},i}({{\boldsymbol{x}}}_{j})}{\sum _{k=1}^{N}{c}_{k}{\rho }_{\mathrm{orb},k}({{\boldsymbol{x}}}_{j})}.\end{eqnarray} \tag{ 11 }$$

The i-th integrals of motion are calculated from the i-th observed star, so the case of i = j realizes for any i. Therefore, the probability P_ij should be maximized at i = j for arbitrary j. We define the weights c_i so that ${\prod }_{i=1}^{N}{P}_{{ii}}$ is maximized (Sommer-Larsen & Zhen 1990; Chiba & Beers 2000, 2001).

The orbit-weighting factors c_i are designed to be proportional to the inverse of the time that the i-th star spends at the currently observed position, or the inverse of the i-th orbital density at the position where the i-th star is observed, ρ_orb,i ( x _i ).

$$\begin{eqnarray}&&{c}_{i}\propto \displaystyle \frac{1}{{\rho }_{\mathrm{orb},i}({{\boldsymbol{x}}}_{i})}\end{eqnarray} \tag{ 12 }$$

This procedure means that a higher weight is given to a star located at the position of lower ρ_orb,i ( x _i ) in the course of its orbital motion: since the probability of finding such a star at its current position is low, the fact that it is actually observed there leads to the assignment of a higher orbit-weighting factor c_i . Conversely, a lower c_i is assigned to a star currently located at higher ρ_orb,i ( x _i ), e.g., near the edge of its orbit, where a star stays for a longer time.

In addition to the global density structure of metal-poor halo stars, we also construct their global rotational velocity distribution. We calculate the values of mean rotational velocity 〈v_ϕ 〉 as follows:

$$\begin{eqnarray}&&\langle {v}_{\phi }({\boldsymbol{x}})\rangle =\displaystyle \frac{1}{\rho ({\boldsymbol{x}})}\sum _{i}{c}_{i}{\rho }_{\mathrm{orb},i}({\boldsymbol{x}}){v}_{\phi ,i}({\boldsymbol{x}})\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\mathrm{where}\,{v}_{\phi ,i}({\boldsymbol{x}})=\pm \sqrt{\displaystyle \frac{2{I}_{2,i}}{{R}^{2}}}.\end{eqnarray} \tag{ 14 }$$

The sign in Equation (14) is determined to match the rotational direction of the observational i-th star.

The reproduced global density distribution of the halo sample expressed by Equation (10) is shown in Figure 7, where the total amplitude of the density is normalized by the value at (R, z) = (10 kpc, 0). The density contours hold nearly ellipsoidal shapes in all metallicity ranges, although the detailed contour pattern looks like a somewhat boxy or peanut-like shape. In order to characterize these global density distributions quantitatively, we fit each of them to a simple ellipsoidal profile:

$$\begin{eqnarray}&&\rho ({\boldsymbol{x}})={\left[\displaystyle \frac{{R}_{c}^{2}}{{R}^{2}+{\left(z/q\right)}^{2}}\right]}^{\alpha /2},\end{eqnarray} \tag{ 15 }$$

where R_c (= 10 kpc) denotes the radius of the density normalization.

Zoom In Zoom Out Reset image size

Figure 8 shows the reproduced density distribution along the Galactic plane (solid lines) and fitting curves (dashed lines). We find that all of the distributions have a discontinuity at R ∼ 8 kpc, and the density below this radius falls short of the fitting curves in all stellar metallicity ranges. The cause of this continuity is due to an observational bias (Sommer-Larsen & Zhen 1990): a star moving only within the region inside the position of the Sun, r ∼ 8 kpc, i.e., a star with apocentric radius r_apo ≲ 8 kpc, cannot reach the observable region near the Sun. Therefore, the observational data are devoid of these stars. This bias is also seen in Figure 3, where stars with E < − 1.5 ×10⁵ km² s⁻², i.e., those orbiting inside the solar position, are lacking in the current local sample (see also Myeong et al. 2018b; Carollo & Chiba 2021). We thus avoid the innermost distribution r < 8 kpc for the following fittings.

Zoom In Zoom Out Reset image size

We first determine the value of the slope α by the root-mean-square fitting of the radial distributions in Figure 8 to the power-law function ρ ∝ R^−α , i.e., the density distribution at z = 0 in Equation (15). We obtain α = 3.2 − 3.4 in the entire metallicity range. Specifically, the values of α in each metallicity range are ${3.42}_{-0.07}^{+0.17}$, ${3.23}_{-0.08}^{+0.11}$, ${3.34}_{-0.11}^{+0.11}$, and ${3.42}_{-0.13}^{+0.09}$, from the metal-rich to metal-poor ranges.

Next, we determine the value of an axial ratio, q, at each position r and metallicity range. Figure 9 shows the reproduced density distribution along θ = 0° − 90° at r = 8, 15, and 25 kpc (solid lines). The dashed lines correspond to the ellipsoidal fitting curves to these density distributions with a parameter q: while these lines are constant along θ in the case of a spherically symmetric case of q = 1, a smaller q, i.e., a more flattened axial ratio, leads to a more negative density gradient with increasing θ. We find that the stellar halo in a more metal-poor range has a rounder shape, and is more flattened in a metal-rich range. This result is in agreement with previous works before Gaia using much smaller numbers of sample stars (Chiba & Beers 2000; Carollo et al. 2007).

Zoom In Zoom Out Reset image size

Figure 10 shows the radial distribution of q derived in this manner over r = 8–30 kpc. We find that more metal-rich halo stars tend to have a smaller axial ratio q over the entire radii (except for a somewhat reversed trend in the ranges of −1.8 < [Fe/H] < − 1.0 at r > 20 kpc), and that the axial ratio of the most metal-poor halo ([Fe/H] <−2.2) is nearly 1 over the entire radii, i.e., having a nearly spherical shape.

Zoom In Zoom Out Reset image size

It is worth remarking from Figures 7 and 9 that the actual shapes of the equidensity contours deviate from the ellipsoidal profiles in Equation (15): they have somewhat boxy or peanut-like shapes, so that the derived density distributions depicted in Figure 9 have a peak around θ ∼ 15°, whereas simple ellipsoidal profiles should decrease monotonically with θ. We note that these properties of the density shapes are related to the lack of halo stars with a small I₃ or small orbital angle ${\theta }_{\max }$ shown in Figures 4 and 6 (see Section 4.1.2). We quantify the degree of the derivation from ellipsoidal shapes by Δ, defined as

$$\begin{eqnarray}&&{\rm{\Delta }}(r)=\sqrt{\displaystyle \frac{1}{{N}_{\theta }}\sum _{{\theta }_{i}}\displaystyle \frac{{[{\mathrm{log}}_{10}{\rho }_{\mathrm{obs}}(r,{\theta }_{i})-{\mathrm{log}}_{10}{\rho }_{\mathrm{fit}}(r,{\theta }_{i})]}^{2}}{{\sigma }_{{\mathrm{log}}_{10}[{\rho }_{\mathrm{obs}}(r,{\theta }_{i})]}^{2}}},\end{eqnarray} \tag{ 16 }$$

where ρ_obs is the density calculated from Equation (10; the solid lines in Figure 9), ρ_fit is the fitted ellipsoid given in Equation (15; the dashed lines), and N_θ = 31 is the number of bins for θ_i . Figure 11 shows the distributions of Δ having a sharp peak at r = 10 kpc, except for [Fe/H] <−2.2; the latter range shows spherical shapes with little deviation from the ellipsoidal distribution over a wide range of r. At r > 10 kpc, Δ is decreasing with r, where Δ is largest in −1.8 < [Fe/H] < − 1.4.

Zoom In Zoom Out Reset image size

These behaviors in Figure 11 suggest that some substructures that perturb the ellipsoidal profile exist at around r = 10 kpc in [Fe/H] >−2.2 and r > 10 kpc in −1.8 < [Fe/H] < − 1.4. This may include the presence of the metal-weak thick disk (MWTD; e.g., Beers et al. 2002; Carollo et al. 2019) at r < 10 kpc and/or the presence of the Monoceros Stream (Newberg et al. 2002; Ibata et al. 2003; Yanny et al. 2003). Alternatively, a large Δ may suggest that the shape of the stellar halo is intrinsically boxy or peanut-like, rather than spheroid. We discuss the origin of this characteristic shape of the halo in more detail in 4.1.2.

Next, using Equation (13), we obtain the global rotational velocity distribution, 〈v_ϕ 〉. The results for the relatively metal-rich ranges are shown in Figure 12. We find that the stars in −1.2 < [Fe/H] < − 1.0 have a relatively large rotation, whose maximum is over 100 km s⁻¹. This disk-like structure is local, R ≲ 10 kpc and z ≲ 3 kpc, and decays with lowering stellar metallicity. It is consistent with the properties of the MWTD. Outside the disk-like region, the velocity field is almost 0, with no significant prograde and retrograde rotation.

Zoom In Zoom Out Reset image size

Figure 13 shows the rotational velocity distribution along the Galactic plane z = 0. The mean rotational velocity in −1.2 < [Fe/H] < − 1.0 is more than twice those of the other metallicity ranges at 8 < R < 15 kpc. On the other hand, over R > 20 kpc, 〈v_ϕ 〉 is almost zero in all metallicity ranges.

Zoom In Zoom Out Reset image size

Note that Figure 13 shows the mixture distribution of the stellar disk and the stellar halo. Carollo et al. (2019) showed that the rotational velocity of the MWTD reaches ∼ 150 km s⁻¹, while Figure 13 shows that the maximum velocity is ∼ 110 km s⁻¹. This difference is due to the contamination of the halo stars, which move relatively randomly with less rotation. Thus, Figure 13 underestimates the velocity distribution of the pure MWTD. For the same reason, we do not claim that the metallicity of the MWTD is only over [Fe/H] > −1.2.

Additionally, Figures 12 and 13 show that for the metal-poor stars, especially in [Fe/H] <−2.2, the velocity field in the outer region, R > 15 − 20 kpc, is slightly but systematically negative below zero about 1σ uncertainty, but this feature is not thought to be real for the MW's stellar disk. This rotational property for very metal-poor stars in the outer region is in agreement with the retrograde rotation of the outer-halo sample stars (e.g., Carollo et al. 2010). While the samples of previous research show a significant retrograde rotation at the solar position, the rotational value at R ∼ 30 kpc in this work is reduced at the currently derived one under the angular-momentum conservation.

As clearly shown in Figure 10, the global shapes of the relatively metal-rich stellar halos in the range of [Fe/H] >−1.8 have smaller axial ratios q, i.e., are more flattened, than the more metal-poor halos with [Fe/H] <−1.8 over all radii. Also, at radii below 18 kpc, we find that a more metal-rich halo shows a more flattened shape for all of the metallicity ranges below [Fe/H] = − 1.

The flattened structures in the relatively metal-rich ranges may, at least partly, be affected by the presence of the MWTD at [Fe/H] <−1. The disk-like motion of the MWTD is appreciable in the 〈v_ϕ 〉 distribution of the −1.2 < [Fe/H] < − 1.0 range (Figures 12 and 13), with radial and vertical sizes of ∼14 kpc and ∼4 kpc, respectively. If we adopt a condition, 〈v_ϕ 〉/σ_ϕ > 1, with σ_ϕ = 50 km s⁻¹ as a typical velocity dispersion of the thick disk in ϕ direction (Chiba & Beers 2000), these figures also show that such a disk-like motion quickly becomes weak in −1.4 < [Fe/H] < − 1.2, where 〈v_ϕ 〉/σ_ϕ is smaller than 1 at all radii, so the effect of the MWTD may be sufficiently weak at [Fe/H] between −1.4 and −1.2. The fact that the derived relatively metal-rich halo sample shows a flattened shape over the regions beyond where this MWTD-like feature permeates suggests that this flattened structure is a real halo structure. These flattened parts of the halo may correspond to the in situ halo, in which stars have been formed inside a main parent halo.

The formation of this flattened halo part may be explained within the framework of the cold dark matter (CDM) model (White & Rees 1978; Bekki & Chiba 2000): a flattened structure of a relatively metal-rich halo was formed by dissipative merging of a few massive clumps at an early stage of the MW's evolution. Alternatively, this flattened structure can be a dynamically heated, preexisting disk, driven by falling satellites (e.g., Benson et al. 2004; Purcell et al. 2010; McCarthy et al. 2012; Tissera et al. 2013; Cooper et al. 2015).

On the contrary, the more metal-poor halos have a structure closer to a sphere, and this characteristic suggests that it is the ex situ halo component. The presence of this metal-poor, spherical halo is consistent with the MW's formation scenario under the CDM model: the outer halo was formed by the later accretion of small clumps resembling metal-poor dwarf spheriodals, which are stripped by tidal force from the MW's dark halo (Bekki & Chiba 2001; Carollo et al. 2007).

Recent numerical simulations based on CDM models provide similar chemodynamical properties of stellar halos to those reported here. Monachesi et al. (2019) simulated the formation process of MW-like stellar halos with accretions of stars. The simulation showed that the axial ratios of the halos are typically 0.6 to 0.9 at R = 20–40 kpc. This simulated range is consistent with our result shown in Figure 10. The simulation also showed that the shapes of MW-like stellar halos can be prolate in outer regions, where the ex situ component dominates. Our result of [Fe/H] <−2.2 range has the same shape, so it is considered that the ex situ component is dominant in this metallicity range.

Regarding the radial density slope, α, the derived values in this work are almost the same, 3.2–3.4, in the different metallicity ranges, though there may exist varieties in the formation processes. Rodriguez-Gomez et al. (2016) and Monachesi et al. (2019) showed that the accreted halo components have a shallower density slope with α < 3.0, in the radius range 15 < r < 100 kpc, than the in situ components. On the other hand, in the range 8 < r < 30 kpc, with which the current work is concerned, the ex situ halo should be less dominant than in the range of 15 < r < 100 kpc. Thus, the values of α may not be strongly affected by the later accretion process of small clumps, and this work obtains almost the same values of α for the different metallicity ranges.

Here we discuss the origin of the nonellipsoidal profile peaked at around r ∼ 10 kpc, as shown in Figure 7. We note that the presence of the MWTD is insufficient to explain this feature, because the nonellipsoidal structure is most prominent within −1.8 < [Fe/H] < − 1.4, whereas the MWTD dominates in −1.4 < [Fe/H] < − 1.0 (Figures 12 and 13). We consider two possibilities for explaining this feature.

First, the intrinsic structure of the inner part of the MW's stellar halo forms a boxy/peanut-like shape rather than the ellipsoidal shape at r < 20 kpc. Carollo et al. (2007, 2010) showed that the inner halo component is dominant up to r = 10 − 15 kpc, and the outer halo is dominant at r > 20 kpc. These two regions, r = 10 − 15 kpc and r > 20 kpc, correspond to those where the current deviation measure from an ellipsoidal shape, Δ, shows a peak and is relatively small, respectively.

Figure 9 shows that the density distributions along r = 8 and 15 kpc are curved upward at around θ ∼ 15° and downward over θ > 60°. In fact, the density deficiency at large θ also appeared in the results of Sommer-Larsen & Zhen (1990), Chiba & Beers (2000), and Carollo et al. (2007). It was interpreted due to the observational bias: there is a small probability that such stars orbiting at large θ are represented in the solar neighborhood, and this effect may have been nonnegligible in previous works using only the small number of sample stars available before Gaia. However, our reconstruction using the more numerous stars cataloged in Gaia leads to a similar shape of the density distribution of the MW's stellar halo. There are indeed no discontinuities in the observational inclination distribution shown in Figure 6. Thus, we regard the boxy/peanut-shaped structure as a real feature, not biased.

It is also worth remarking in Figure 6 that there is a sharp drop in the number of orbits at ${\theta }_{\max }\lt 10^\circ$. Carollo & Chiba (2021) showed that the lack of these low-θ (or low-I₃) stars is not artificial: if stars with small orbital inclinations really exist in the MW, such stars surely pass by the solar neighborhood and thus can be included in the current local sample. This fact means that since such low-θ stars possess highly eccentric, radial orbits, they can surely be identified, if they exist, in the current survey volume. As a result, the lack of such stars near the Galactic plane leads to the reconstructed global density distribution having a boxy shape.

Second, the nonellipsoidal structure of the halo may reflect the Monoceros Stream (Newberg et al. 2002; Ibata et al. 2003; Yanny et al. 2003). This stream feature is a ring-shaped overdensity with a heliocentric radius of 6 kpc in the south and 9 kpc in the north (Morganson et al. 2016), and a height of 3 < z < 4 kpc (Ivezić et al. 2008). These spatial distributions generally agree with our results that the density contours in Figure 7 are curved upward to a similar vertical range, and the nonellipsoidal structure, as shown in Figure 11, is dominant over the similar Galactocentric distance.

However, Ivezić et al. (2008) and Morganson et al. (2016) showed that the Monoceros Stream is centered at [Fe/H] = − 0.95 and spreads between −1.1 < [Fe/H] < − 0.7. This metallicity is quite different from the range where the nonellipsoidal structure dominates in Figure 11, −1.8 < [Fe/H] < − 1.4. This suggests that the nonellipsoidal structure of the stellar halo obtained here is not due to the Monoceros Stream.

We thus conclude that the nonellipsoidal structure of the stellar halo, which is mainly revealed in the inner parts of the halo, is a physically real feature, and this shape of the stellar halo is considered to reflect the MW's merger history. It is known that a merger between early-type gas-poor galaxies with the same masses (a dry and major merger) forms a boxy-shaped elliptical galaxy (Bell et al. 2006; Cox et al. 2006; Naab et al. 2006). Even including gas in the merging process, although gas dissipates into the equatorial plane and forms more metal-rich stars with disk-like kinematics, the preexisting stars in the progenitor satellites can end up with a boxy structure after merging, consisting of 3D box orbits. Also, the dynamical heating of this formed disk through the later merging of satellites is another source of halo stars, where the heating process can be more effective in the outer disk regions due to weaker vertical gravitational force, so that the heated debris would have a peanut-shaped structure, like the edge-on view of long-lived bars (Athanassoula 2005). Indeed, the simulation works under CDM scenarios show some boxy/peanut-shaped stellar halos (e.g., Bullock & Johnston 2005; Cooper et al. 2010; Monachesi et al. 2019). Further theoretical works will be needed to quantify this nonellipsoidal structure of the simulated halos and clarify what merging histories of progenitor galaxies are associated with this characteristic structure.

We have mentioned that the reproduced density distribution has remarkable discontinuities at R ∼ 8 kpc (see Figures 8 and 15). These discontinuities are caused by the observational bias that stars orbiting only in the inner region of R ≲ 8 kpc never reach the observable region near the Sun, and are not contained in the currently adopted stellar catalog. However, in principle, it is also conceivable that another discontinuity may appear, caused by stars staying far from the Sun, but we do not find such discontinuities. This subsection explains the properties of the observational bias caused by stars staying inside and outside the locally observable region.

In the following, "the inner region" is defined as $| z| \lt (R\,-{R}_{\odot })\tan 165^\circ$ and "the outer region" is $| z| \lt (R-{R}_{\odot })\tan 15^\circ$, which, respectively, are the X-shape vacant regions in Figure 1. Stars staying in either of these regions are never included in the observational sample.

Figure 18 shows the E–I₃ space distribution, similar to Figure 4, but the color of the points indicates the range of ∣L_z ∣. The shaded region indicates the area where stars staying in the inner region are located, depending on their ∣L_z ∣. For example, stars staying in the inner region with ∣L_z ∣ > 1.0 × 10³kpc km s⁻¹ are located in the shaded region below the green line. In other words, the stars shown with green and blue points (the observed stars with ∣L_z ∣ > 1.0 × 10³kpc km s⁻¹) cannot be located in the shaded region below the green line due to the bias. It is clear that the shaded region in Figure 18 demonstrates the observational bias in the inner region.

Zoom In Zoom Out Reset image size

Figure 19 also shows the E–I₃ space distribution, but the shaded region indicates the area where the stars stay only in the outer region. In this case, for example, red-pointed stars (with ∣L_z ∣ > 2.0 × 10³ kpc km s⁻¹) cannot be located in the shaded region below the red line; this prohibited region is so small. For another example, stars with ∣L_z ∣ > 7.0 × 10³ kpc km s⁻¹ are greatly constrained, such that they cannot be located in the wide shaded region below the blue line, but such stars with high ∣L_z ∣ do not exist at least within R ∼ 35 kpc, provided that the circular rotation velocity there remains about 200 km s⁻¹ (see Figures 2 and 3)). Additionally, the black points (∣L_z ∣ < 2.0 × 10³ kpc km s⁻¹; 96% of the sample) are not limited to the E–I₃ plane, which means that they cannot stay in the outer region. This suggests that stars staying only in the outer region are few, so the corresponding discontinuity is not expected to appear.

Zoom In Zoom Out Reset image size

The recent study of Carollo & Chiba (2021) supports the current argument. It showed that the $\sqrt{2{I}_{3}}$ distribution quickly damps at around $\sqrt{2{I}_{3}}\sim 0.5\times {10}^{3}\mathrm{kpc}\ \mathrm{km}\ {{\rm{s}}}^{-1}$. Stars staying in the outer region need to have a small inclination θ₀. Figure 5 shows that θ₀ strongly correlates with I₃. Thus, the damping of the $\sqrt{2{I}_{3}}$ distribution suggests that the stars staying in the outer region are very few. On the other hand, the stars staying in the inner region need to have low E. In Figure 3, the dotted line descends beyond the frame around ∣L_z ∣ ∼ 0 − 1.0 ×10³ kpc km s⁻¹, but the data points do not follow the lower-limit lines. The gap between the dotted line and the fiducial points suggests that there are many inner stars not captured by observations. Therefore, the discontinuities at around R ∼ 8 kpc appear, while there are no discontinuities in the outer region.