# The global structure of the Milky Way’s stellar halo based on the orbits of
local metal-poor stars

Genta Sato  Astronomical Institute, Tohoku University, Aoba-ku, Sendai
980-8578, Japan  Masashi Chiba  Astronomical Institute, Tohoku University,
Aoba-ku, Sendai 980-8578, Japan

###### Abstract

We analyze the global structure of the Milky Way (MW)’s stellar halo including
its dominant subcomponent, Gaia-Sausage-Enceladus (GSE). The method to
reconstruct the global distribution of this old stellar component is to employ
the superposition of the orbits covering over the large MW’s space, where each
of the orbit-weighting factor is assigned following the probability that the
star is located at its currently observed position. The selected local, metal-
poor sample with ${\rm[Fe/H]}<-1$ using Gaia EDR3 and SDSS DR16 shows that the
global shape of the halo is systematically rounder at all radii in more metal-
poor ranges, such that an axial ratio, $q$, is nearly 1 for ${\rm[Fe/H]}<-2.2$
and $\sim 0.7$ for $-1.4<{\rm[Fe/H]}<-1.0$. It is also found that a halo in
relatively metal-rich range of ${\rm[Fe/H]}>-1.8$ actually shows a
boxy/peanut-like shape, suggesting a major merger event. The distribution of
azimuthal velocities shows a disk-like flattened structure at
$-1.4<{\rm[Fe/H]}<-1.0$, which is thought to be the metal-weak thick disk. For
the subsample of stars showing GSE-like kinematics and at ${\rm[Fe/H]}>-1.8$,
its global density distribution is more spherical with $q\sim 0.9$ than the
general halo sample, having an outer ridge at $r\sim 20$ kpc. This spherical
shape is consistent with the feature of accreted halo components and the ridge
suggests that the orbit of GSE’s progenitor has an apocenter of $\sim 20$ kpc.
Implications for the formation of the stellar halo are also presented.

Milky Way — stellar halo — metallicity — formation history

##  1 Introduction

It is generally thought that the stellar halo of the Milky Way (MW) contains
the fossil record on the MW’s formation history imprinted in its kinematical
and chemical properties. Since Eggen et al. (1962) showed a model of the MW’s
formation process from the kinematic analysis of halo stars in the solar
neighborhood, many studies of this fossil Galactic component have been carried
out using a large number of newly available halo sample, based on considerable
observational efforts for their assembly over the past decades. These include,
for example, the Hipparcos Catalog obtained from the first astrometry
satellite (Perryman et al., 1997), and spectroscopic catalogs such as RAVE
(Steinmetz et al., 2006), SEGUE (Yanny et al., 2009), LAMOST (Cui et al.,
2012; Zhao et al., 2012), GALAH (De Silva et al., 2015), APOGEE (Majewski et
al., 2017), H3 (Conroy et al., 2019), and more. Perhaps the most significant
impacts on this field of research have been brought by the second astrometry
satellite, Gaia. Gaia catalog (Gaia Collaboration et al., 2016, 2018, 2021)
provides trigonometric parallaxes and proper motions for billions of Galactic
stars with unprecedented high accuracy. Based on these astrometry data of
stars in Gaia catalog combined with spectroscopic catalogs, the MW’s new
dynamical maps have been drawn (e.g. Helmi et al., 2018; Belokurov et al.,
2018; Myeong et al., 2018a, b; Beane et al., 2019; Hagen et al., 2019; Iorio &
Belokurov, 2019; Anguiano et al., 2020; Cordoni et al., 2021; Koppelman et
al., 2021), and new aspects of the MW’s formation and evolution history have
been revealed (e.g. Antoja et al., 2018; Sestito et al., 2019; Wyse, 2019). In
particular, Gaia catalog has enabled to discover new substructures in the MW’s
stellar halo (e.g. Koppelman et al., 2019; Myeong et al., 2019; Li et al.,
2020), which are remnants of past merging/accretion events associated in the
MW’s formation history (e.g. Fernández-Trincado et al., 2020; Naidu et al.,
2020, 2021).

These observational information of the MW’s stellar halo have been
investigated in comparison with various numerical simulations, which produce a
number of theoretically predicted halo structures in the MW-sized galaxies,
including their global density profiles and dynamical motions over a large
halo space (e.g. Rodriguez-Gomez et al., 2016; Liang et al., 2021; Monachesi
et al., 2019) as well as the merging/accretion history of satellite galaxies
(e.g. Cooper et al., 2010; Mackereth et al., 2019; Kruijssen et al., 2020;
Santistevan et al., 2020).

One of the most significant substructures in the MW’s stellar halo revealed in
Gaia catalog is Gaia-Sausage-Enceladus (GSE) (Helmi et al., 2018; Belokurov et
al., 2018). GSE shows an elongated distribution in a velocity space with
azimuthal velocities of $v_{\phi}\sim 0$ and is thought as a tracer of a
merger event with a SMC-class galaxy during the MW’s evolution about 10 Gyr
ago (Helmi et al., 2018). The further details of GSE’s dynamical and chemical
properties have been determined observationally (e.g. Feuillet et al., 2020;
Monty et al., 2020), and several simulations have also been explored to
constrain the properties of GSE’s progenitor (e.g. Bignone et al., 2019;
Koppelman et al., 2020; Kim et al., 2021). Evans et al. (2020) showed that a
merger event with a GSE-like progenitor is rare for MW-like galaxies. Thus,
GSE is a unique feature of the MW’s stellar halo, thereby suggesting that the
merger should be essential for the MW’s evolution history.

In previous observational studies, a stellar sample selected and used in the
dynamical analysis is generally confined in the local volume of the Galactic
space, so the global structure of the stellar halo is hardly derived. For
example, the region where Gaia catalog obtains the accurate parallax data is
limited near the Sun, within heliocentric distances of $\sim 4$ kpc (see
Subsection 2.1). To solve this issue and derive the spatial distribution of
halo stars over a large Galactic space, we adopt the method based on stellar
orbits (May & Binney, 1986; Sommer-Larsen & Zhen, 1990; Chiba & Beers, 2000,
2001). This method is designed, based on the superposition of stellar orbits,
for reconstructing the stellar distribution over a large region, up to the
apocentric distances where stars can reach. Chiba & Beers (2000) used about
1200 metal-poor stars mostly selected from Hipparcos catalog for this method.
To step further in this work, we adopt Gaia Early Data Release 3rd (EDR3)
(Gaia Collaboration et al., 2021) in combination with the spectroscopic
information from SDSS DR16 SEGUE, which allows us to have a much larger number
of stars with better astrometric accuracy. Thus, it is expected that our
reconstruction of the global halo distribution, over $r\sim 30$ kpc, is more
accurate and statistically more significant than previous studies.

Moreover, we also apply this method to the stars showing GSE-like kinematics
and derive GSE’s global structure in the Galactic halo. In contrast to
previous works (e.g. Feuillet et al., 2020; Monty et al., 2020), which
investigated GSE’s properties based on their local observational information,
we adopt the strategy using the superposition of orbits for GSE-like member
stars to construct the GSE’s global distribution for the first time.

This paper is organized as follows. Section 2 explains the observational data
with the selection conditions we adopt and the method to reconstruct the
distribution of the MW’s stellar halo as well as the GSE-like subsample. In
Section 3, we show the structure of the reproduced global halo distribution
and the results of the fitting to an ellipsoidal profile characterized by an
index of the power-law radial profile, $\alpha$, and an axial ratio, $q$. In
Section 4, we discuss the implication from the derived distribution and
structures of the stellar halo, and infer, based on the current results, about
the MW’s formation and evolution history. Section 5 summarizes our work and
conclusions.

##  2 Method

###  2.1 Data

We adopt Gaia EDR3 Catalog (Gaia Collaboration et al., 2021) for astrometry
data and SDSS DR16 SEGUE catalog (Jurić et al., 2008; Ivezić et al., 2008;
Bond et al., 2010) for spectroscopic data. We extract the common stellar data
from these catalogs satisfying the following five conditions to adopt for our
analysis.

The first condition is that the stellar metallicity ${\rm[Fe/H]}$ is more
metal-poor than $-1$. This condition is to select stars belonging to the
fiducial stellar halo having low metals and possibly old ages.

The second and third conditions are that the effective temperature and the
surface gravity of a star are within the range $4,500<T_{\rm eff}<7,000$ K and
$3.0<\log{g}<5.2$, respectively. These constraints are to extract the main
sequence stars, following the work of Ivezić et al. (2008), which also derived
the range of the observational quantities that fulfil main sequence stars.

The fourth and fifth conditions are for ensuring the accuracy of the data,
i.e., to avoid stars having large uncertainties in parallaxes and
metallicities, respectively. We select the stars having the relative
uncertainty of the parallax $\sigma_{\varpi}/\varpi$ smaller than 0.2 and the
uncertainty of ${\rm[Fe/H]}$ smaller than 0.1 dex.

Next, we carry out cross-match between these two catalogs using the CDS Xmatch
Service 111http://cdsxmatch.u-strasbg.fr, where the cross-matching condition
for the stellar position is given to be smaller than 1 arcsec. Setting only
this positional condition however leaves some cases that a few number of SDSS
stars match with a single Gaia star. In this case, we impose an additional
condition for the photometry of stars, following the relation between the $G$
band used in Gaia and $(g,i)$ bands in SDSS (Busso et al., 2021):

$\begin{split}G-g=&-0.1064-0.4964(g-i)\\\
&-0.09339(g-i)^{2}+0.004444(g-i)^{3}\end{split}$ (1)

We adopt the star having smallest magnitude difference among the position-
matched candidate stars between the two catalogs.

Figure 1: Upper panel: The meridian distribution of fiducial metal-poor halo
stars adopted in this work selected from the cross-match of Gaia EDR3 with
SDSS DR16 SEGUE catalogs. Different colors correspond to different ranges of
stellar metallicities: red for $-1.4<{\rm[Fe/H]}<-1.0$, orange for
$-1.8<{\rm[Fe/H]}<-1.4$, green for $-2.2<{\rm[Fe/H]}<-1.8$, and blue for
${\rm[Fe/H]}<-2.2$. This correspondence is applied in following figures.
Bottom panel: The cumulative number distribution of the heliocentric distance,
$D$, for the adopted stellar sample divided by the metallicity ranges. The
histograms are normalized so that the total number is 1 for each range.

This spectroscopic catalog (SDSS SEGUE) is used to extract line-of-sight
velocities and spectroscopic metallicities of sample stars. Stars containing
these spectroscopic information with good accuracy are limited to only bright
ones, so the range of $G$-band absolute magnitude of our sample is
$M_{G}\lesssim 8$ mag.

We have extracted 25,093 stars as the fiducial members of the stellar halo
222In comparison, Kim & Lépine (2021) recently adopt $\sim$ 500,000 thick-disk
and halo-like stars with $M_{G}<13$ mag to analyze their reduced proper
motions without the limitation of the availability of spectroscopic data.. The
distributions of these stars are shown in Figure 1. In the following, we
divide the sample stars mainly into four subsamples based on their
metallicities, $-1.4<{\rm[Fe/H]}<-1.0$, $-1.8<{\rm[Fe/H]}<-1.4$,
$-2.2<{\rm[Fe/H]}<-1.8$, and ${\rm[Fe/H]}<-2.2$ (shown with different colors
in Figure 1 and following figures). The number of stars in each subsample is
10,214, 8,484, 4,518, and 1,877, respectively.

Additionally, we also consider the likely member stars belonging to the Gaia-
Sausage-Enceladus (GSE). The GSE stars are distributed around the azimutal
velocities, $v_{\phi}$, near 0 (Belokurov et al., 2018; Myeong et al., 2018a),
which corresponds to the vertical angular momentum $L_{z}$ as small as 0. We
note that strongly-bound stars, namely stars with low orbital energy, are
dominated by the so-called in-situ halo component (e.g. Naidu et al., 2020),
and such stars are not GSE’s members. Thus, we extract the GSE-like stars,
which fulfil the conditions about the binding energy $E$ and the vertical
angular momentum $L_{z}$; $E>-150,000~{}{\rm km^{2}~{}s^{-2}}$ and
$-500<L_{z}<500~{}{\rm kpc~{}km~{}s^{-1}}$ (see Subsection 2.2.1 for the
definition of $E$ and $L_{z}$). The number of the extracted stars is 9,402.

These stars belonging to the GSE-like members are also divided into four
subsamples depending on the metallicity ranges of $-1.4<{\rm[Fe/H]}<-1.0$,
$-1.8<{\rm[Fe/H]}<-1.4$, $-2.2<{\rm[Fe/H]}<-1.8$, and ${\rm[Fe/H]}<-2.2$,
where the number of stars are 3,269, 3,607, 1,827, and 699, respectively.
Although the real GSE may dominate in the intermediate metallicity range,
e.g., $-1.66<{\rm[Fe/H]}<-1.33$ (Belokurov et al., 2018), covering our
subsample of $-1.8<{\rm[Fe/H]}<-1.4$, we also consider other metallicity
ranges to analyze the metallicity dependence of the GSE-like structure and to
compare with the general distribution of the stellar halo.

###  2.2 Dynamical Model

We adopt the method based on the superposition of the discrete stellar orbits
to reproduce and analyze the global distribution of the MW’s stellar halo
following the seminal works by May & Binney (1986) and Sommer-Larsen & Zhen
(1990). Although Gaia catalog provides billions of astrometric data, stars
with accurate parallaxes of $\sigma_{\varpi}/\varpi<0.2$ are available only in
the local volume near the Sun, which is small compared with the entire MW’s
halo space (Figure 1). By calculating the stellar orbits and their
superposition with appropriate weighting, we can construct the MW’s halo
distribution beyond the region where we can observe stars locally and
individually. This subsection presents the adopted gravitational potential,
the properties of stellar orbits in it, and then the method to construct the
global distribution of the stellar halo from the orbits.

####  2.2.1 Gravitational potential

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

$\psi(\lambda,\nu)=-\frac{(\lambda+\gamma)G(\lambda)-(\nu+\gamma)G(\nu)}{\lambda-\nu}$
(2)

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$,

$\frac{R^{2}}{\tau+\alpha}+\frac{z^{2}}{\tau+\gamma}=1$ (3)

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.

Figure 2: The rotation curve of the mass model adopted in this work, that is
basically in agreement with the observed one in the MW.

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,

$G(\tau)=\frac{G_{\rm grav}kM}{\sqrt{\tau}+\sqrt{-\gamma}}+\frac{G_{\rm
grav}(1-k)M}{\sqrt{\tau+b}+\sqrt{-\gamma+b}}$ (4)

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$.

Figure 3: The phase-space $(E,L_{z})$ distribution of fiducial halo stars
adopted in this work. The black dotted line corresponds to a circular orbit,
i.e., the lowest orbital energy, $E$, for a given vertical component of
angular momentum, $L_{z}$. The rectangle indicates the selection region of
stars we adopted as GSE-like subsample.

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.

####  2.2.2 Stellar orbits

Figure 4: The distribution of fiducial halo stars in the phase-space
$(\sqrt{2I_{3}},E)$. We note that $\sqrt{2I_{3}}$ has the dimension of angular
momentum and becomes $\sqrt{L_{x}^{2}+L_{y}^{2}}$ in the case of a spherically
symmetric gravitational potential.

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

$\displaystyle E$ $\displaystyle=\frac{1}{2}|{\bm{v}}|^{2}+\psi({\bm{x}})$ (5)
$\displaystyle I_{2}$ $\displaystyle=\frac{1}{2}L_{z}^{2}~{}{\rm
with}~{}L_{z}=Rv_{\phi}$ (6)

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

$I_{3}=\frac{1}{2}(L_{x}^{2}+L_{y}^{2})+(\gamma-\alpha)\left[\frac{1}{2}v_{z}^{2}-z^{2}\frac{G(\lambda)-G(\nu)}{\lambda-\nu}\right]$
(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_{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}$).

Figure 5: The correlation between orbital angles measured from the Galactic
plane, $\theta_{\rm max}$, and $\sqrt{I_{3}}$ in fiducial halo stars.

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}$.

Figure 6: The distribution of orbital angle, $\theta_{\rm max}$. While the
number of stars with $\theta_{\rm max}>15^{\circ}$ shows a monotonous decrease
with increasing $\theta_{\rm max}$, the stars with $\theta_{\rm
max}<15^{\circ}$ are deficient from the extrapolation from this distribution.

Figure 7: The constructed global density distribution calculated from the
entire halo sample for the four different metallicity ranges,
$-1.4<{\rm[Fe/H]}<-1.0$ (upper left), $-1.8<{\rm[Fe/H]}<-1.4$ (upper right),
$-2.2<{\rm[Fe/H]}<-1.8$ (lower left), and ${\rm[Fe/H]}<-2.2$ (lower right
panel). The distribution is normalized at $(R,z)=(10~{}{\rm kpc},0)$.

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)

$\displaystyle\rho_{orb,i}$
$\displaystyle({\bm{x}})=\frac{2\sqrt{2}}{R}\frac{1}{\sqrt{N_{i}(\lambda)(\lambda+\gamma)}\sqrt{-N_{i}(\nu)(\nu+\gamma)}\sqrt{I_{2,i}}}$
(8) $\displaystyle{\rm
where}~{}N_{i}(\tau)=G(\tau)-\frac{I_{2,i}}{\tau+\alpha}-\frac{I_{3,i}}{\tau+\gamma}+E_{i}$
(9)

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

$\rho({\bm{x}})=\sum_{i=1}^{N}c_{i}\rho_{orb,i}({\bm{x}})$ (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 ${\bm{x}}_{j}$ is expressed as

$P_{ij}=\frac{c_{i}\rho_{orb,i}({\bm{x}}_{j})}{\sum_{k=1}^{N}c_{k}\rho_{orb,k}({\bm{x}}_{j})}$
(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 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})$.

$c_{i}\propto\frac{1}{\rho_{orb,i}({\bm{x}}_{i})}$ (12)

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.

$\displaystyle\langle v_{\phi}$
$\displaystyle({\bm{x}})\rangle=\frac{1}{\rho({\bm{x}})}\sum_{i}c_{i}\rho_{orb,i}({\bm{x}})v_{\phi,i}({\bm{x}})$
(13) $\displaystyle{\rm
where}~{}v_{\phi,i}({\bm{x}})=\pm\sqrt{\frac{2I_{2,i}}{R^{2}}}$ (14)

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

##  3 Results

In what follows, we show our results on the meridian plane using the
coordinates ($r,\theta$) and also the cylindrical coordinates ($R,z$). The
relation between these coordinates is $R=r\cos\theta,z=r\sin\theta$, so
$(r,\theta)$ correspond to the 2-dimensional polar coordinates where $\theta$
is measured from the Galactic plane.

###  3.1 The entire halo sample

####  3.1.1 Density distribution

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:

$\rho(\bm{x})=\left[\frac{R_{c}^{2}}{R^{2}+(z/q)^{2}}\right]^{\alpha/2}$ (15)

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

Figure 8: The density distribution along the Galactic plane, extracted from
Figure 7. The solid lines show the reproduced distribution and the dashed
lines present the fitting curves with a best fitted power-law slope, $\alpha$.
The colors of lines correspond to the range of metallicity. The light colors
indicate the region at $R<8$ kpc, which are excluded in the fitting to the
profile $\rho\propto R^{-\alpha}$. Each line is shifted vertically for the
purpose of comparing between different metallicity ranges. The error bars show
the uncertainties propagating from the observational errors associated with
the position on the sky, parallax, proper motion and line-of-sight velocity.
The error bars in the following figures also include all of these
observational errors.

Figure 9: The density distribution along $\theta=0^{\circ}-90^{\circ}$ at
$r=8,15,25$ kpc (solid lines) and the fitted curves with the ellipsoidal
profile, Eq. (15) (dashed lines) for each metallicity range. Each curve at a
different $r$ is shifted vertically for the comparison between different
metallicity ranges.

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.

Figure 10: The axial ratio $q$ of the ellipsoidal profile (15) as a function
of the Galactocentric distance $r$ for each metallicity.

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

$\Delta(r)=\sqrt{\frac{1}{N_{\theta}}\sum_{\theta_{i}}\frac{[\log_{10}\rho_{\rm
obs}(r,\theta_{i})-\log_{10}\rho_{\rm
fit}(r,\theta_{i})]^{2}}{\sigma^{2}_{\log_{10}[\rho_{\rm
obs}(r,\theta_{i})]}}}$ (16)

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$.

Figure 11: The degree of derivation, $\Delta$, defined in Eq. (16), between
the reproduced density distribution (10) and the fitting ellipsoidal profile
(15), normalized by each uncertainty.

Figure 12:  The global distribution of $\left<v_{\phi}\right>$, defined in Eq.
(13) and (14), for the halo sample with $-1.2<{\rm[Fe/H]}<-1.0$ (top),
$-1.4<{\rm[Fe/H]}<-1.2$ (middle), and $-1.6<{\rm[Fe/H]}<-1.4$ (bottom panel).
As is clear, there exists a region with disk-like kinematics (as realized from
red color) corresponding to the metal-weal thick disk.

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.

####  3.1.2 Velocity Distribution

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: The mean rotational velocity along the Galactic plane for each
metallicity range. The dotted line indicates the zero velocity.

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.

Figure 14: The contour of the global density distribution of the GSE-like
subsample. The normalization and the metallicity ranges are same as in Figure
7.

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.

###  3.2 The GSE-like subsample

We carry out the same analysis for the subsample stars showing GSE-like
kinematics. The derived density distributions for this subsample are shown in
Figure 14, for the different metallicity ranges, which are the same as in
Figure 7. We fit this distribution to the ellipsoidal profile given in Eq.
(15).

Figure 15 shows that the density distribution along the Galactic plane (solid
lines), from which we can determine the slope $\alpha$ given in Eq. (15) for
the ellipsoidal profile (dashed lines). We note that in Figure 15 there are
discontinuities at $R\sim 8$ kpc, i.e., at the same position in Figure 8.
Since these discontinuities can be due to the observational bias, as is the
case in the entire halo sample, we do not use these inner radii to determine
the value of $\alpha$.

The determined values of $\alpha$ are $2.94^{+0.16}_{-0.05}$,
$2.89^{+0.13}_{-0.00}$, $3.28^{+0.15}_{-0.03}$, and $3.47^{+0.24}_{-0.08}$,
respectively, from metal-rich to metal-poor range. Compared with the case of
the entire halo sample, the slopes in relatively metal-rich range,
${\rm[Fe/H]}>-1.8$, are systematically smaller, and those in more metal-poor
range, ${\rm[Fe/H]}<-1.8$, are almost the same. This result is consistent with
the fact that GSE appears only in a relatively metal-rich range (Belokurov et
al., 2018). The small $\alpha$ means that the density distribution of the GSE-
like subsample decreases with Galactocentric distance more slowly the general
stellar halo.

It is worth noting in Figure 15 that the derived density distribution falls
short of a single power-law distribution at around $R=20$ kpc in the
relatively metal-rich range, ${\rm[Fe/H]}>-1.8$. This feature does not appear
in the entire halo sample shown in Figure 8. We note that we do not
artificially remove the region at $R>20$ kpc when we fit the global
distribution and determine $\alpha$. Thus we conclude that the density
distribution of the GSE-like sample shows a drop at around 20 kpc, which may
correspond to the GSE’s radial boundary (see Subsection 4.2).

To clarify that the drops are real properties of GSE-like subsample, we
calculate the value indicating the deviation of the reconstructed density
distribution from the ellipsoidal profile, defined as

$\Delta^{2}_{\theta}(r)=\frac{[\log_{10}\rho_{\rm
obs}(r,\theta)-\log_{10}\rho_{\rm
fit}(r,\theta)]^{2}}{\sigma^{2}_{\log_{10}[\rho_{\rm obs}(r,\theta)]}}$ (17)

This is the same as the factors of the sum in Eq. (16). We calculate both
$\Delta_{\theta=0^{\circ}}(r)$ for entire halo sample and GSE-like subsample
and compare them in Figure 16. We find that $\Delta_{\theta=0^{\circ}}(r)$ of
GSE-like subsample in -1.8$<$[Fe/H]$<$-1.0 gets a remarkably large value at
$10\lesssim r\lesssim 15$ kpc and $r\gtrsim 20$ kpc. This figure also suggests
that the GSE-like subsample in -1.8$<$[Fe/H]$<$-1.0 range does not follow the
power-law profile at $r\gtrsim 20$ kpc, so this is GSE’s characteristic
length. On the other hand, the reason of large $\Delta_{\theta}(r)$ at
$10\lesssim r\lesssim 15$ kpc is due to the non-ellipsoidal shape like boxy or
peanut-like shape, shown in Figure 11.

Figure 15: The distribution of the GSE-like subsample along the Galactic plane
(solid lines) and their fitting curves with a power law (dashed lines), to be
compared with Figure 8 for the entire halo sample. This comparison indicates
that the density distribution of the GSE-like subsample shows a noticeable
drop at $R\sim 20$ kpc in the metallicity range of ${\rm[Fe/H]}>-1.8$.  Figure
16: The degree of the deviation of the reconstructed density distribution from
the ellipsoidal fitting profile along the Galactic plane,
$\Delta_{\theta=0^{\circ}}(r)$ defined by Eq. (17).

Next, we determine the value of the axial ratio $q$ for this subsample. The
result is shown in Figure 17, compared with the case of the entire halo
sample. We find that in the relatively metal-rich range ${\rm[Fe/H]}>-1.8$,
the axial ratios of the GSE-like subsample are systematically larger than
those of the general stellar halo at $r<15$ kpc, whereas $q$ of the GSE-like
subsample are smaller beyond $r\sim 20$ kpc. This is consistent with the
Figure 15, implying the presence of the radial boundary of the GSE-like
members as 20 kpc. On the other hand, the axial ratios for more metal-poor
ranges are larger than or almost the same as compared with ones of the general
stellar halo. Their larger axial ratios may come from the dynamically
selection $L_{z}\sim 0$.

Finally, we remark on the $\left<v_{\phi}\right>$ distribution of this GSE-
like sample. We also find no remarkable structures like the MWTD even in
metal-rich ranges, because of the pre-selection of this sample with $L_{z}$
around 0.

##  4 Discussion

Derived global structures of the current halo sample based on the
superposition of their orbits suggest several implications for the formation
and evolution histories of the MW’s stellar halo. We discuss here what
insights are inferred from the current analysis in comparison with recent
numerical simulations for the formation of stellar halos.

###  4.1 Global halo structure from the entire sample

####  4.1.1 Implication from axial ratios and radial profiles

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$.

Figure 17: The axial ratios $q$ of the GSE-like subsample compared with those
of the entire halo sample as a function of the Galactocentric distance $r$ for
each metallicity range.

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.

####  4.1.2 Implication from a boxy/peanut-shape configuration

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.

####  4.1.3 On the discontinuities by the observational bias

Figure 18: The distribution of sample stars in the phase space
$(\sqrt{2I_{3}},E)$. The color of the points indicates the vertical angular
momentum $|L_{z}|$ of stars; red: $0.0<|L_{z}|<0.5$, yellow:
$0.5<|L_{z}|<1.0$, green: $1.0<|L_{z}|<1.5$, blue: $|L_{z}|>1.5$ ( in units of
$L_{z}$ is $\times 10^{3}{\rm kpc\ km\ s^{-1}}$). The shaded region is the
area, where the stars stay only in ”the inner region” relative to the survey
volume, depending on the values of $|L_{z}|$ depicted with the color-coded
lines: such stars having $|L_{z}|$ larger than a specific value locate below
the corresponding colored-line within the shaded region.

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.

Figure 19: The phase-space distribution of sample stars as in Fig. 18, but the
color of the points corresponds to the different range for $|L_{z}|$; black:
$0.0<|L_{z}|<2.0$, red: $2.0<|L_{z}|<3.0$, pink: $3.0<|L_{z}|<4.0$, yellow:
$|L_{z}|>4.0$ (in units of $L_{z}$ is $\times 10^{3}\ {\rm kpc\ km\ s^{-1}}$).
The shaded region is the area, where the stars staying only in ”the outer
region” relative to the survey volume, depending on the values of $|L_{z}|$
depicted with the color-coded lines: such stars having $|L_{z}|$ larger than a
specific value locate below the corresponding-colored line within the shaded
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.

###  4.2 Global structure from the GSE-like subsample

As shown in Figure 15, the distributions of the GSE-like subsample in
$-1.8<{\rm[Fe/H]}<-1.0$ fall short from a single power-law profile at $R\sim
20$ kpc on the Galactic plane. This metallicity range is just consistent with
the reported metallicities for the likely GSE member stars (Belokurov et al.,
2018; Koppelman et al., 2019). By combining with the result in Fig. 16, we
regard 20 kpc as GSE’s likely radial scale.

Figure 20: The degree of the deviation $\Delta_{\theta}(r)$ of GSE-like sample
along $\theta=0^{\circ},30^{\circ},45^{\circ}$ and $60^{\circ}$, which
correspond to solid, dashed, dashed-dotted, and dotted line, respectively.

The result that GSE’s scale is 20 kpc is consistent with some previous works.
Naidu et al. (2020) showed the occupancy fraction of each substructure as a
function of position based on the observation of H3 survey, and revealed that
GSE’s fraction peaks at $r\sim 20$ kpc and rapidly declines with increasing
$r$. Deason et al. (2018) found the sausage-like structure in the stellar halo
has its apocenter at $25-30$ kpc. Lancaster et al. (2019) also indicated that
GSE’s density rapidly decreases at 20 kpc.

This scale is considered to indicate the position of the apocenter of GSE’s
progenitor. Because it is thought that GSE is a trace of a merger event
between the MW and a large dwarf galaxy about 10 Gyr ago, the radius 20 kpc
should be the characteristic scale of the event. That “characteristic scale”
is typically expected the apocenter of the merged galaxy (Deason et al., 2018;
Lancaster et al., 2019).

Furthermore, the distance of 20 kpc is consistent with the end where the inner
halo is dominant (Carollo et al., 2007, 2010), so our result suggests that GSE
is the dominant substructure within the inner halo. This inference is
supported by the behavior in $r=10-15$ kpc of the GSE-like sample. We have
mentioned that the large values of $\Delta(r)$ might be the feature of the
inner halo (Section 4.1.2), and Fig. 16 shows that
$\Delta_{\theta=0^{\circ}}(r)$ in $r=10-15$ kpc and $-1.8<{\rm[Fe/H]}<-1.0$ of
the GSE-like sample is much larger than one of the entire halo sample. Namely,
the GSE-like sample reflects the basic properties of the inner halo
significantly.

We also attempt to find the characteristic inclination $\theta$ of the GSE-
like sample. Fig. 20 shows the deviation measure $\Delta_{\theta}(r)$ of the
GSE-like sample for various $\theta$. It shows that three lines except for
$\theta=0^{\circ}$ are similar and have no remarkable difference for each
metallicity range. Exceptionally, $\theta=30^{\circ}$ line of
$-2.2<{\rm[Fe/H]}<-1.8$ has a sharp peak at $R\sim 14$ kpc, but it is likely
to be an outlier. Thus, the GSE’s characteristic inclination is not available
here.

Next, we have shown that the slope $\alpha$ of the GSE-like subsample in
$-1.8<{\rm[Fe/H]}<-1.0$ is 2.9, being systematically smaller than slopes of
the general halo, $\alpha=3.2-3.4$. On the other hand, in more metal-poor
ranges ${\rm[Fe/H]}<-1.8$, the values of $\alpha$ are almost the same for both
samples. As shown in Subsection 4.1.1, it is known that accreted halo
components have smaller $\alpha$ with $\alpha<3.0$ (Rodriguez-Gomez et al.,
2016; Monachesi et al., 2019). Thus smaller $\alpha$ of GSE-like subsample in
our work is consistent with its origin.

Additionally, Figure 17 shows that the axial ratios of the GSE-like subsample
in $-1.8<{\rm[Fe/H]}<-1.0$ are systematically larger than the general stellar
halo within $r\lesssim 20$ kpc, suggesting that GSE’s shape is closer to a
sphere than the general halo. This property of the GSE, which is a merging
origin, was also suggested in the simulation work by Monachesi et al. (2019).
It showed that a stellar halo which experienced accretion of a satellite
recently forms more spherical shape than a general halo.

###  4.3 The completeness of the observational data

This work analyzes the global distribution of the MW’s stellar halo by using
local observational data. There might be some concern about the effect of the
completeness in the result. We carry out the additional analysis below to get
insights into this issue.

Figure 21: The result for the axial ratio of subsamples of $D<1$ kpc (upper
panel) and $D<2$ kpc (lower panel), compared with the main sample (dark-
colored lines, which are the same as in Fig. 10.)

We create other two subsamples, whose additional selection condition is that
the heliocentric distance, $D$, is smaller than 1 kpc and 2 kpc, respectively.
We carry out the same analysis for these subsamples.

The axial ratio of these subsample is shown in Fig. 21 as the representative
result. In the case of $D<1$ kpc the reconstructed distribution is very bumpy
with large uncertainty (the upper panel of Fig. 21), so this result has poor
reliability. Two reasons of the poor reliability are considered. The first
reason is that the $D<1$ kpc sample has too few small number of stars (1,895
stars, 7.6% of the main sample) to reproduce the global distribution. The
second is that the region within 1 kpc from the Sun is too small to represent
a global distribution. In any case, this result has insufficient completeness.

On the other hand, for $D<2$ kpc subsample, we obtain the result almost the
same as the main result (the lower panel of Fig. 21). It is suggested that the
completeness of $D<2$ kpc subsample is sufficient to reproduce the global
distribution up to $\sim 30$ kpc. Therefore, it suggests that our main result,
using stellar $D\lesssim 4$ kpc, is insensitive to the sample completeness.

##  5 Conclusion

We have constructed and analyzed the global distribution of the MW’s stellar
halo using Gaia EDR3 combined with SDSS DR16 SEGUE catalogs. For this purpose,
we adopt the method developed by Sommer-Larsen & Zhen (1990) following May &
Binney (1986)’s strategy of reconstructing the global halo based on the
superposition of stars’ orbits. We calculate the orbital density of each star
and assign a proper weight to it through a maximum likelihood approach. We
then sum up all of the orbital densities with these weights and derive the
global stellar distribution from observational sample stars currently located
within a local volume near the Sun.

We have found that more metal-rich halo stars show more flattened distribution
as have been derived in previous studies using much smaller numbers of sample
stars (e.g. Carollo et al., 2007). In the metal-rich range of
$-1.4<{\rm[Fe/H]}<-1.0$, we have identified, through the reconstruction of
global rotational velocities, $\left<v_{\phi}\right>$, the presence of the
MWTD at $z<4$ kpc and $R<15$ kpc. The detailed configuration of the global
halo structures in relatively metal-rich ranges of ${\rm[Fe/H]}>-1.8$ is found
to have boxy/peanut-shape rather than a simple ellipsoidal shape and this may
reflect a past merging event. On the other hand, in most metal-poor range of
${\rm[Fe/H]}<-2.2$, the halo structure is close to a spherical shape.

We have selected the subsample of the GSE-like stars with $E>-150,000~{}{\rm
km^{2}~{}s^{-2}}$ and $-500<L_{z}<500~{}{\rm kpc~{}km~{}s^{-1}}$. In contrast
to the entire halo sample, the global density profile of these GSE-like stars
show a noticeable drop at $r\sim 20$ kpc over $0^{\circ}<\theta<45^{\circ}$,
implying a apocentric boundary of debris stars from a merging progenitor
galaxy of the GSE. Additionally, the behaviors of $\alpha$ and $q$ of GSE-like
subsample indicate the features of accreted halo components: its density
distribution decreases more slowly with Galactocentric radius and has a more
spherical shape over $r<20$ kpc than the general stellar halo sample.

We are grateful to the referee for invaluable comments, which help revise the
current paper. We also thank Daniela Carollo for her many useful comments on
the original manuscript. We acknowledge support from the JSPS and MEXT Grant-
in-Aid for Scientific Research (No. 17H01101, 18H04434 and 18H05437). Funding
for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan
Foundation, the U.S. Department of Energy Office of Science, and the
Participating Institutions. SDSS-IV acknowledges support and resources from
the Center for High Performance Computing at the University of Utah. The SDSS
website is www.sdss.org. SDSS-IV is managed by the Astrophysical Research
Consortium for the Participating Institutions of the SDSS Collaboration
including the Brazilian Participation Group, the Carnegie Institution for
Science, Carnegie Mellon University, Center for Astrophysics — Harvard &
Smithsonian, the Chilean Participation Group, the French Participation Group,
Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli
Institute for the Physics and Mathematics of the Universe (IPMU) / University
of Tokyo, the Korean Participation Group, Lawrence Berkeley National
Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-
Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik
(MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE),
National Astronomical Observatories of China, New Mexico State University, New
York University, University of Notre Dame, Observatário Nacional / MCTI, The
Ohio State University, Pennsylvania State University, Shanghai Astronomical
Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma
de México, University of Arizona, University of Colorado Boulder, University
of Oxford, University of Portsmouth, University of Utah, University of
Virginia, University of Washington, University of Wisconsin, Vanderbilt
University, and Yale University. This work has made use of data from the
European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia),
processed by the Gaia Data Processing and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

## References

  * Anguiano et al. (2020) Anguiano, B., Majewski, S. R., Hayes, C. R., et al. 2020, AJ, 160, 43, doi: 10.3847/1538-3881/ab9813
  * Antoja et al. (2018) Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, Nature, 561, 360, doi: 10.1038/s41586-018-0510-7
  * Athanassoula (2005) Athanassoula, E. 2005, MNRAS, 358, 1477, doi: 10.1111/j.1365-2966.2005.08872.x
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147, doi: 10.3847/1538-3881/abd806
  * Batsleer & Dejonghe (1994) Batsleer, P., & Dejonghe, H. 1994, A&A, 287, 43
  * Beane et al. (2019) Beane, A., Sanderson, R. E., Ness, M. K., et al. 2019, ApJ, 883, 103, doi: 10.3847/1538-4357/ab3d3c
  * Beers et al. (2002) Beers, T. C., Drilling, J. S., Rossi, S., et al. 2002, AJ, 124, 931, doi: 10.1086/341377
  * Bekki & Chiba (2000) Bekki, K., & Chiba, M. 2000, ApJ, 534, L89, doi: 10.1086/312636
  * Bekki & Chiba (2001) —. 2001, ApJ, 558, 666, doi: 10.1086/322300
  * Bell et al. (2006) Bell, E. F., Naab, T., McIntosh, D. H., et al. 2006, ApJ, 640, 241, doi: 10.1086/499931
  * Belokurov et al. (2018) Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611, doi: 10.1093/mnras/sty982
  * Benson et al. (2004) Benson, A. J., Lacey, C. G., Frenk, C. S., Baugh, C. M., & Cole, S. 2004, MNRAS, 351, 1215, doi: 10.1111/j.1365-2966.2004.07870.x
  * Bignone et al. (2019) Bignone, L. A., Helmi, A., & Tissera, P. B. 2019, ApJ, 883, L5, doi: 10.3847/2041-8213/ab3e0e
  * Bond et al. (2010) Bond, N. A., Ivezić, Ž., Sesar, B., et al. 2010, ApJ, 716, 1, doi: 10.1088/0004-637X/716/1/1
  * Bovy (2015) Bovy, J. 2015, ApJS, 216, 29, doi: 10.1088/0067-0049/216/2/29
  * Bullock & Johnston (2005) Bullock, J. S., & Johnston, K. V. 2005, ApJ, 635, 931, doi: 10.1086/497422
  * Busso et al. (2021) Busso, G., Cacciari, C., Bellazzini, M., et al. 2021, Gaia EDR3 documentation Chapter 5: Photometric data, Gaia EDR3 documentation
  * Carollo & Chiba (2021) Carollo, D., & Chiba, M. 2021, ApJ, 908, 191, doi: 10.3847/1538-4357/abd7a4
  * Carollo et al. (2007) Carollo, D., Beers, T. C., Lee, Y. S., et al. 2007, Nature, 450, 1020, doi: 10.1038/nature06460
  * Carollo et al. (2010) Carollo, D., Beers, T. C., Chiba, M., et al. 2010, ApJ, 712, 692, doi: 10.1088/0004-637X/712/1/692
  * Carollo et al. (2019) Carollo, D., Chiba, M., Ishigaki, M., et al. 2019, ApJ, 887, 22, doi: 10.3847/1538-4357/ab517c
  * Chiba & Beers (2000) Chiba, M., & Beers, T. C. 2000, AJ, 119, 2843, doi: 10.1086/301409
  * Chiba & Beers (2001) —. 2001, ApJ, 549, 325, doi: 10.1086/319068
  * Conroy et al. (2019) Conroy, C., Bonaca, A., Cargile, P., et al. 2019, The Astrophysical Journal, 883, 107, doi: 10.3847/1538-4357/ab38b8
  * Cooper et al. (2015) Cooper, A. P., Parry, O. H., Lowing, B., Cole, S., & Frenk, C. 2015, MNRAS, 454, 3185, doi: 10.1093/mnras/stv2057
  * Cooper et al. (2010) Cooper, A. P., Cole, S., Frenk, C. S., et al. 2010, Monthly Notices of the Royal Astronomical Society, 406, 744–766, doi: 10.1111/j.1365-2966.2010.16740.x
  * Cordoni et al. (2021) Cordoni, G., Da Costa, G. S., Yong, D., et al. 2021, MNRAS, 503, 2539, doi: 10.1093/mnras/staa3417
  * Cox et al. (2006) Cox, T. J., Dutta, S. N., Di Matteo, T., et al. 2006, ApJ, 650, 791, doi: 10.1086/507474
  * Cui et al. (2012) Cui, X.-Q., Zhao, Y.-H., Chu, Y.-Q., et al. 2012, Research in Astronomy and Astrophysics, 12, 1197, doi: 10.1088/1674-4527/12/9/003
  * De Silva et al. (2015) De Silva, G. M., Freeman, K. C., Bland-Hawthorn, J., et al. 2015, MNRAS, 449, 2604, doi: 10.1093/mnras/stv327
  * de Zeeuw & Lynden-Bell (1985) de Zeeuw, P. T., & Lynden-Bell, D. 1985, MNRAS, 215, 713, doi: 10.1093/mnras/215.4.713
  * de Zeeuw (1985) de Zeeuw, T. 1985, MNRAS, 216, 273, doi: 10.1093/mnras/216.2.273
  * Deason et al. (2018) Deason, A. J., Belokurov, V., Koposov, S. E., & Lancaster, L. 2018, ApJ, 862, L1, doi: 10.3847/2041-8213/aad0ee
  * Deason et al. (2019) Deason, A. J., Fattahi, A., Belokurov, V., et al. 2019, MNRAS, 485, 3514, doi: 10.1093/mnras/stz623
  * Deason et al. (2021) Deason, A. J., Erkal, D., Belokurov, V., et al. 2021, MNRAS, 501, 5964, doi: 10.1093/mnras/staa3984
  * Dejonghe & de Zeeuw (1988) Dejonghe, H., & de Zeeuw, T. 1988, ApJ, 333, 90, doi: 10.1086/166727
  * Eadie & Jurić (2019) Eadie, G., & Jurić, M. 2019, ApJ, 875, 159, doi: 10.3847/1538-4357/ab0f97
  * Eggen et al. (1962) Eggen, O. J., Lynden-Bell, D., & Sandage, A. R. 1962, ApJ, 136, 748, doi: 10.1086/147433
  * Evans et al. (2020) Evans, T. A., Fattahi, A., Deason, A. J., & Frenk, C. S. 2020, MNRAS, 497, 4311, doi: 10.1093/mnras/staa2202
  * Fernández-Trincado et al. (2020) Fernández-Trincado, J. G., Beers, T. C., & Minniti, D. 2020, A&A, 644, A83, doi: 10.1051/0004-6361/202039434
  * Feuillet et al. (2020) Feuillet, D. K., Feltzing, S., Sahlholdt, C. L., & Casagrande, L. 2020, MNRAS, 497, 109, doi: 10.1093/mnras/staa1888
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2016, A&A, 595, A2, doi: 10.1051/0004-6361/201629512
  * Gaia Collaboration et al. (2018) —. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gaia Collaboration et al. (2021) —. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Gillessen et al. (2009) Gillessen, S., Eisenhauer, F., Trippe, S., et al. 2009, ApJ, 692, 1075, doi: 10.1088/0004-637X/692/2/1075
  * Grand et al. (2019) Grand, R. J. J., Deason, A. J., White, S. D. M., et al. 2019, MNRAS, 487, L72, doi: 10.1093/mnrasl/slz092
  * Hagen et al. (2019) Hagen, J. H. J., Helmi, A., de Zeeuw, P. T., & Posti, L. 2019, A&A, 629, A70, doi: 10.1051/0004-6361/201935264
  * Helmi et al. (2018) Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Nature, 563, 85, doi: 10.1038/s41586-018-0625-x
  * Ibata et al. (2003) Ibata, R. A., Irwin, M. J., Lewis, G. F., Ferguson, A. M. N., & Tanvir, N. 2003, Monthly Notices of the Royal Astronomical Society, 340, L21–L27, doi: 10.1046/j.1365-8711.2003.06545.x
  * Iorio & Belokurov (2019) Iorio, G., & Belokurov, V. 2019, MNRAS, 482, 3868, doi: 10.1093/mnras/sty2806
  * Ivezić et al. (2008) Ivezić, Ž., Sesar, B., Jurić, M., et al. 2008, ApJ, 684, 287, doi: 10.1086/589678
  * Jurić et al. (2008) Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864, doi: 10.1086/523619
  * Kim & Lépine (2021) Kim, B., & Lépine, S. 2021, MNRAS, doi: 10.1093/mnras/stab3671
  * Kim et al. (2021) Kim, Y. K., Lee, Y. S., Beers, T. C., & Koo, J.-R. 2021, ApJ, 911, L21, doi: 10.3847/2041-8213/abf35e
  * Koppelman et al. (2020) Koppelman, H. H., Bos, R. O. Y., & Helmi, A. 2020, A&A, 642, L18, doi: 10.1051/0004-6361/202038652
  * Koppelman et al. (2021) Koppelman, H. H., Hagen, J. H. J., & Helmi, A. 2021, A&A, 647, A37, doi: 10.1051/0004-6361/202039390
  * Koppelman et al. (2019) Koppelman, H. H., Helmi, A., Massari, D., Price-Whelan, A. M., & Starkenburg, T. K. 2019, A&A, 631, L9, doi: 10.1051/0004-6361/201936738
  * Kruijssen et al. (2020) Kruijssen, J. M. D., Pfeffer, J. L., Chevance, M., et al. 2020, MNRAS, 498, 2472, doi: 10.1093/mnras/staa2452
  * Lancaster et al. (2019) Lancaster, L., Koposov, S. E., Belokurov, V., Evans, N. W., & Deason, A. J. 2019, MNRAS, 486, 378, doi: 10.1093/mnras/stz853
  * Li et al. (2020) Li, H., Du, C., Yang, Y., et al. 2020, ApJ, 895, 23, doi: 10.3847/1538-4357/ab8733
  * Liang et al. (2021) Liang, X.-L., Chen, Y.-Q., Zhao, J.-K., & Zhao, G. 2021, Research in Astronomy and Astrophysics, 21, 128, doi: 10.1088/1674-4527/21/5/128
  * Mackereth et al. (2019) Mackereth, J. T., Schiavon, R. P., Pfeffer, J., et al. 2019, MNRAS, 482, 3426, doi: 10.1093/mnras/sty2955
  * Majewski et al. (2017) Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, The Astronomical Journal, 154, 94, doi: 10.3847/1538-3881/aa784d
  * May & Binney (1986) May, A., & Binney, J. 1986, MNRAS, 221, 857, doi: 10.1093/mnras/221.4.857
  * McCarthy et al. (2012) McCarthy, I. G., Font, A. S., Crain, R. A., et al. 2012, MNRAS, 420, 2245, doi: 10.1111/j.1365-2966.2011.20189.x
  * Monachesi et al. (2019) Monachesi, A., Gómez, F. A., Grand, R. J. J., et al. 2019, MNRAS, 485, 2589, doi: 10.1093/mnras/stz538
  * Monty et al. (2020) Monty, S., Venn, K. A., Lane, J. M. M., Lokhorst, D., & Yong, D. 2020, MNRAS, 497, 1236, doi: 10.1093/mnras/staa1995
  * Morganson et al. (2016) Morganson, E., Conn, B., Rix, H.-W., et al. 2016, ApJ, 825, 140, doi: 10.3847/0004-637X/825/2/140
  * Mróz et al. (2019) Mróz, P., Udalski, A., Skowron, D. M., et al. 2019, ApJ, 870, L10, doi: 10.3847/2041-8213/aaf73f
  * Myeong et al. (2018a) Myeong, G. C., Evans, N. W., Belokurov, V., Sanders, J. L., & Koposov, S. E. 2018a, ApJ, 863, L28, doi: 10.3847/2041-8213/aad7f7
  * Myeong et al. (2018b) —. 2018b, ApJ, 856, L26, doi: 10.3847/2041-8213/aab613
  * Myeong et al. (2019) Myeong, G. C., Vasiliev, E., Iorio, G., Evans, N. W., & Belokurov, V. 2019, MNRAS, 488, 1235, doi: 10.1093/mnras/stz1770
  * Naab et al. (2006) Naab, T., Khochfar, S., & Burkert, A. 2006, ApJ, 636, L81, doi: 10.1086/500205
  * Naidu et al. (2020) Naidu, R. P., Conroy, C., Bonaca, A., et al. 2020, ApJ, 901, 48, doi: 10.3847/1538-4357/abaef4
  * Naidu et al. (2021) —. 2021, arXiv e-prints, arXiv:2103.03251.  https://arxiv.org/abs/2103.03251
  * Newberg et al. (2002) Newberg, H. J., Yanny, B., Rockosi, C., et al. 2002, ApJ, 569, 245, doi: 10.1086/338983
  * Perryman et al. (1997) Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, A&A, 500, 501
  * Posti & Helmi (2019) Posti, L., & Helmi, A. 2019, A&A, 621, A56, doi: 10.1051/0004-6361/201833355
  * Purcell et al. (2010) Purcell, C. W., Bullock, J. S., & Kazantzidis, S. 2010, MNRAS, 404, 1711, doi: 10.1111/j.1365-2966.2010.16429.x
  * Rodriguez-Gomez et al. (2016) Rodriguez-Gomez, V., Pillepich, A., Sales, L. V., et al. 2016, MNRAS, 458, 2371, doi: 10.1093/mnras/stw456
  * Santistevan et al. (2020) Santistevan, I. B., Wetzel, A., El-Badry, K., et al. 2020, MNRAS, 497, 747, doi: 10.1093/mnras/staa1923
  * Sestito et al. (2019) Sestito, F., Longeard, N., Martin, N. F., et al. 2019, MNRAS, 484, 2166, doi: 10.1093/mnras/stz043
  * Sofue (2020) Sofue, Y. 2020, Galaxies, 8, 37, doi: 10.3390/galaxies8020037
  * Sommer-Larsen & Zhen (1990) Sommer-Larsen, J., & Zhen, C. 1990, MNRAS, 242, 10, doi: 10.1093/mnras/242.1.10
  * Steinmetz et al. (2006) Steinmetz, M., Zwitter, T., Siebert, A., et al. 2006, AJ, 132, 1645, doi: 10.1086/506564
  * Tissera et al. (2013) Tissera, P. B., Scannapieco, C., Beers, T. C., & Carollo, D. 2013, MNRAS, 432, 3391, doi: 10.1093/mnras/stt691
  * White & Rees (1978) White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341, doi: 10.1093/mnras/183.3.341
  * Wyse (2019) Wyse, R. 2019, Science, 365, 979, doi: 10.1126/science.aay0628
  * Yanny et al. (2003) Yanny, B., Newberg, H. J., Grebel, E. K., et al. 2003, The Astrophysical Journal, 588, 824–841, doi: 10.1086/374220
  * Yanny et al. (2009) Yanny, B., Rockosi, C., Newberg, H. J., et al. 2009, The Astronomical Journal, 137, 4377, doi: 10.1088/0004-6256/137/5/4377
  * Zhao et al. (2012) Zhao, G., Zhao, Y.-H., Chu, Y.-Q., Jing, Y.-P., & Deng, L.-C. 2012, Research in Astronomy and Astrophysics, 12, 723, doi: 10.1088/1674-4527/12/7/002


# Near Field Optimization Algorithm for Reconfigurable Intelligent Surface

E. Colellaref1ref2  L. Bastianelli ref1ref2  F. Moglie ref1ref2 and V. Mariani
Primianiref1ref2  ref1Università Politecnica delle Marche ref2Consorzio
Nazionale Interuniversitario delle Telecomunicazioni (CNIT)

###### Abstract

Reconfigurable intelligent surface (RIS) is a type of wireless communication
technology that uses a reconfigurable surface, such as a wall or building that
is able to adjust its properties by an integrated optimization algorithm in
order to optimize the signal propagation for a given communication scenario.
As a reconfiguration algorithm the multidimensional optimization of the GNU
scientific library was analyzed to evaluate the performance of the smart
surface in the quality of signal reception. This analysis took place by means
of electrodynamic simulations based on the finite difference time domain
method. Through these simulations it was possible to observe the efficiency of
the algorithm in the reconfiguration of the RIS, managing to focus the
electromagnetic waves in a remarkable way towards the point of interest.

##  1 Introduction

A metasurface is a planar destruction of metamaterial, a man-made material
engineered to have electromagnetic properties not found in nature [1, 2]. This
material is made up of elementary units, arranged periodically in repeating
patterns with dimensions and spacings much smaller than the wavelength [4].
Thanks to these microscopic characteristics, chemical due to the nature of the
resonant elements and geometric due to the dimensional nature, the metasurface
is able to control the propagation of the electromagnetic waves with which it
interacts, manipulating the distribution of the fields in space. Their ability
to control the response of incident waves lies in the interaction of the
electric and magnetic fields of the waves with the periodic structure of the
surface which leads to resonances capable of changing the oscillatory nature
of the electromagnetic fields [7]. These interactions lead to resonant effects
that lead to a wide range of applications such as broadband focusing [6],
peculiar reflection/refraction [7], orbital angular momentum creation [8],
space–surface wave manipulations [9] and other applications in both microwave
and optical frequencies. Although the metasurface has the ability to focus
electromagnetic waves at various points in space, it remains a static
structure, allowing it to always respond to the incident wave in the same way.
Nowadays, therefore, the greatest technological interest has turned towards
dynamic structures of metasurfaces capable of dynamically controlling the
electromagnetic response of the surface. These surfaces are called
Reconfigurable Intelligent Surfaces (RISs). The RIS is a type of surface that
can change its configuration to affect signal propagation with low energy
consumption managed by optimization software for the dynamic control of the
propagation of electromagnetic waves [10]. This can be done through the use of
antenna elements which can be switched or adjusted to change the signal
reflection and transmission behaviour [14]. This allows to improve the
performance of the communication system and to obtain greater flexibility in
the design of the networks [11]. The control of electromagnetic waves takes
place through an intelligent reconfiguration of the wireless propagation
environment which allows to improve the capacity and coverage of wireless
networks [12]. Their goal is to overcome the destructive effect of multi-path
fading, which attenuates the communication signal strength, and improve the
signal quality. Therefore, the RISs appear to be an advanced form of smart
surface for the emerging hardware technology for the sixth-generation (6G) of
radio communication networks [13]. For RIS reconfiguration there are several
optimization algorithms such as Several integrated software for metasurface
reconfiguration such as Particle Swarm Optimization (PSO) [14], Network
Slicing Resource Optimization (NSRO) [15], Salp Swarm Algorithm (SSA) [16],
Gradient-based Optimization Algorithms which can be used to maximize signal
reception in space. In this paper we analyze an open source multidimensional
minimization algorithm from GNU Scientific Library (GSL) as the multivariable
optimization algorithm for the RIS reconfiguration for smart radio
environment.

##  2 Simulation Set-Up

The electromagnetic simulations are based on the analysis of the RIS focusing
performance by its dynamic re-configurations in electromagnetic simulations.
In order to analyze the EM response of the RIS reconfiguration, numerical
simulations were run according to the finite difference time domain (FDTD)
method to reproduce computational electrodynamics. The entire FDTD code has
been implemented following the standard formulation written in C-language
code. The electromagnetic simulations were performed in a domain of
140x140x200 cells of 1 mm cell size. The domain includes the RIS, two
antennas, one transmitting antenna (Tx) and one receiving antenna (Rx)
separated by a perfect electric conductor (PEC) metallic barrier. The antennas
are dipole antennas of 10 cells working at frequency band 2500-5000 MHz. The
barrier, instead, consists of 1x60x60 cells placed in the middle of the two
antennas in order to avoid the direct communication between them. The main
characteristic of the barrier is the PEC boundary condition that consists of
zero tangent electric field, that does not allow to transmit signal through
the barrier. In this way, the only contribution that will allow the signal
from the Tx to reach the Rx antenna will reside in the RIS. The simulated RIS
is a periodic structure made up of 5x5 perfectly conductive resonant units
(PECs) connected to each other by means of varactor diodes. The structure has
dimensions equal to $56\times 1\times 56$ cells, in which the patches are
separated by 1 cell. Each patch communicates with two patches via variable
capacitance varactor diodes for a total of 40 diodes in a dielectric film of 1
cell with $\epsilon_{r}=3$, separated from margins of 1 cell. This RIS is
positioned in front of the dividing barrier, representing the only
contribution for the receiving antenna. The entire configuration is shown in
Figure 1. The dielectric material of the interior of the chamber is air while
the walls of the simulation space are perfect matched layers (PML).

Figure 1: Representation of the simulated 56 $\times$ 56 mm 2 RIS. It consists
of 5 $\times$ 5 PEC patches of 1 cm 2 connected by 40 varactor diodes,
vertically and horizontally polarized.The dielectric support of
$\epsilon_{r}=3$ and 1 mm of thickness. Figure 2: Representation of the
simulation configuration: the RIS is positioned in the upper center of the
workspace, the PEC barrier in the lower center while the two Tx and Rx
antennas are located on the sides of the barrier. All the measures are in mm.

##  3 Multidimensional Minimization

The GNU Scientific Library is a set of open source numerical routines for
scientific computing written in C and C++. It includes different kind of
optimization algorithms that can be implemented in order to find the minimum
or maximum of a given function. Among the optimization algorithm available in
GSL, the multidimensional minimization has been chosen. This choice lies in
the number of vactor diodes each of which is characterized by its own
capacity. The goal of the algorithm is to iteratively find the ideal
configuration of the varactor capacities in order to have the maximum towards
the receiving antenna. The optimization function called multidimensional
minimization looks for a point where the function to be optimized
$f(x_{1},..x_{n})$ assumes the lowest value with respect to the neighboring
points. To find the maximum, however, the objective is reversed to find
${1-\text{min}(f(x_{1},..x_{n}))}$. The algorithm begins with a guess and
descends using a search technique. A one-dimensional line is minimized using
the gradient of the function until the lowest point is located within a
reasonable tolerance. The process is repeated until the true n-dimensional
minimum is found, at which point the search direction is updated using local
information from the function and its derivatives. Until the simplex is small
enough, the iteration process is continued. The process entails three steps:
setting up the minimizer state s for algorithm T, updating s with iteration T,
determining whether s has reached convergence, and repeating the iteration if
necessary. Initializing the minimizer s to minimize the function fdf, the pre-
built function initializes the multidimensional minimizer by starting from the
initial point x. The tol parameter controls the accuracy of the line
minimization, while the step-size parameter controls the size of the first
trial step. If the gradient of the function g is orthogonal to the current
search direction p and has a relative accuracy of tol, line minimization is
frequently regarded as successful, with:

$p\cdot g<tol|p||g|$ (1)

Figure 3: Electromagnetic field distribution in case of not reconfigured RIS
(on the left) and in case of configured RIS (on the right).

A tolerance value of 0.1 is chosen for the majority of tasks because the
multidimensional minimization algorithm’s line reduction process only needs to
be performed roughly. A parametric function with n variables
$f(x,\text{params})$ must be defined for the minimizers, along with two
additional routines: one for calculating the function’s gradient and the
second for calculating both the function value and the gradient. The iteration
functions are used to update the minimizer’s state, and because the same
function is used by all minimizers, different techniques can be swapped at
runtime without modifying the code. The minimization process should stop when
one of the following conditions is met: a minimum has been found within the
user-specified precision, a user-specified maximum number of iterations has
been reached, or an error has occurred. The minimization process compares the
gradient g norm to the absolute tolerance epsabs. At a minimum, the gradient
of a multidimensional function is zero. The test returns success if this
condition is met:

$|g|<epsabs$ (2)

At this point the algorithm stops and the optimum set of parameters are then
run on the last FDTD simulation for the best RIS reconfiguration.

##  4 Results

This section shows the results of the electromagnetic simulations according to
the FDTD method considering the set-up shown in Figure 1. The reported results
refer to the electromagnetic response of the RIS before and after the
metasurface reconfiguration. On the left in Figure 2, it is shown the
electromagnetic field distribution in case of not re-configured RIS. In this
case the diodes are all switched on at fixed 1 pF capacitance, and how it is
possibile to see there is not focusing toward the receiving antenna. On the
right in Figure 2, however, it is shown the electromagnetic field distribution
in case of re-configured RIS.In this case, instead, there is greater focus
near the reception point.

Figure 4: Received signal amplitude in time domain in case of not configured
RIS (on top - red line) and configured RIS (top - blue line). The maximum
value of the not configured RIS is 0.015 V while the configured one is 0.02 V.
On the bottom is reported the $s_{21}$ in frequency domain in case of not
configured RIS (red line) and configured RIS (blue line) from 2.5 GHz up to
5.0 GHz.

This distribution of fields focused towards the receiving antenna
significantly improves the received signal. To quantify the signal received in
the two configurations it is possible to observe the Figure 3. On the left,
the received signal in both reconfiguration is reported. The signal amplitude
in time domain is shown on the left. As can be seen, in case of not
reconfigured RIS (red line) the received signal amplitude is lower than in
case of reconfigured RIS (blue line). On the right, instead, is reported the
$s_{21}$ in frequency domain, in order to analyze the spectrum of the received
signal. From 2.5 GHz up to 3.25 GHz, the received signal of the reconfigured
RIS is higher as from 3.25 GHz up to 5 GHz with respect to the not
reconfigured RIS.

## Acknowledgements

This work has been supported by EU H2020 RISE-6G project under the grant
number 101017011.

##  5 Conclusion

In conclusion, the GNU Scientific Library’s multidimensional optimization
algorithm is effective at finding the optimal configuration for a specific
point of focus. The algorithm has also two variations, the Fletcher-Reeves
conjugate gradient algorithm and the Polak-Ribiere conjugate gradient
algorithm. Both methods use a sequence of search directions to approximate the
curvature of the function near the minimum. The Fletcher-Reeves method uses
one approach for the coefficient beta, while the Polak-Ribiere method uses a
different approach. Next steps foresee the analysis of the two versions of the
multidimensional algorithm to guarantee a better reconfiguration of the
metasurface. Upcoming activities are focused on analyzing different
optimization algorithms and setups for the best reconfiguration of the RIS in
smart radio environments.

## References

  * [1] C. Caloz, T. Itoh, Electromagnetic Metamaterials, Transmission Line Theory and Microwave Applications, Wiley-IEEE Press 2005
  * [2] K. Achouri, B. A. Khan, S. Gupta, G. Lavigne, M. A. Salem, and C. Caloz, “Synthesis of electromagnetic metasurfaces: Principles and illustrations,” EPJ Appl. Metamater., vol. 2, 2015, Art. no. 12.
  * [3] K. Achouri and C. Caloz, “Design, concepts and applications of electromagnetic metasurfaces,” arXiv:1712.00618, Dec. 2017.
  * [4] T. J. Cui, D. Smith, and R. Liu, Metamaterials: Theory, Design, and Applications, 1st ed. New York, NY, USA: Springer, 2009.
  * [5] C. L. Holloway, M. A. Mohamed, E. F. Kuester, and A. Dienstfrey, “Reflection and transmission properties of a metafilm: With an application to a controllable surface composed of resonant particles,” IEEE Trans. Electromagn. Compat., vol. 47, no. 4, pp. 853–865, Nov. 2005.
  * [6] A. Pors, M. G. Nielsen, R. L. Eriksen, and S. I. Bozhevolnyi, “Broadband focusing flat mirrors based on plasmonic gradient metasurfaces,” Nano Lett., vol. 13, no. 2, pp. 829–834, 2013.
  * [7] N. Yu and F. Capasso, “Flat optics with designer metasurfaces,” Nature Mater., vol. 13, no. 2, pp. 139–150, 2014.
  * [8] E. Karimi, S. A. Schulz, I. De Leon, H. Qassim, J. Upham, and R. W. Boyd, “Generating optical orbital angular momentum at visible wavelengths using a plasmonic metasurface,” Light, Sci. Appl., vol. 3, no. 5, p. e167, 2014.
  * [9] K. Achouri and C. Caloz. (Dec. 2016). “Space-wave routing via surface waves using a metasurface system.” https://arxiv.org/abs/1612.05576
  * [10] E. Basar, M. Di Renzo, J. De Rosny, M. Debbah, M. S. Alouini, and R. Zhang, “Wireless communications through reconfigurable intelligent surfaces,” IEEE Access, vol. 7, pp. 116753-116773, 2019.
  * [11] Y.-C. Liang, R. Long, Q. Zhang, J. Chen, H. Victor Cheng, and H. Guo, ”Large intelligent Surface Antennas (LISA): Making reflective radios smart,” 2019, http://arxiv.org/abs/1906.06578.
  * [12] W. Tang, X. Li, J. Y. Dai, S. Jin, Y. Zeng, Q. Cheng, and T. J. Cui, ”Wire- less communications with programmable metasurface: Transceiver design and experimental results,” China Commun., vol. 16, no. 5, pp. 46–61, May 2019.
  * [13] W. Tang, M. Zheng Chen, J. Yan Dai, Y. Zeng, X. Zhao, S. Jin, Q. Cheng, and T. Jun Cui, “Wireless communications with programmable metasurface: New paradigms, opportunities, and challenges on transceiver design,” 2019, http://arxiv.org/abs/1907.01956
  * [14] Wang, D., Tan, D. Liu, L. “Particle swarm optimization algorithm: an overview”. Soft Comput 22, 387–408 (2018). https://doi.org/10.1007/s00500-016-2474-6
  * [15] X. Shen, Z. Zeng, X. Liu,f “RIS-Assisted Network Slicing Resource Optimization Algorithm for Coexistence of eMBB and URLLC,” Electronics 2022, 11, 2575. https://doi.org/10.3390/electronics11162575
  * [16] E. Hegazy, M.A. Makhlouf, Gh. S. El-Tawel, “Improved salp swarm algorithm for feature selection”, Journal of King Saud University - Computer and Information Sciences, https://doi.org/10.1016/j.jksuci.2018.06.003.