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, $\psi$, is generally expressed as

where $G(\tau)$ is an arbitrary function (de Zeeuw, 1985; Dejonghe & de Zeeuw, 1988). $\lambda$ and $\nu$ are the spheroidal coordinates system ($\lambda,\nu$), defined as the solution to the following equation for $\tau$,

where $\alpha$ and $\gamma$ are constants and $(R,z)$ are the Galactocentric cylindrical coordinates, assuming that the position of the Sun is $(R_{\odot},z_{\odot})\equiv(8.3$ kpc, 0) (Gillessen et al., 2009). The contours of $\lambda,\nu=$ constant describe ellipses and hyperbolas, respectively.

For the arbitrary function $G(\tau)$ 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,

where $G_{\rm 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 to make this two-component model in the Stäckel form.

In this work, we adopt $M=10^{12}M_{\odot}$, which is nearly in agreement with several recent measurements, e.g., $0.7\times 10^{12}M_{\odot}$ (Eadie & Jurić, 2019), $1.3\times 10^{12}M_{\odot}$ (Posti & Helmi, 2019; Grand et al., 2019), $1.0\times 10^{12}M_{\odot}$ (Deason et al., 2019), and $1.2\times 10^{12}M_{\odot}$ (Deason et al., 2021). For other parameters, $k$, $\alpha$, $\gamma$, 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 $\simeq 220$ km s^-1, the local mass density at $R=R_{\odot}$ is $0.1\simeq 0.2$ $M_{\odot}$ pc^-3, and the surface mass density at $R=R_{\odot}$ is 70 $M_{\odot}$ pc^-2 within $z\lesssim 1.1$ kpc. We adopt $k=0.09$ for the disk fraction, and $\alpha$ and $\gamma$ 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$. 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 feature that the rotational speed near the Sun, $R\sim 8$ kpc, is about 220 km s^-1 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 in high-$z$ region widely, we consider that the difference of the potentials near the disk does not affect the result critically.

In this work, we define the heliocentric distance as the inverse of parallax. Bailer-Jones et al. (2021) showed the more accurate distance considering the effect of extinction, magnitude and color. We find that the difference between $\varpi^{-1}$ and the distance of Bailer-Jones et al. (2021) lies within 1$\sigma$ for 75% of our sample and within 2$\sigma$ for 90%, under the condition of $\sigma_{\varpi}/\varpi<$ 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_{2},I_{3})$. $E$ and $I_{2}$ are defined as

where $v_{\phi}$ 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_{3}$ is expressed as

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_{2},I_{3})$ is shown in Figure 3 and 4. The rectangle in Figure 3 means that stars inside it belong to the GSE-like subsample.

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

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

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 $i$-th star at position ${\bm{x}}$ is expressed as (de Zeeuw, 1985; Dejonghe & de Zeeuw, 1988; Sommer-Larsen & Zhen, 1990; Chiba & Beers, 2001)

We reconstruct the density distribution of the sample halo stars, $\rho({\bm{x}})$, by overlaying the orbital densities with proper weights, $c_{i}$.

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 ${\bm{x}}_{j}$ is expressed as

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 factor $c_{i}$ is 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, $\rho_{orb,i}({\bm{x}}_{i})$.

This procedure means that a higher weight is given to a star located at the position of lower $\rho_{orb,i}({\bm{x}}_{i})$ in the course of its orbital motion: since the probability to find such a star is low at its current position, 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 $\rho_{orb,i}({\bm{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 $\langle v_{\phi}\rangle$ as follows.

The sign in Eq. (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 Eq.(10) is shown in Figure 7, where the total amplitude of the density is normalized by the value at $(R,z)=(10~{}{\rm 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:

where $R_{c}(=10$ kpc$)$ denotes the radius of the density normalization.

Figure 8 shows the reproduced density distribution along the Galactic plane (solid lines) and fitting curves (dashed lines). We find that all the distributions have a discontinuity at $R\sim 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\sim 8$ kpc, i.e., a star with apocentric radius $r_{\rm apo}\lesssim 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\times 10^{5}~{}{\rm km^{2}~{}s^{-2}}$, 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.

We first determine the value of the slope $\alpha$ by the root-mean-squares fitting of the radial distributions in Figure 8 into the power-law function, $\rho\propto R^{-\alpha}$, i.e. the density distribution at $z=0$ in Eq. (15). We obtain $\alpha=3.2-3.4$ in the entire metallicity ranges. Specifically, the value of $\alpha$ in each metallicity range is $3.42^{+0.17}_{-0.07}$, $3.23^{+0.11}_{-0.08}$, $3.34^{+0.11}_{-0.11}$, and 3$.42^{+0.09}_{-0.13}$, from 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 $\theta=0^{\circ}-90^{\circ}$ at $r=8,15,25$ kpc (solid lines). The dashed lines correspond to the ellipsoidal fitting curves to these density distributions with a parameter, $q$: while these lines be constant along $\theta$ in the case of a spherically symmetric case of $q=1$, a smaller $q$, i.e., more flattened axial ratio, leads to a more negative density gradient with increasing $\theta$. We find that the stellar halo in a more metal-poor range has a rounder shape, and metal-rich range is more flattened. This result is in agreement with previous works before Gaia using much smaller numbers of sample stars (Chiba & Beers, 2000; Carollo et al., 2007).

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

It is worth remarking from Figure 7 and 9 that the actual shapes of equidensity contours deviate from ellipsoidal profiles in Eq. (15): they appear somewhat boxy or peanut-like shapes, so that the derived density distributions depicted in Figure 9 have a peak around $\theta\sim 15^{\circ}$, whereas simple ellipsoidal profiles should decrease monotonically with $\theta$. We note that these properties of the density shapes are related to the lack of halo stars with small $I_{3}$ or small orbital angle $\theta_{\rm max}$ shown in Figure 4 and 6 (see Subsection 4.1.2). We quantify the degree of the derivation from ellipsoidal shapes by $\Delta$ defined as

where $\rho_{\rm obs}$ is the density calculated from Eq. (10) (solid lines in Figure 9), $\rho_{\rm fit}$ is the fitted ellipsoid given in Eq. (15) (dashed lines), and $N_{\theta}=31$ is the number of bins for $\theta_{i}$. Figure 11 shows the distributions of $\Delta$ having a peak sharply at $r=10$ kpc except for ${\rm[Fe/H]}<-2.2$; the latter range shows spherical shapes with little deviation from the ellipsoidal distribution over wide range of $r$. At $r>10$ kpc, $\Delta$ is decreasing with $r$, where $\Delta$ is largest in $-1.8<{\rm[Fe/H]}<-1.4$.

These behaviors in Figure 11 suggest that some substructures that perturb the ellipsoidal profile exist at around $r=10$ kpc in ${\rm[Fe/H]}>-2.2$ and $r>10$ kpc in $-1.8<{\rm[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 Monoceros Stream (Newberg et al., 2002; Yanny et al., 2003; Ibata et al., 2003). Alternatively, a large $\Delta$ may suggest that the shape of the stellar halo is intrinsically boxy or peanut-like shaped rather than spheroid. We discuss the origin of this characteristic shape of the halo in 4.1.2 in more detail.

Next, using Eq. (13), we obtain the global rotational velocity distribution, $\left<v_{\phi}\right>$. The results for relatively metal-rich ranges are shown in Figure 12. We find that stars in $-1.2<{\rm[Fe/H]}<-1.0$ have the relatively large rotation, whose maximum is over $100~{}{\rm km~{}s^{-1}}$. This disk-like structure is local, $R\lesssim 10$ kpc and $z\lesssim 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.

Figure 13 shows the rotational velocity distribution along the Galactic plane $z=0$. The mean rotational velocity in $-1.2<{\rm[Fe/H]}<-1.0$ is more than twice as ones of other metallicity ranges at $8<R<15$ kpc. On the other hand, over $R>20$ kpc, $\left<v_{\phi}\right>$ is almost zero in any metallicity ranges.

Note that Fig. 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 $\sim 150~{}{\rm km\ s^{-1}}$, while Fig. 13 shows that the maximum velocity is $\sim 110~{}{\rm km\ s^{-1}}$. This difference is due to the contamination of halo stars, which move relatively randomly with less rotation. Thus Fig. 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, Fig. 12 and 13 show that for the metal-poor stars, especially in ${\rm[Fe/H]}<-2.2$, the velocity field in the outer region, $R>15-20$ kpc, is slightly but systematically negative below zero about $1\sigma$ 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 sample of previous researches shows a significant retrograde rotation at the solar position, the rotational value at $R\sim 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 relatively metal-rich stellar halos in the range of ${\rm[Fe/H]}>-1.8$ have smaller axial ratios $q$, i.e. more flattened, than more metal-poor halos with ${\rm[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 the metallicity ranges below ${\rm[Fe/H]}=-1$.

The flattened structures in relatively metal-rich ranges may, at least partly, be affected by the presence of the MWTD at ${\rm[Fe/H]}<-1$. The disk-like motion of the MWTD is appreciable in the $\left<v_{\phi}\right>$ distribution of $-1.2<{\rm[Fe/H]}<-1.0$ range (Figure 12 and 13), with the radial and vertical sizes of $\sim 14$ kpc and vertical of $\sim 4$ kpc, respectively. If we adopt a condition, $\left<v_{\phi}\right>/\sigma_{\phi}>1$, with $\sigma_{\phi}=50~{}{\rm km~{}s^{-1}}$ as a typical velocity dispersion of the thick disk in $\phi$ direction (Chiba & Beers, 2000), these figures also show that such a disk-like motion becomes weak quickly in $-1.4<{\rm[Fe/H]}<-1.2$, where $\left<v_{\phi}\right>/\sigma_{\phi}$ is smaller than 1 at all radii, so the effect of the MWTD may be sufficiently weak at ${\rm[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 this MWTD-like feature permeates 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 parts may be explained within the framework of the CDM model (White & Rees, 1978; Bekki & Chiba, 2000): a flattened structure of a relatively metal-rich halo was formed by dissipative merging of few massive clumps at early stage of the MW’s evolution. Alternatively, this flattened structure can be a dynamically heated, pre-existing 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, more metal-poor halos hold a structure closer to a sphere, and this characteristics 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 dSphs, 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 accretion of stars. It 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. It 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 ${\rm[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, $\alpha$, the derived values in this work are almost the same, $3.2-3.4$, in 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 $\alpha<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 that 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 $\alpha$ may not be strongly affected by later accretion process of small clumps and this work obtain the almost same values of $\alpha$ for different metallicity ranges.

Here we discuss the origin of the non-ellipsoidal profile peaked at around $r\sim 10$ kpc, as shown in Figure 7. We note that the presence of the MWTD is insufficient to explain this feature, because the non-ellipsoidal structure is most prominent within $-1.8<{\rm[Fe/H]}<-1.4$, whereas the MWTD dominates in $-1.4<{\rm[Fe/H]}<-1.0$ (Figure 12 and 13). We consider two possibilities to explain 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 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 regions correspond to those, where the current deviation measure from an ellipsoidal shape, $\Delta$, shows a peak and is relatively small, respectively.

Figure 9 shows that the density distributions along $r=8$ and 15 kpc is curved upward at around $\theta\sim 15^{\circ}$ and downward over $\theta>60^{\circ}$. In fact, the density deficiency at large $\theta$ 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 $\theta$ are represented in the solar neighborhood and this effect may have been non-negligible in previous works using only a small number of sample stars available before Gaia. However, our reconstruction using much numerous stars cataloged in Gaia leads the 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_{\rm max}<10^{\circ}$. Carollo & Chiba (2021) showed that the lack of these low-$\theta$ (or low-$I_{3}$) is not artificial: if stars with small orbital inclination 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-$\theta$ stars possess highly eccentric, radial orbits, they can surely be identified, if they exist, in the current survey volume. Then, the lack of such stars near the Galactic plane leads to the reconstructed global density distribution having a boxy shape.

Second, the non-ellipsoidal structure of the halo may reflect the Monoceros Stream (Newberg et al., 2002; Yanny et al., 2003; Ibata et al., 2003). This stream feature is a ring-shaped overdensity with heliocentric radius of 6 kpc in the south and 9 kpc in the north (Morganson et al., 2016) and 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 the similar vertical range and the non-ellipsoidal 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 Monoceros Stream is centered at ${\rm[Fe/H]}=-0.95$ and spreads $-1.1<{\rm[Fe/H]}<-0.7$. This metallicity is quite different from the range where the non-ellipsoidal structure dominates in Figure 11, $-1.8<{\rm[Fe/H]}<-1.4$. This suggests that the non-ellipsoidal structure of the stellar halo obtained here is not due to the Monoceros Stream.

We thus conclude that the non-ellipsoidal structure of the stellar halo, which mainly reveals 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 same masses (dry and major merger) forms a boxy shaped elliptical galaxy (Naab et al., 2006; Bell et al., 2006; Cox et al., 2006). Even including gas in merging process, while gas dissipates into the equatorial plane and form more metal-rich stars with disk-like kinematics, pre-existing stars in progenitor satellites can end up with a boxy structure after merging, which consists of 3D box orbits. Also, dynamical heating of this formed disk through later merging of satellites is another source of halo stars, where heating process can be more effective in the outer disk regions due to weaker vertical gravitational force, so that the heated debris would look like a peanut-shaped structure, like 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 non-ellipsoidal 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\sim 8$ kpc (cf. Fig. 8 and 15). These discontinuities are caused by the observational bias that stars orbiting only in the inner region of $R\lesssim 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|<(R-R_{\odot})\tan 165^{\circ}$ and “the outer region” is $|z|<(R-R_{\odot})\tan 15^{\circ}$, which respectively are the X-shape vacant regions in Fig. 1. Stars staying either of these regions are never included the observational sample.

Fig. 18 shows the $E$-$I_{3}$ space distribution, similar to Fig. 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 locate, depending on their $|L_{z}|$. For example, stars staying in the inner region with $|L_{z}|>1.0\times 10^{3}{\rm\ kpc\ km\ s^{-1}}$ locate in the shaded region below the green line. In other words, stars shown with green and blue points (observed stars with $|L_{z}|>1.0\times 10^{3}{\rm\ kpc\ km\ s^{-1}}$) cannot locate in the shaded region below the green line due to the bias. It is clear that the shaded region in Fig. 18 demonstrates the observational bias in the inner region.

Fig. 19 also shows the $E$-$I_{3}$ space distribution, but the shaded region indicates the area where stars staying only in the outer region. In this case, for example, red-pointed stars (with $|L_{z}|>2.0\times 10^{3}{\rm\ kpc\ km\ s^{-1}}$) cannot locate in the shaded region below the red line; this prohibited region is so small. For another example, stars with $|L_{z}|>7.0\times 10^{3}{\rm\ kpc\ km\ s^{-1}}$ are greatly constrained such that they cannot locate in the wide shaded region below the blue line, but such stars with high-$|L_{z}|$ do not exist at least within $R\sim 35$ kpc, provided a circular rotation velocity there remains about $200~{}{\rm km\ s^{-1}}$ (cf. Fig. 2 and 3)). Additionally, the black points ($|L_{z}|<2.0\times 10^{3}{\rm\ kpc\ km\ s^{-1}}$, 96% of the sample) are not limited in the $E$-$I_{3}$ plane, which means that they cannot stay in the outer region. It suggests that stars staying only in the outer region are few, so the corresponding discontinuity is not expected to appear.

Carollo & Chiba (2021) supports the current argument. It showed that the $\sqrt{2I_{3}}$ distribution quickly damps at around $\sqrt{2I_{3}}\sim 0.5\times 10^{3}{\rm\ kpc\ km\ s^{-1}}$. Stars staying in the outer region need to have small inclination $\theta_{0}$. Fig. 5 shows that $\theta_{0}$ strongly correlates with $I_{3}$. Thus the damp of $\sqrt{2I_{3}}$ distribution suggests that stars staying in the outer region is very few. On the other hand, stars staying in the inner region need to have low $E$. In Fig. 3, the dotted line descend beyond the frame around $|L_{z}|\sim 0-1.0\times 10^{3}{\rm\ kpc\ km\ s^{-1}}$, but the data points do not follow this lower-limit lines. The gap between the dotted line and the fiducial points suggests that there are many inner stars not captured by the observation. Therefore, the discontinuities at around $R\sim 8$ kpc appear while there are no discontinuities in the outer region.