††thanks<EMAIL_ADDRESS>
# Effect of the Large Magellanic Cloud on the kinematics of Milky Way
satellites
and virial mass estimate

Andrey Kravtsov 1,2,3,⋆ 0000-0003-4307-634X Sophia Winney 1
0009-0004-9634-4808 1Department of Astronomy & Astrophysics, The University of
Chicago, Chicago, IL 60637 USA  2Kavli Institute for Cosmological Physics, The
University of Chicago, Chicago, IL 60637 USA  3Enrico Fermi Institute, The
University of Chicago, Chicago, IL 60637 USA

###### Abstract

We present a study illustrating the effects of the passage of a Large
Magellanic Cloud (LMC) mass satellite on the distance and velocity
distributions of satellites in $\Lambda+$Cold Dark Matter simulations of Milky
Way (MW) sized halos. In agreement with previous studies, we find that during
such a passage the velocity distribution develops a high-velocity tail due to
the reflex motion of the inner part of the halo, which can bias velocity-based
virial halo mass estimates. When the velocity distribution of MW satellites is
corrected for effects of the LMC passage, it is consistent with the
distributions in halos of masses as low as $M_{\rm 200c}=8\times
10^{11}\,\,M_{\odot}$ and as high as $1.5\times 10^{12}\,\,M_{\odot}$. We
present a new halo mass estimator $M_{\rm 200c}=c\sigma^{2}_{\rm 3D}r_{\rm
med}$, where $c$ is a constant calibrated using satellite systems in the
simulated MW-sized halos, $\sigma^{2}_{\rm 3D}$ is the variance of 3D
velocities taken with the sign of the radial velocity of each satellite, and
$r_{\rm med}$ is the median halocentric distance of the satellites. We show
that the estimator has only $s=8\%$ scatter around the median relation of the
estimated and true halo masses and deviates by $<2s$ from the median during
the pericentric passage of an LMC-like subhalo. This is because
$\sigma^{2}_{\rm 3D}$ and $r_{\rm med}$ deviate in the opposite directions
during such passages. We apply the estimator to the MW satellite system and
estimate the virial mass of the Milky Way of $M_{\rm 200c}=9.96\pm 1.45\times
10^{11}\,\,M_{\odot}$, in good agreement with several recent estimates using
other methods.

###### keywords:

galaxies: evolution, galaxies: Large Magellanic Cloud, galaxies: dwarf,
galaxies: halos

##  1 Introduction

Velocities of galaxies have been long used to infer information about the mass
of virialized objects, such as galaxy groups and clusters (Zwicky, 1933, 1937;
Smith, 1936), unvirialized outskirts of such objects (e.g., Diaferio et al.,
2005) and gravitationally interacting unvirialized pairs of galaxies (Kahn &
Woltjer, 1959; Li & White, 2008; González et al., 2014; Peñarrubia et al.,
2016).

Velocities of galaxies moving in and around the Milky Way (MW) galaxy have
also been used to probe its mass under the assumption that such galaxies have
negligible tangential velocity or are bound (Sandage, 1986; Zaritsky et al.,
1989; Kochanek, 1996). Boylan-Kolchin et al. (2013, cf. also ), for example,
used the velocity of motion of the Leo I galaxy in the outskirts of the MW
halo to constrain the virial mass of the Milky Way to $>10^{12}\,M_{\odot}$ at
the 95% confidence level. At the same time, modeling of the Sagittarius dwarf
stream indicated a Milky Way halo mass of $<10^{12}\,M_{\odot}$ (Gibbons et
al., 2014).

Such differences can be expected because the use of instantaneous velocity and
position of individual satellite galaxies to estimate host halo mass is prone
to biases if the satellite is in an atypical part of its orbit (Patel et al.,
2017a, b). Using the velocities of an ensemble of satellite galaxies can
mitigate some of the biases. Watkins et al. (2010) present an estimator of
halo mass and use radial velocities of the MW satellites to estimate MW halo
mass within 300 kpc distance, but find that the mass estimate is very
sensitive to the assumption about velocity anisotropy of the satellite system.
At the time, only a limited number of tangential velocity estimates existed
and thus estimates of the MW halo mass focused on the Bayesian methods
incorporating constraints from individual satellites using the phase-space
distribution function calibrated with cosmological simulations of MW-sized
halos (Patel et al., 2017a, b, 2018; Eadie et al., 2017).

The dramatic increase in the number of tangential velocity measurements for
many MW dwarf satellites and globular clusters with the advent of the Gaia
satellite data (McConnachie & Venn, 2020c, a; Li et al., 2021; Battaglia et
al., 2022) allowed statistical analyses of the velocity distributions using
much larger ensembles of objects (Callingham et al., 2019; Eadie & Jurić,
2019; Vasiliev, 2019; Posti & Helmi, 2019; Fritz et al., 2020; Li et al.,
2020; Rodriguez Wimberly et al., 2022). Similarly, accurate measurements of
phase-space information for stars in the Sagittarius dwarf spheroidal stream
and measurements of the LMC velocity allowed for the identification of a clear
misalignment between the stream track and the direction of the reflex-
corrected proper motions in the leading arm of the stream. Vasiliev et al.
(2021) showed that this misalignment is due to the interaction between the LMC
and MW, and its modeling allowed them to accurately estimate the total halo
masses of both objects.

The improved velocity accuracy and sophistication of the analyses led to
reduced uncertainties in halo mass estimates and reduced discrepancies between
different studies. Nevertheless, recent MW halo virial mass measurements still
allow a fairly wide range of values $\approx 6-20\times 10^{11}\,M_{\odot}$
(e.g., Vasiliev et al., 2021; Correa Magnus & Vasiliev, 2022; Koposov et al.,
2023; Zhou et al., 2023; Roche et al., 2024).

One systematic uncertainty in the inferences of halo mass from the positions
and velocities of dwarf satellite galaxies and globular clusters is the effect
of the Large Magellanic Cloud (LMC) on the orbits and velocities of such
tracers (e.g., Vasiliev et al., 2021; Pace et al., 2022, see Vasiliev 2023 for
a review). Effects of the LMC on the nearby stellar streams and dwarf galaxy
orbits indicate its halo mass of $1.3\pm 0.3\times 10^{11}\,M_{\odot}$ (Erkal
et al., 2019; Vasiliev et al., 2021; Shipp et al., 2021; Correa Magnus &
Vasiliev, 2022; Koposov et al., 2023; Watkins et al., 2024; Sheng et al.,
2024) or $\approx 15\%$ of the MW halo mass.

Such a large LMC mass results in perturbations of orbits of individual MW
satellites (Erkal et al., 2020; Pace et al., 2022; Makarov et al., 2023), a
significant reflex motion of the inner Milky Way, and asymmetries of the
potential in which dwarf satellites and orbit (Gómez et al., 2015; Erkal et
al., 2020; Garavito-Camargo et al., 2019; Garavito-Camargo et al., 2021;
D’Souza & Bell, 2022). These effects have been detected using different
tracers in the Milky Way halo (Petersen & Peñarrubia, 2020, 2021; Erkal et
al., 2021; Conroy et al., 2021). In particular, the reflex motion of the Milky
Way affects the velocities of the outer satellites, such as Antlia II (Ji et
al., 2021) and Leo I (Erkal et al., 2020). The use of the latter in particular
tended to produce higher halo mass estimates in prior analyses. Erkal et al.
(2020) also showed that the mass estimator of Watkins et al. (2010) applied to
the velocities of satellites and stellar tracers overestimates the enclosed
mass $M(<r)$ by up to $50\%$ at large $r$ if the LMC-induced reflex motion is
neglected. The presence of the bound Magellanic Clouds pair also affects
statistical estimates of the Milky Way halo mass (González et al., 2013).
Finally, the approach of the LMC has a large effect on the spatial
distribution of MW satellites (Nadler et al., 2020; Manwadkar & Kravtsov,
2022; Garavito-Camargo et al., 2023).

Large effects of the LMC on the phase-space distribution of the potential
Milky Way tracers need to be accounted for in theoretical inferences of halo
mass that use these tracers. This is challenging, however, because only
$\approx 5-10\%$ of the MW-sized halos have a pair of LMC-like satellite with
a pericenter close to $\approx 50$ kpc (Boylan-Kolchin et al., 2010; Busha et
al., 2011). Thus, large samples of halos need to be simulated with a high
resolution sufficient to resolve and model satellite evolution reliably.
Alternatively, the analyses can focus on the properties of tracers that are
not significantly affected by the dynamical perturbations due to the LMC, such
as specific angular momentum (Patel et al., 2018).

In this study, we explore the effects of the LMC on the velocity distribution
of satellite systems using a Caterpillar suite of high-resolution simulations
of MW-sized halos (Griffen et al., 2016, see Section 2). We demonstrate that
similarly to the real LMC, LMC-like objects in simulations induce a reflex
motion onto their host halo that creates a high-velocity tail in the
distribution of halocentric 3D satellite velocities (Section 3.1). We also
show how accounting for this tail affects the conclusion about MW halo mass
based on the velocity distribution. We then present a new halo mass estimator
that has only $\approx 8\%$ scatter around the true mass when applied to
simulated halos (Section 3.2) and which is affected by the LMC-like objects by
only $\lesssim 10\%$. We use the estimator calibrated with the simulated halos
to infer the MW halo mass of $M_{\rm 200c}=9.93\pm 1.45\times
10^{11}\,M_{\odot}$. We summarize our results and conclusions in Section 4.

##  2 Modeling Milky Way satellite system

We model the population of Milky Way dwarf satellite galaxies around the Milky
Way using tracks of haloes and subhaloes from the Caterpillar (Griffen et al.,
2016) suite of $N$-body simulations111https://www.caterpillarproject.org of 32
MW-sized haloes. We use the highest resolution suite LX14 to maximize the
dynamic range of halo masses probed by our modelling.

The haloes were identified using the modified version of the Rockstar halo
finder and the Consistent Trees Code (Behroozi et al., 2013), with
modification improving recovery of the subhaloes with high fraction of unbound
particles Griffen et al. (see discussion in Section 2.5 of 2016). As was shown
in Manwadkar & Kravtsov (2022, see their Fig. 1), the subhalo peak mass
function in the LX14 simulations is complete at $\mu=M_{\rm peak}/M_{\rm
host}\gtrsim 4\times 10^{-6}$, or $M_{\rm peak}\approx 4\times
10^{6}\,M_{\odot}$ for the host halo mass $M_{\rm host}\approx
10^{12}\,M_{\odot}$, even in the innermost regions of the host ($r<50$ kpc).
This is sufficient to model the full range of luminosities of observed Milky
Way satellites, as faintest ultrafaint dwarfs are hosted in haloes of $M_{\rm
peak}\gtrsim 10^{7}\,M_{\odot}$ in our model (Kravtsov & Manwadkar, 2022;
Manwadkar & Kravtsov, 2022).

The mass evolution tracks of subhaloes of MW-sized host haloes are used as
input for the GRUMPY galaxy formation model, which evolves various properties
of gas and stars of the galaxies they host. It is a regulator-type galaxy
formation framework (e.g., Krumholz & Dekel, 2012; Lilly et al., 2013;
Feldmann, 2013) designed to model galaxies of $\lesssim L_{\star}$ luminosity
(Kravtsov & Manwadkar, 2022) which follows the evolution of several key galaxy
properties by solving a system of coupled different equations. The model
accounts for UV heating after reionization and associated gas accretion
suppression onto small mass haloes, galactic outflows, a model for gaseous
disk and its size, molecular hydrogen mass, star formation, etc. The evolution
of the half-mass radius of the stellar distribution is also modeled. The
galaxy model parameters used in this study are identical to those used in
Manwadkar & Kravtsov (2022). Importantly for this study, the model was shown
to reproduce the luminosity function and radial distribution of the Milky Way
satellites and the size-luminosity relation of observed dwarf galaxies
(Manwadkar & Kravtsov, 2022). Here we use the model to predict luminosities
and stellar masses of satellite galaxies around the MW-sized haloes from the
Caterpillar suite.

Figure 1: The cumulative distribution (CDF) of distances of satellites to the
center of their parent halo for the Milky Way satellites (thick solid line).
The dark and light-shaded regions show the $68.2$ and $95.6$ percentiles of
the CDFs for the bootstrap resamples of the MW satellite distances at a given
distance $d$. The dot-dashed line shows the cumulative distance distribution
of subhalos of a peak mass of $M_{\rm peak}>10^{8}\,M_{\odot}$ in the main
host halos in the Caterpillar simulations. The light long-dashed line shows
the cumulative distribution of the distances for subhalos selected using the
distance matching procedure described in the text (Section 3.1). The figure
shows that the matching procedure selects subhalos with a distance
distribution statistically consistent with that of the Milky Way satellites.
Figure 2: The cumulative distribution of 3D velocities of satellites in the
Milky Way and simulated halos. The velocity distribution of the Milky Way
satellites is shown by the thick solid line. The dark and light-shaded regions
show the $68.2$ and $95.6$ percentiles of the CDFs for the bootstrap resamples
of the MW satellite velocities at a given velocity $v$. The dot-dashed line
shows the cumulative distance distribution of velocities for subhalos of the
peak mass of $M_{\rm peak}>10^{8}\,M_{\odot}$ in the main host halos in the
Caterpillar simulations. The light long-dashed line shows the cumulative
distribution of the velocities for subhalos selected using the distance
matching procedure described in the text (Section 3.1). The figure shows that
the matching procedure selects subhalos with a velocity distribution
statistically consistent with that of the Milky Way satellites. Figure 3: The
cumulative distribution of 3D velocities of satellites in the Milky Way and
simulated halos of different halo mass ranges. The velocity distribution of
the Milky Way satellites is shown by the thick solid line. The long-dashed,
dot-dashed, and dashed lines show velocity CDFs of satellites in the low,
intermediate, and high mass host halo mass ranges. The satellites have been
selected using the distance matching procedure described in the text (Section
3.1).  Figure 4: Top panel: evolution of distance (in physical kpc) of two
LMC-like satellites to the center of their host halo from $z=0.5$ to $z=0$ in
two different Caterpillar host halos. The satellites are on their first infall
and reach their pericenter of $\approx 40-50$ kpc at $z_{\rm LMC,peri}\approx
0.04-0.05$. Bottom panel: Evolution of halo masses (solid lines) and stellar
masses (dashed lines) of the two LMC-like satellites shown in the top panel.
The satellites reach peak halo masses of $M_{\rm peak}\approx 1.16\times
10^{11}\,M_{\odot}$ and $M_{\rm peak}\approx 1.35\times 10^{11}\,M_{\odot}$ at
$z\approx 0.1$ before they start to lose mass due to tidal stripping. The
stellar mass is computed using our galaxy formation model and continues to
grow due to ongoing star formation and reaches $M_{\star}\approx 3\times
10^{9}\,M_{\odot}$, comparable to the stellar mass of the observed LMC. Figure
5: Distribution of satellite subhalos (dots) in the plane of halocentric
distance and 3D velocity in one of the host halos in which an LMC-like
pericentric passage occurs at $z_{\rm LMC,peri}\approx 0.05$. The velocities
are taken with the sign of the radial velocity of each satellite. The dark and
light-shaded bands show the $15.9-84.1\%$ and $2.3-97.7\%$ ranges of subhalo
velocity distribution in distance bins. The large circle in each panel shows
the location of the LMC-like satellite. The gray dashed lines show escape
velocity profiles for the NFW halo with the same $M_{\rm 200c}$ and
concentration as the simulated halo at each redshift. The lower panel shows
the distribution at $z=z_{\rm LMC,peri}$ and includes the Milky Way dwarf
galaxy satellites (stars) for comparison. For reference, the real LMC is at
$r\approx 50$ kpc and a total velocity of $320$ km/s (e.g., Patel et al.,
2018) and has just passed its pericenter. The upper panel shows the velocity
distribution at $z=0$, which is approximately symmetric around zero, while at
$z_{\rm LMC,peri}$ it significantly shifted towards negative velocities.

Figure 6: Effect of the LMC on the cumulative distribution of 3D velocities
$v=|v_{\rm 3d}|$ (upper panel) and halocentric distances (lower panel) of
satellite subhalos. Dashed lines show the CDFs of satellites of the main host
at $z=0.1$ when the LMC-like satellite was still $\gtrsim 150$ kpc away from
the center, while solid lines show the CDFs of satellites at $z=0$ when the
LMC-like satellite is again at $d\approx 150$ kpc (see Figure 4 where this
LMC-like object is labeled as LMC1). The dot-dashed lines show the CDFs at
$z_{\rm LMC,peri}=0.05$ when the LMC-like satellite was close to the
pericenter. The velocity distributions are for satellites not selected to
match the distance distribution of the MW satellites. Note that the
distribution at $z_{\rm LMC,peri}$ has a significant tail of high-velocity
satellites, which is not present in the velocity CDFs at $z=0.1$ and $z=0$.
The LMC thus significantly skews the velocity distribution towards higher
velocities as it passes near the pericenter of its orbit. At the same time,
the distribution of distances is skewed towards smaller distances at the LMC
pericentric passage.

Figure 7: The distributions of the Kolmogorov-Smirnov test $p$-values for MW-
sized halos in the sample. The vertical lines show the 68 percentile of the
distribution of $p$-values obtained for bootstrap resamples of the satellites
in the simulated halos and the Milky Way satellites. The star symbols show the
medians of the distributions. Left panel shows the $p$-values computed using
velocity distributions not corrected for the effect of the LMC, while right
panel shows the case when simulated satellite velocities have been corrected
to account for the effect of the LMC. The correction significantly diminishes
the $p$-values for halos of $M_{\rm 200c}>1.3\times 10^{12}\,M_{\odot}$, but
they are still sufficiently high for most halos to be statistically consistent
with the Milky Way. In all cases, satellites in simulated halos are selected
to match the distance distribution of the MW satellites.

Figure 8: The distributions of the Mann-Whitney $U$ test $p$-values for MW-
sized halos in the sample. The vertical lines show the 68 percentile of the
distribution of $p$-values obtained for bootstrap resamples of the satellites
in the simulated halos and the Milky Way satellites. The star symbols show the
medians of the distributions. Left panel shows the case without taking into
account the effect of the LMC on the velocity distribution of satellites,
while right panel shows the case when simulated satellite velocities have been
corrected to account for the effect of the LMC. The correction significantly
diminishes the $p$-values for halos of $M_{\rm 200c}>1.3\times
10^{12}\,M_{\odot}$. As for the KS test $p$-values in Figure 7, the correction
significantly diminishes the $p$-values for halos of $M_{\rm 200c}>1.3\times
10^{12}\,M_{\odot}$, but they are still sufficiently high for most halos to be
statistically consistent with the Milky Way. In all cases, satellites in
simulated halos are selected to match the distance distribution of the MW
satellites.

Figure 9: Upper panel: effect of the LMC on the cumulative distribution
function of specific angular momentum of satellites
$j=|\mathbf{v}\times\mathbf{r}|$. The distributions are shown at $z=0.1$
(dashed line), when the LMC-like satellite was still $\gtrsim 150$ kpc away
from the center, and at $z=0$, when the LMC-like satellite is again at
$d\approx 150$ kpc (solid line). The dot-dashed line shows the CDF at $z_{\rm
LMC,peri}=0.05$ when the LMC-like satellite was close to the pericenter.
Unlike velocity distribution, the specific angular momentum distribution
changes little during pericentric passage of the LMC-sized halo. Lower panel:
specific angular momentum distribution of satellites in the MW-sized host
halos in the low- ($M_{\rm 200c}\leq 10^{12}\,M_{\odot}$), intermediate
($M_{\rm 200c}\in(1-1.2]\times 10^{12}\,M_{\odot}$), and high-mass bins
($M_{\rm 200c}>1.210^{12}\,M_{\odot}$). The dependence of the $j$ distribution
on halo mass is weaker than the dependence of the velocity distribution.
Figure 10: True halo mass $M_{\rm 200c}$ for the 32 host halos in the
Caterpillar sample vs the mass $M_{\rm 200c,est}$ estimated using the median
satellite distance within 300 kpc, the velocity dispersion of their 3D
velocity distribution, and the constant determining using the true $M_{\rm
200c}$ of the Caterpillar MW-sized host halos. The white long-dashed line is a
one-to-one relation, while dark and light-shaded bands show the $1\sigma$
($\approx 0.08$) and $2\sigma$ scatter of the sample estimated around this
relation. The two thin lines show the evolution of the halo mass and its
estimate from $z=0.15$ to $z=0$ for the two host halos in which the LMC-like
satellite passes through its pericenter at $z_{\rm LMC,peri}\approx
0.04-0.05$. Masses estimated for the MW-sized halos using eq. 1 in the sample
exhibit a tight linear relation with their true masses with only $\approx 8\%$
scatter.  Figure 11: Upper panel: effect of the LMC on the median halocentric
distance of satellites within 300 kpc of the host (dot-dashed line), variance
of the 3D velocity $\sigma^{2}$, and the halo mass estimate from $z=0.15$ to
$z=0$ in the halo, in which pericentric passage of the LMC-sized object on the
first infall occurs at $z\approx 0.04-0.05$ when distance to the host is
$\approx 50$ kpc (see Fig. 4). All quantities are normalized to their values
at $z=0$. Note that the LMC induces an upward fluctuation of $\sigma^{2}$ but
a downward fluctuation of $d_{\rm med}$. These fluctuations thus offset each
other in the estimate $M_{\rm 200c,est}$ using eq. 1, which exhibits a
fluctuation of only $\approx 10\%$. Lower panel: the ratio of the halo mass
estimated using equation 1 to its true mass during the LMC-sized object’s
pericentric passage. Prior to the pericentric passage, eq. 1 provided an
unbiased estimate of $M_{\rm 200c}$. During and immediately after the passage,
the estimate is biased by $\approx 5-10\%$.

##  3 Results

###  3.1 Effect of the LMC on the velocity distribution of satellites

For comparison of Caterpillar-simulated subhalos to Milky Way satellites, we
use the distance and velocity measurements found in Table 1 of Rodriguez
Wimberly et al. (2022, hereafter RW22). Galactocentric tangential velocities
and heliocentric radial velocities are taken from McConnachie & Venn (2020b)
and McConnachie & Venn (2020d), while the heliocentric distance and error
measurements are from Simon (2019), Karachentsev & Kashibadze (2006), Weisz et
al. (2016), and Torrealba et al. (2019). The measurements are converted to the
Galactocentric reference frame in RW22.

Comparison to the Milky Way requires a satellite population whose cumulative
Galactocentric distance distribution matches that of the MW satellites. To
obtain this, we draw a set of random distances from the MW satellite distance
distribution using the inverse transform sampling method and spline
approximation of the inverse cumulative distribution function of MW satellite
distances. We then selected subhalos in each simulated MW-sized halo that were
closest to each of the drawn distances. Figure 1 shows that this distance-
matching technique produces a satellite population with cumulative distance
distribution consistent with the distance distribution of MW the satellites.
At the same time, Figure 2 demonstrates that the 3D velocity distribution of
the distance-matched subhalos is statistically consistent with the velocity
distribution of Milky Way satellites.

Following RW22, we divided the simulated host halos into three mass bins: low
mass ($<10^{12}\ \,M_{\odot}$), intermediate mass ($1-1.2\times 10^{12}\
\,M_{\odot}$), and high mass ($>1.2\times 10^{12}\ \,M_{\odot}$). Figure 3
shows the cumulative velocity distribution of satellites in halos in each mass
bin and compares them to the velocity distribution of the Milky Way
satellites. The intermediate and high mass bins are closest to the
distribution of the Milky Way, while the low mass bin is shifted towards lower
velocities. Our results are somewhat different from the results shown in
Figure 3 of RW22 where the distribution of velocities in halos in the high
mass bin is the least consistent with the Milky Way. This is likely due to the
difference in the mass distribution of halos in the high-mass bin. The high-
mass halos in our sample contain masses in the range $1.2-1.6\times
10^{12}\,M_{\odot}$, while the sample of RW22 lacks halos in the
$1.2-1.4\times 10^{12}\ \,M_{\odot}$ range but contains halos of mass up to
$1.9\times 10^{12}\ \,M_{\odot}$. The velocity distribution of satellites in
their high-mass bin is thus shifted to higher velocities compared to our
sample.

Hosts 4 and 53 from the Caterpillar host halo tracks include the passage of
subhalos that reach peak masses of $M_{\rm peak}=1.16\times
10^{11}\,\,M_{\odot}$ and $M_{\rm peak}=1.35\times 10^{11}\,\,M_{\odot}$,
respectively. These are comparable to the mass of the LMC with $M_{\rm
peak}=1.3\pm 0.3\times 10^{11}\,M_{\odot}$ (Erkal et al., 2019; Vasiliev et
al., 2021; Shipp et al., 2021; Watkins et al., 2024). The time evolution of
the halocentric distance and mass of these subhalos are shown in Figure 4. The
satellites reach their peak masses $M_{\rm peak}\approx 1.16\times
10^{11}\,M_{\odot}$ and $M_{\rm peak}\approx 1.35\times 10^{11}\,M_{\odot}$ at
$z\approx 0.1$ and begin to lose mass due to tidal stripping during their
pericentric passage at $z\approx 0.05$. As the LMC-like satellites continue to
be tidally stripped, stellar mass steadily increases due to continued star
formation and reaches a peak of $M_{\star}\approx 3\times 10^{9}\,M_{\odot}$,
comparable to the observational estimates of the LMC’s stellar mass.

For one of the host halos, Figure 5 shows the satellite velocity distribution
at $z=0$ and at the time of pericentric passage $z\approx 0.05$, as well as
the 15.9 - 84.1 and 2.3 - 97.7 percentile intervals of the velocity
distribution in bins of distance. The effect of the LMC-like satellite is
greatest near the pericenter, where the subhalo velocity distribution becomes
skewed toward larger negative velocities. At $z=0$, the velocities return to a
symmetric distribution. The effect of the LMC-like subhalo passage on the
cumulative velocity distribution is shown in Figure 6, which shows that the
satellite distribution has a tail of higher velocities near the pericentric
passage. The cumulative velocity distribution at $z=0$, on the other hand, is
similar to the distribution before the LMC-like subhalo entered the host halo
at $z=0.1$. The lower panel of Figure 5 shows that the cumulative distance
distribution of satellites becomes more centrally concentrated during the
pericentric passage of the LMC-like subhalo. Like the velocity distribution,
the distance distribution at $z=0$ becomes closer to the distribution before
the passage at $z=0.1$.

The physical origin of the large change in the velocity and distance
distribution of satellites during the pericentric passage of an LMC-like
subhalo is the reflex motion of the central region of the MW-sized host halo
in response to the gravitational pull of the massive subhalo. During the
passage, the mean motion of the host halo relative to the subhalo system is
$\approx 40-50\,\rm km/s$, comparable to the measured magnitude of the Milky
Way disk reflex motion (e.g., Petersen & Peñarrubia, 2021).

The large effect of the LMC on the velocity and distance distribution needs to
be taken into account in the estimates of MW halo mass (Erkal et al., 2020).
Although the effect becomes small well after the passage, the LMC in the Milky
Way has just passed its pericenter and therefore effects on the MW satellite
velocities are expected to be significant. To illustrate the effect of change
in velocity distribution on halo mass inference, we approximate the difference
in the velocity distributions at the pericentric passage of the LMC-like
subhalo and at $z=0$ as a function of 3D velocity and apply the difference as
the correction to the velocity of each observed satellite of the Milky Way.

Figure 7 shows the Kolmogorov-Smirnov (KS) $p$-values for the consistency
between the velocity distributions of distance-matched satellite samples in
simulated host halos of different masses and the Milky Way satellites with and
without applying this correction. The figure shows the medians of the
$p$-values obtained for the bootstrap resamples of the velocities of the
distance-match subhalo and MW satellite samples, while vertical lines show the
extent of the 68th percentiles of the $p$-value distribution of these samples.
The figure shows that if the effect of the LMC is neglected, low halo masses
are disfavored by the KS test, while intermediate and high-mass halos are
consistent with the MW distribution. When the velocity correction is applied,
the $p$-values of the high-mass halos decrease markedly and the distribution
of $p$-values instead peaks in the intermediate halo masses.

The KS test is sensitive to the differences of both the mean shift and the
overall shapes of the distribution. RW22 used the Mann-Whitney $U$ test, which
is primarily sensitive to differences in the median of distributions. Figure 8
shows the results of the Mann-Whitney test with and without correction for the
LMC effect. The points and error bars show medians and 68th percentiles of the
$p$-values estimated for the bootstrap resamples of the velocities of MW
satellites and distance-matched satellite subhalos in the simulated host
halos. The figure shows results similar to the KS test results in Figure 7.
Without the correction, the $p$-values are high for halos with $M_{\rm
200c}\gtrsim 10^{12}\,M_{\odot}$, while smaller mass halos are disfavored.
When the LMC correction is applied to velocities, the $p$-values for halos
with $M_{\rm 200c}\gtrsim 1.3\times 10^{12}\,M_{\odot}$ decrease
substantially, although the overall range of $p$-values still does not exclude
these virial masses.

We note that when the distribution of the $p$-values for bootstrap resamples
is considered, the velocity distribution does not exclude any halos in the
mass range of the Caterpillar suite regardless of whether the LMC correction
to velocities is applied. The correction only lowers the $p$-values of halos
with $M_{\rm 200c}\gtrsim 1.3\times 10^{12}\,M_{\odot}$. Thus, the velocity
distribution is not a sensitive indicator of halo mass.

Patel et al. (2018) used the specific angular momenta of satellites to
estimate MW halo mass and showed that they are not sensitive to deviations
from dynamic equilibrium. The upper panel of Figure 9 shows the cumulative
distribution of specific angular momentum $j=|\mathbf{r}\times\mathbf{v}_{\rm
3D}|$ at three different epochs, $z=0.1$, $z=z_{\rm LMC,peri}$ and $z=0$, for
one of the host halos with an LMC-like satellite passing through its
pericenter at $z_{\rm LMC,peri}\approx 0.05$. The effect of the LMC-like
satellite on the $j$ distribution is indeed very small, which is consistent
with the conclusion of Patel et al. (2018). Thus, the use of the angular
momentum to estimate halo mass would not require a correction for the effect
of the LMC.

At the same time, the lower panel of Figure 9 shows that the cumulative $j$
distributions for halos in the three mass bins defined above differ
considerably less than the distributions of velocity. Specific angular
momentum distribution is thus not a very sensitive probe of halo mass. Indeed,
the halo mass is proportional to the second moment of the velocity
distribution $\langle v^{2}\rangle$, while Patel et al. (2018) find that the
halo mass depends on the median $j$ of satellites as $M_{\rm 200c}\propto
j^{1.3}$, and this relation has a large scatter. This motivated us to search
for an alternative proxy of $M_{\rm 200c}$ which is both sensitive to the
value of the mass and robust against effects of the LMC.

###  3.2 Robust halo mass estimator

The mass of a system in equilibrium should be proportional to a measure of the
size of the system and the velocity variance of some tracers of the halo
potential. We therefore searched for an estimator of this kind, taking into
account that during the LMC passage, the distance distribution of the
satellites shifts to smaller distances, while velocity distribution shifts to
larger velocities. We find that an accurate estimator of the halo mass
enclosing the density contrast of 200 relative to the critical mass, $M_{\rm
200c}$, is given by

$M_{\rm 200c,est}=c\,r_{\rm med}\,\sigma^{2}_{\rm 3D},$ (1)

where the constant $c=2.85\pm 0.028\times 10^{5}\,\,M_{\odot}$ estimated using
the true values of $M_{\rm 200c}$ of the 32 Caterpillar host halos, $r_{\rm
med}$ is the median halocentric distance of the satellites within $r_{\rm
max}=300$ kpc in kiloparsecs, and $\sigma^{2}_{\rm 3D}$ is the variance of
$v_{\rm 3D}\times\mathbf{v}_{r}/|v_{r}|$ for satellites within $r_{\rm max}$
in units of kilometers per second. Note that the variance of velocities with
the sign of radial velocities is used here, not just the variance of the
absolute values of $v_{\rm 3D}$. This variance measures the width of the 3D
velocity distribution of halo satellites shown in Figure 5.

Figure 10 demonstrates that the accuracy of this halo mass estimator is not
significantly impacted by the passage of the LMC. The ratio of the halo mass
obtained using the estimator to the halo mass given by Caterpillar for each
host shows only $s\approx 8\%$ scatter around the 1:1 ratio. The light lines
show the evolution of the estimated versus true mass for the two host halos
with an LMC-like satellite from before its infall at $z=0.15$ to $z=0$, well
after its pericentric passage. The estimator stays within the $2s$ range
through the duration of the satellites’ passages, indicating that the
estimator is not overly sensitive to the effects of an LMC-sized satellite.

The upper panel of Figure 11 shows the evolution of the factors of eq. 1
compared to their values at $z=0$ as a function of redshift. At $z=0.04-0.05$,
when the LMC-sized satellite is closest to the host, the median halocentric
distance of the satellites decreases while the variance of the 3D velocity
distribution increases. This reflects the changes in the velocity distribution
and distance cumulative distribution shown in Figures 5 and 6. As a result,
the value of the estimated mass deviates from its $z=0$ estimate by less than
10% during the pericentric passage of the LMC-like satellite at $z\approx
0.04$.

The lower panel of Figure 11 shows the evolution of the ratio of the estimated
halo mass to the true mass from $z=0.15$ to $z=0$. It shows that during and
after the pericentric passage, the estimated mass is $\approx 10\%$ higher
than the true mass. Thus, given that the scatter around the relation is
$\approx 8\%$, the LMC passage biases the mass by $1.25\sigma$. This means
that the mass estimator of eq. 1 is robust against effects of the LMC passage.
The small 10% correction due to the effect of the LMC can be applied to the
mass estimate of the Milky Way.

With the proportionality constant obtained from the simulated halos, using eq.
1 with all of the Milky Way satellites with measured 3D velocities and
applying 10% correction for the effect of the LMC yields a Milky Way halo mass
of $M_{\rm 200c}=9.96\pm 1.45\times 10^{11}\,M_{\odot}$, where uncertainty is
estimated using bootstrap resampling. If we use only the classical satellites
of the Milky Way ($M_{V}<-8$) we get a consistent estimate $M_{\rm 200c}=8\pm
2.4\times 10^{11}\,\,M_{\odot}$, albeit with a larger uncertainty.

##  4 Summary and conclusions

We presented a study illustrating the effects of the passage of LMC-mass
satellites on the distance and velocity distributions of satellites in the
Milky Way-sized halos. We also considered the effect of such passages on their
mass estimates. Our results and conclusions are as follows.

  * 1.

In agreement with results of Erkal et al. (2020), we find that the effect of
the passage of an LMC-like satellite on the velocity distribution is
substantial and this effect can bias mass estimates based on velocities. We
show that during the pericentric passage, the velocity distribution develops a
high-$v$ tail due to the reflex motion of the inner part of the halo, while
the satellite spatial distribution becomes more centrally concentrated.

  * 2.

We show that correcting for the effect of the LMC on the 3D velocity
distribution of satellites when using this distribution to estimate halo mass
favors halo masses in the range $1-1.3\times 10^{12}\,\,M_{\odot}$. However,
we also show that the velocity distribution is a rather low accuracy estimator
of halo mass as velocity distributions in halos of masses as low as $M_{\rm
200c}=8\times 10^{11}\,\,M_{\odot}$ and as high as $1.5\times
10^{12}\,\,M_{\odot}$ are still statistically consistent with the velocity
distribution of the MW satellites.

  * 3.

We considered the distribution of specific angular momenta $j$ as an estimator
of mass and confirmed that the $j$ distribution is not sensitive to the
passage of LMC-mass satellites. At the same time, we find that the $j$
distribution is a weak function of host halo mass and is therefore a
relatively poor mass indicator.

  * 4.

We present a new halo mass estimate that is proportional to the variance of 3D
velocities taken with the sign of the radial velocity of each satellite and
the median halocentric distance of the satellites (eq. 1). The velocity
variance is a measure of the width of the velocity distribution shown in
Figure 5. The estimator has only $8\%$ scatter around the mean $M_{\rm
200c,est}-M_{\rm 200c}$ relation.

  * 5.

We show that the velocity variance and median distance deviate in opposite
directions during the pericentric passage of the LMC-like satellites (Fig.
11). For this reason, the estimator deviates by $<2s$, where $s$ is the
scatter around the median $M_{\rm 200c,est}-M_{\rm 200c}$ relation during the
pericentric passage of the LMC-like subhalos.

  * 6.

We apply the estimator to the MW satellite system and estimate the virial mass
of the Milky Way to be $M_{\rm 200c}=9.96\pm 1.45\times 10^{11}\,\,M_{\odot}$,
in good agreement with several recent estimates using other methods.

## Acknowledgements

We are grateful to Vasily Belokurov and the UChicago structure formation group
for useful discussions during this project. We thank Alexander Ji and the
Caterpillar collaboration for providing halo tracks of the Caterpillar
simulations used in this study. AK was supported by the National Science
Foundation grants AST-1714658 and AST-1911111 and NASA ATP grant
80NSSC20K0512. SW was supported by the University of Chicago CCRF’s Quad
Research Scholarship program. Analyses presented in this paper were greatly
aided by the following free software packages: NumPy (Van Der Walt et al.,
2011), SciPy (Jones et al., 01 ), Matplotlib (Hunter, 2007), and GitHub. We
have also used the Astrophysics Data Service (ADS) and arXiv preprint
repository extensively during this project and the writing of the paper.

## Data Availability

Halo catalogs from the Caterpillar simulations are available at
https://www.caterpillarproject.org/. The GRUMPY model pipeline is available at
https://github.com/kibokov/GRUMPY. The data used in the plots within this
article are available on request to the corresponding author.

## References

  * Battaglia et al. (2022) Battaglia G., Taibi S., Thomas G. F., Fritz T. K., 2022, A&A, 657, A54
  * Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Wu H.-Y., Busha M. T., Klypin A. A., Primack J. R., 2013, ApJ, 763, 18
  * Boylan-Kolchin et al. (2010) Boylan-Kolchin M., Springel V., White S. D. M., Jenkins A., 2010, MNRAS, 406, 896
  * Boylan-Kolchin et al. (2013) Boylan-Kolchin M., Bullock J. S., Sohn S. T., Besla G., van der Marel R. P., 2013, ApJ, 768, 140
  * Busha et al. (2011) Busha M. T., Wechsler R. H., Behroozi P. S., Gerke B. F., Klypin A. A., Primack J. R., 2011, ApJ, 743, 117
  * Callingham et al. (2019) Callingham T. M., et al., 2019, MNRAS, 484, 5453
  * Conroy et al. (2021) Conroy C., Naidu R. P., Garavito-Camargo N., Besla G., Zaritsky D., Bonaca A., Johnson B. D., 2021, Nature, 592, 534
  * Correa Magnus & Vasiliev (2022) Correa Magnus L., Vasiliev E., 2022, MNRAS, 511, 2610
  * D’Souza & Bell (2022) D’Souza R., Bell E. F., 2022, MNRAS, 512, 739
  * Diaferio et al. (2005) Diaferio A., Geller M. J., Rines K. J., 2005, ApJ, 628, L97
  * Eadie & Jurić (2019) Eadie G., Jurić M., 2019, ApJ, 875, 159
  * Eadie et al. (2017) Eadie G. M., Springford A., Harris W. E., 2017, ApJ, 835, 167
  * Erkal et al. (2019) Erkal D., et al., 2019, MNRAS, 487, 2685
  * Erkal et al. (2020) Erkal D., Belokurov V. A., Parkin D. L., 2020, MNRAS, 498, 5574
  * Erkal et al. (2021) Erkal D., et al., 2021, MNRAS, 506, 2677
  * Feldmann (2013) Feldmann R., 2013, MNRAS, 433, 1910
  * Fritz et al. (2020) Fritz T. K., Di Cintio A., Battaglia G., Brook C., Taibi S., 2020, MNRAS, 494, 5178
  * Garavito-Camargo et al. (2019) Garavito-Camargo N., Besla G., Laporte C. F. P., Johnston K. V., Gómez F. A., Watkins L. L., 2019, ApJ, 884, 51
  * Garavito-Camargo et al. (2021) Garavito-Camargo N., Besla G., Laporte C. F. P., Price-Whelan A. M., Cunningham E. C., Johnston K. V., Weinberg M., Gómez F. A., 2021, ApJ, 919, 109
  * Garavito-Camargo et al. (2023) Garavito-Camargo N., et al., 2023, arXiv e-prints, p. arXiv:2311.11359
  * Gibbons et al. (2014) Gibbons S. L. J., Belokurov V., Evans N. W., 2014, MNRAS, 445, 3788
  * Gómez et al. (2015) Gómez F. A., Besla G., Carpintero D. D., Villalobos Á., O’Shea B. W., Bell E. F., 2015, ApJ, 802, 128
  * González et al. (2013) González R. E., Kravtsov A. V., Gnedin N. Y., 2013, ApJ, 770, 96
  * González et al. (2014) González R. E., Kravtsov A. V., Gnedin N. Y., 2014, ApJ, 793, 91
  * Griffen et al. (2016) Griffen B. F., Ji A. P., Dooley G. A., Gómez F. A., Vogelsberger M., O’Shea B. W., Frebel A., 2016, ApJ, 818, 10
  * Hunter (2007) Hunter J. D., 2007, Computing In Science & Engineering, 9, 90
  * Ji et al. (2021) Ji A. P., et al., 2021, ApJ, 921, 32
  * Jones et al. (01 ) Jones E., Oliphant T., Peterson P., et al., 2001--, SciPy: Open source scientific tools for Python, http://www.scipy.org/
  * Kahn & Woltjer (1959) Kahn F. D., Woltjer L., 1959, ApJ, 130, 705
  * Karachentsev & Kashibadze (2006) Karachentsev I. D., Kashibadze O. G., 2006, Astrophysics, 49, 3
  * Kochanek (1996) Kochanek C. S., 1996, ApJ, 457, 228
  * Koposov et al. (2023) Koposov S. E., et al., 2023, MNRAS, 521, 4936
  * Kravtsov & Manwadkar (2022) Kravtsov A., Manwadkar V., 2022, MNRAS, 514, 2667
  * Krumholz & Dekel (2012) Krumholz M. R., Dekel A., 2012, ApJ, 753, 16
  * Li & White (2008) Li Y.-S., White S. D. M., 2008, MNRAS, 384, 1459
  * Li et al. (2020) Li Z.-Z., Qian Y.-Z., Han J., Li T. S., Wang W., Jing Y. P., 2020, ApJ, 894, 10
  * Li et al. (2021) Li H., Hammer F., Babusiaux C., Pawlowski M. S., Yang Y., Arenou F., Du C., Wang J., 2021, ApJ, 916, 8
  * Lilly et al. (2013) Lilly S. J., Carollo C. M., Pipino A., Renzini A., Peng Y., 2013, ApJ, 772, 119
  * Makarov et al. (2023) Makarov D., Khoperskov S., Makarov D., Makarova L., Libeskind N., Salomon J.-B., 2023, MNRAS, 521, 3540
  * Manwadkar & Kravtsov (2022) Manwadkar V., Kravtsov A. V., 2022, MNRAS, 516, 3944
  * McConnachie & Venn (2020a) McConnachie A. W., Venn K. A., 2020a, Research Notes of the American Astronomical Society, 4, 229
  * McConnachie & Venn (2020b) McConnachie A. W., Venn K. A., 2020b, Research Notes of the AAS, 4, 229
  * McConnachie & Venn (2020c) McConnachie A. W., Venn K. A., 2020c, AJ, 160, 124
  * McConnachie & Venn (2020d) McConnachie A. W., Venn K. A., 2020d, The Astronomical Journal, 160, 124
  * Nadler et al. (2020) Nadler E. O., et al., 2020, ApJ, 893, 48
  * Pace et al. (2022) Pace A. B., Erkal D., Li T. S., 2022, ApJ, 940, 136
  * Patel et al. (2017a) Patel E., Besla G., Sohn S. T., 2017a, MNRAS, 464, 3825
  * Patel et al. (2017b) Patel E., Besla G., Mandel K., 2017b, MNRAS, 468, 3428
  * Patel et al. (2018) Patel E., Besla G., Mandel K., Sohn S. T., 2018, ApJ, 857, 78
  * Peñarrubia et al. (2016) Peñarrubia J., Gómez F. A., Besla G., Erkal D., Ma Y.-Z., 2016, MNRAS, 456, L54
  * Petersen & Peñarrubia (2020) Petersen M. S., Peñarrubia J., 2020, MNRAS, 494, L11
  * Petersen & Peñarrubia (2021) Petersen M. S., Peñarrubia J., 2021, Nature Astronomy, 5, 251
  * Posti & Helmi (2019) Posti L., Helmi A., 2019, A&A, 621, A56
  * Roche et al. (2024) Roche C., Necib L., Lin T., Ou X., Nguyen T., 2024, arXiv e-prints, p. arXiv:2402.00108
  * Rodriguez Wimberly et al. (2022) Rodriguez Wimberly M. K., et al., 2022, MNRAS, 513, 4968
  * Sandage (1986) Sandage A., 1986, ApJ, 307, 1
  * Sheng et al. (2024) Sheng Y., Ting Y.-S., Xue X.-X., Chang J., Tian H., 2024, arXiv e-prints, p. arXiv:2404.08975
  * Shipp et al. (2021) Shipp N., et al., 2021, ApJ, 923, 149
  * Simon (2019) Simon J. D., 2019, ARA&A, 57, 375
  * Smith (1936) Smith S., 1936, ApJ, 83, 23
  * Sohn et al. (2013) Sohn S. T., Besla G., van der Marel R. P., Boylan-Kolchin M., Majewski S. R., Bullock J. S., 2013, ApJ, 768, 139
  * Torrealba et al. (2019) Torrealba G., et al., 2019, MNRAS, 488, 2743
  * Van Der Walt et al. (2011) Van Der Walt S., Colbert S. C., Varoquaux G., 2011, ArXiv:1102.1523, p. 1
  * Vasiliev (2019) Vasiliev E., 2019, MNRAS, 484, 2832
  * Vasiliev (2023) Vasiliev E., 2023, Galaxies, 11, 59
  * Vasiliev et al. (2021) Vasiliev E., Belokurov V., Erkal D., 2021, MNRAS, 501, 2279
  * Watkins et al. (2010) Watkins L. L., Evans N. W., An J. H., 2010, MNRAS, 406, 264
  * Watkins et al. (2024) Watkins L. L., van der Marel R. P., Bennet P., 2024, ApJ, 963, 84
  * Weisz et al. (2016) Weisz D. R., et al., 2016, The Astrophysical Journal, 822, 32
  * Zaritsky et al. (1989) Zaritsky D., Olszewski E. W., Schommer R. A., Peterson R. C., Aaronson M., 1989, ApJ, 345, 759
  * Zhou et al. (2023) Zhou Y., Li X., Huang Y., Zhang H., 2023, ApJ, 946, 73
  * Zwicky (1933) Zwicky F., 1933, Helvetica Physica Acta, 6, 110
  * Zwicky (1937) Zwicky F., 1937, ApJ, 86, 217


11institutetext: Leiden University, Leiden, the Netherlands
11email: {j.shi, d.m.pelt<EMAIL_ADDRESS>

# SR4ZCT: Self-supervised Through-plane Resolution Enhancement for CT Images
with Arbitrary Resolution and Overlap

Jiayang Shi 11 Daniël M. Pelt 11 K. Joost Batenburg 1111

###### Abstract

Computed tomography (CT) is a widely used non-invasive medical imaging
technique for disease diagnosis. The diagnostic accuracy is often affected by
image resolution, which can be insufficient in practice. For medical CT
images, the through-plane resolution is often worse than the in-plane
resolution and there can be overlap between slices, causing difficulties in
diagnoses. Self-supervised methods for through-plane resolution enhancement,
which train on in-plane images and infer on through-plane images, have shown
promise for both CT and MRI imaging. However, existing self-supervised methods
either neglect overlap or can only handle specific cases with fixed
combinations of resolution and overlap. To address these limitations, we
propose a self-supervised method called SR4ZCT. It employs the same off-axis
training approach while being capable of handling arbitrary combinations of
resolution and overlap. Our method explicitly models the relationship between
resolutions and voxel spacings of different planes to accurately simulate
training images that match the original through-plane images. We highlight the
significance of accurate modeling in self-supervised off-axis training and
demonstrate the effectiveness of SR4ZCT using a real-world dataset.

###### Keywords:

CT Resolution enhancement Self-supervised learning

##  1 Introduction

CT is a valuable tool in disease diagnosis due to its ability to perform non-
invasive examinations [5]. The diagnostic accuracy of CT imaging is dependent
on the resolution of the images, which can sometimes be insufficient in
practice [8]. For medical CT images, the through-plane resolution is often
inferior to the in-plane axial resolution [17, 1], resulting in anisotropic
reconstructed voxels that can create difficulties in identifying lesions and
consequently lead to inaccurate diagnoses [3, 6, 1]. Furthermore, medical CT
scans are commonly performed in a helical trajectory, which can cause
overlapping slices along the through-plane axis [9, 2]. The overlap caused by
certain combinations of through-plane resolution and spacing introduces extra
blurriness in sagittal and coronal images and potentially leads to
inaccuracies in image interpretation [7, 15, 4].

To improve the through-plane resolution of CT images, various deep learning-
based methods have been proposed. Supervised methods that use pairs of low-
(LR) and high-resolution (HR) volumes have demonstrated promising results [14,
11, 19]. However, HR reference volumes may not always be available. As an
alternative, self-supervised methods have been introduced, training on in-
plane images and applying the learned model to through-plane images [18, 20].
These methods simulate LR training images from HR in-plane images that match
the through-plane images. However, they have limitations in handling CT images
with overlap or complicated resolution/overlap combinations. For example, [18]
can only handle CT images without overlap and integer resolution ratios, while
the use of convolutional operations in [20] restricts its usage to certain
specific combinations of overlap and resolution.

In this work, we propose a self-supervised method for enhancing the through-
plane resolution of CT images by explicitly modeling the relation between the
resolutions and spacings of different planes. Our method builds upon the same
idea of off-axis training introduced in [18, 20], but with the ability to
handle CT images with arbitrary resolution/overlap. We use accurate
interpolation techniques to simulate LR images from HR images, considering the
specific resolution and overlap parameters. Our contributions are: (1)
proposing a self-supervised method that can enhance resolution for CT images
with arbitrary resolution and overlap, (2) demonstrating the importance of
accurate modeling in off-axis training, and (3) applying our method to real-
world datasets with complicated resolution/overlap combinations, showcasing
its practical applicability.

##  2 Method

Problem Statement and Method. We present in this part a generalized problem
statement for enhancing the resolution of CT images and describe our proposed
method SR4ZCT as illustrated in Figure 2. We denote the reconstructed CT
volume by $\mathbf{I}(x,y,z)\in\mathbb{R}^{X\times Y\times Z}$, where $X,Y,Z$
represent the pixel numbers along $x,y,z$ axes, respectively. CT images can be
viewed along three different orientations: axial, coronal, and sagittal. We
denote the axial image, the coronal image, and the sagittal image at the
position $z^{\prime}$, $y^{\prime}$, and $x^{\prime}$ along the $z$ axis, the
$y$ axis, and the $x$ axis as
$\mathbf{a}_{z^{\prime}}=\mathbf{I}(:,:,z^{\prime})$,
$\mathbf{c}_{y^{\prime}}=\mathbf{I}(:,y^{\prime},:)$, and
$\mathbf{s}_{x^{\prime}}=\mathbf{I}(x^{\prime},:,:)$, respectively. The voxel
size of the CT volume is defined as $r^{x}\times r^{y}\times r^{z}$, where
often $r^{x}=r^{y}=r^{xy}<r^{z}$ for medical CT images. The spacing distance
$d$ between the centers of neighboring voxels along the $z$-axis can be
smaller than $r^{z}$, resulting in the overlap $o^{z}=r^{z}-d$ in the through-
plane direction. Given the target voxel size $r^{x}\times r^{y}\times
r^{z}_{tar}$ and overlap $o^{z}_{tar}$, the goal of resolution enhancement for
CT through-plane images is to reduce their voxel size and overlap in $z$
direction, hereby increasing the through-plane resolution.

Figure 1: Overview of our SR4ZCT method.

Figure 2: $\mathcal{F}^{ver}_{\downarrow}$ to simulate abitrary resolution and
overlap. The function linearly interpolates at various points (indicated by
different colors) inside the pixels, these interpolated points are then
averaged to be the values of the current pixels.

SR4ZCT is based on generating a training dataset that consists of HR axial
images as targets and their corresponding virtual through-plane-like axial
images as inputs. To achieve this, we define the downscaling function
$\mathcal{F}^{ver}_{\downarrow}$ as shown in Figure 2, which downscales an
axial image $\mathbf{a}_{i}$ with pixel sizes $r^{xy}\times r^{xy}$ by
employing linear interpolation in the vertical direction. The output is a
virtual through-plane-like image with pixel sizes $r^{z}\times r^{xy}$ and
overlaps $o^{z}$ between image rows. We similarly define the upscaling
function $\mathcal{F}^{ver}_{\uparrow}$ that upscales the virtual through-
plane-like image with pixel sizes $r^{z}\times r^{xy}$ and overlap $o^{z}$
using linear interpolation to produce a degraded axial image with pixel sizes
$r^{xy}\times r^{xy}$ and zero overlap. In addition, we define scaling
functions $\mathcal{F}^{hor}_{\downarrow}$ and $\mathcal{F}^{hor}_{\uparrow}$
that apply the scaling in the horizontal direction instead of the vertical
direction.

Training data is produced by first creating two network input images for each
axial image $\mathbf{a}_{i}$:
$\mathbf{a}_{i}^{ver}=\mathcal{F}^{ver}_{\uparrow}(\mathcal{F}^{ver}_{\downarrow}(\mathbf{a}_{i}))$
and
$\mathbf{a}_{i}^{hor}=\mathcal{F}^{hor}_{\uparrow}(\mathcal{F}^{hor}_{\downarrow}(\mathbf{a}_{i}))$.
Next, the simulated images are fed into the neural network $f_{\theta}$ to
learn the mapping from degraded images $\mathbf{a}_{i}^{ver}$ and
$\mathbf{a}_{i}^{hor}$ to their corresponding original axial image
$\mathbf{a}_{i}$. As the resolution enhancement is along the $z$ axis, which
is always in the vertical direction for coronal and sagittal images, the
horizontally degraded images and their corresponding axial images are rotated
$90\degree$ using rotation function $\mathcal{R}$. The training is performed
on all $Z$ axial images and their corresponding degraded images using the loss
function $L$, resulting in $2\times Z$ training pairs. The weights $\theta$ of
the neural network $f_{\theta}$ are determined by minimizing the total loss as

$\theta^{*}=\min_{\theta}\sum_{i=1}^{Z}L(f_{\theta}(\mathbf{a}_{i}^{ver}),\mathbf{a}_{i})+L(f_{\theta}(\mathcal{R}(\mathbf{a}_{i}^{hor})),\mathcal{R}(\mathbf{a}_{i})).$
(1)

After training, the network $f_{\theta}$ represents the learned mapping from
virtual through-plane-like axial images with resolution $r^{z}\times r^{xy}$
and overlap $o^{z}$ between image rows to HR axial image with $r^{xy}\times
r^{xy}$ resolution. With the assumption that CT images often share similar
features from different orientations, we can extend the learned mapping from
axial images to coronal and sagittal images. To enhance the resolution, we
first use $\mathcal{F}^{ver}_{\uparrow}$ to upscale coronal and sagittal
images with voxel sizes $r^{z}\times r^{xy}$ and overlap $o^{z}$ between image
rows to images with resolution $r^{xy}\times r^{xy}$. Subsequently, we apply
the trained neural network $f_{\theta}$ directly to the upscaled coronal and
sagittal images. The outputs of the neural network correspond to the improved
images with enhanced resolution.

In general, any image-to-image neural network could be used as $f_{\theta}$.
In this work, we use 2D MS-D networks [13] with 100 layers, trained for 200
epochs using L2 loss and ADAM optimizer [10]. The code is avaiable on GitHub.
111https://github.com/jiayangshi/SR4ZCT

##  3 Experiments and Results

Comparison with Supervised Learning. We conducted an experiment to evaluate
the effectiveness of SR4ZCT in improving the resolution of CT images,
comparing it to the supervised learning method. We selected four volumes (nr.
06, 14, 22, 23) from the Task06 Lung dataset of the Medical Segmentation
Decathlon [16]. These volumes had the same slice thickness ($0.625mm$) but
varying axial resolutions ranging from $0.74mm$ to $0.97mm$, resulting in a
total of 2170 axial images of size $512\times 512$ pixels. We downscaled the
axial images to a resolution of $1mm$, creating volumes with voxel sizes of
$1mm\times 1mm\times 0.65mm$. The through-plane resolution of the simulated
volumes ranged from $2.5mm$ with $1.25mm$ overlap to $6.25mm$ with $3.125mm$
overlap. We trained the supervised learning method using three LR and HR
volume pairs (numbers 06, 14, 22), and evaluated it on the remaining LR CT
volume (nr. 23). We utilized two distinct 2D MS-D networks to train the
coronal and sagittal images separately in a supervised manner, thereby
achieving the best possible results through supervised learning. We applied
SR4ZCT to the LR CT volume (nr. 23) and compared the results with those
obtained using the supervised learning method.

resolution | overlap | view | original | supervised | ours
---|---|---|---|---|---
PSNR | SSIM | PSNR | SSIM | PSNR | SSIM
$2.5mm$ | $1.25mm$ | cor | 37.71 | 0.972 | 51.94 | 0.998 | 46.48 | 0.995
sag | 39.12 | 0.977 | 52.50 | 0.998 | 47.68 | 0.995
$3.75mm$ | $1.875mm$ | cor | 33.65 | 0.934 | 42.96 | 0.987 | 42.36 | 0.985
sag | 35.08 | 0.945 | 43.73 | 0.989 | 43.39 | 0.987
$5mm$ | $2.5mm$ | cor | 31.43 | 0.895 | 37.70 | 0.965 | 38.46 | 0.967
sag | 32.87 | 0.913 | 39.36 | 0.973 | 39.55 | 0.973
$6.25mm$ | $3.125mm$ | cor | 29.91 | 0.9859 | 34.74 | 0.938 | 35.88 | 0.946
sag | 31.32 | 0.882 | 36.02 | 0.949 | 36.96 | 0.955

Table 1: Comparison of SR4ZCT and supervised learning on simulated CT volumes.
The used metrics are PSNR and SSIM calculated from coronal and sagittal images
in the central area, after cropping the empty regions.

Figure 3: The comparison of supervised learning and SR4ZCT on $6.25mm$
resolution and $3.125mm$ overlap dataset. Patches are shown for visualization.
The PSNR and SSIM of each patch are shown in the top- and bottom-left corners.

Table 1 shows that SR4ZCT successfully enhanced the resolution of all four
simulated CT volumes, as indicated by the improvements in both PSNR and SSIM.
While SR4ZCT showed slightly lower PSNR and SSIM values than the supervised
learning method for the $2.5mm$/$1.25mm$ case, its performance approached and
surpassed that of the supervised learning method as the resolution decreased.
Specifically, for the $3.75mm$/$1.875mm$ case, SR4ZCT achieved similar results
to supervised learning, whereas in the $5mm$/$2.5mm$ and $6.25mm$/$3.125mm$
cases, SR4ZCT outperformed supervised learning. Figure 3 presents an example
of the $6.25mm$/$3.125mm$ case, where both the supervised learning method and
SR4ZCT performed similarly in the lung area as SR4ZCT, but supervised learning
created artifacts on the edge area, highlighting one of its disadvantages as
it depends on the quality and quantity of training data. If the training data
slightly differs from the testing data or is insufficient in amount, it may
perform suboptimally. In contrast, SR4ZCT only makes use of the same data for
training and testing, reducing the potential gap between training and testing
data and yielding more robust results. Overall, our results demonstrate the
effectiveness of SR4ZCT in enhancing the resolution of through-plane images,
with similar or even superior performance compared with supervised learning.

Comparison with SMORE. SMORE [20] is a method designed for MRI images and
relies on the fact that MRI acquisition is akin to sampling in the Fourier
domain. It simulates LR images from HR images by removing the high-frequency
components. However, SMORE is not directly applicable to CT imaging.
Nonetheless, we implemented a version of SMORE based on its core idea of
simulating LR images by convolving HR images. While convolution can simulate
certain resolution/overlap combinations by adjusting the filter kernel and
stride, it fails to accurately simulate specific resolution/overlap
combinations, when the target resolution/overlap is not an integer multiple of
the HR resolution/overlap.

To compare our method with SMORE, we used the same volumes of the Medical
Segmentation Decathlon dataset [16] and downscaled them in the vertical
direction to simulate low through-plane resolution with overlaps of $2mm$,
$2.5mm$, and $3mm$ for a through-plane resolution of $5mm$. The results in
Table 2 demonstrate that SMORE performs better when the simulated LR images
using convolution match the through-plane resolution/overlap. For example,
applying 1D convolution with a filter width of 5 and a stride of 3 on the
$1mm\times 1mm$ HR axial images is equivalent to simulating a resolution of
$5mm$ and an overlap of $2mm$. It achieved higher PSNR values on the $5mm/2mm$
volumes than the one using a stride of 2. When the target overlap is $2.5mm$,
no convolution configuration accurately simulates it by adjusting the stride
and resulting in lower PSNRs. The same limitation applies when the ratio of
actual and target resolution is non-integer, making it impossible to find an
accurate convolution filter with a decimal filter width. These limitations of
convolution result in suboptimal performance. In contrast, our method
simulates LR images using accurate interpolation, which can handle arbitrary
resolution/overlap combinations.

|  $5mm$/$2mm$ |  $5mm$/$2.5mm$ |  $5mm$/$3mm$
---|---|---|---
SMORE stride 2 | $37.15$ | $38.39$ | $43.55$
SMORE stride 3 | $41.93$ | $41.86$ | $41.41$
ours | 42.93 | 43.30 | 44.75

Table 2: Results of SMORE and our method applied on the simulated dataset. The
used metric is PSNR, the values are computed based on the average PSNR of the
coronal images of all testing volumes. Output images are given in the
supplementary file.

Figure 4: Results of SR4ZCT applied on L291 from Low Dose CT Grand Challenge
[12]. We show two image patches of coronal and sagittal images. The first row
contains patches of original coronal and sagittal images, and the second row
shows the output of SR4ZCT. Extra results are shown in the supplementary file.

Real-world CT Images without Reference. We present this experiment designed to
evaluate the effectiveness of SR4ZCT on real-world CT images with anisotropic
resolution and overlap. We selected a CT volume of patient L291 from the Low
Dose CT Challenge [12], which has an in-plane resolution of $0.74mm$, through-
plane slice thickness of $3mm$ and $1mm$ overlap. As no reference volume
without overlap was available, we present the result only for visual
assessment.

Figure 4 illustrates the improvement in resolution of the original coronal and
sagittal images using SR4ZCT. The enhanced images exhibit sharper details
compared to the original images, as observed in the visual comparison. This
experiment provides evidence of SR4ZCT’s effectiveness in improving CT image
resolution in real-world scenarios where HR reference images are unavailable.

Correct Modeling is Essential. We present this experiment to demonstrate the
importance of correctly modeling the simulated training images from axial
images, which are used for training the neural network, to match the
resolution and overlap of coronal and sagittal images. To validate this
requirement, we intentionally introduced modeling errors into the training
images, simulating different combinations of resolution and overlap that
deviated by either $0.25mm$ or $0.5mm$ from the actual resolution of $5mm$ and
$2.5mm$ overlap.

Figure 5: The PSNR of cases where training images were inaccurately modeled.
Each block refers to a case where the corresponding resolution/overlap are
used to model the training images. The actual resolution/overlap are
$5mm$/$2.5mm$.

Figure 6: Results of cases where training images were not correctly modeled.
The actual through-plane resolution and overlap are $5mm$ and $2.5mm$. The
second and third column shows cases when modeled training images deviated by
$0.25mm$.

Fig 6 shows SR4ZCT achieved the highest PSNR when the training images
accurately reflected a resolution of $5mm$ and an overlap of $2.5mm$. Even
minor deviations, such as $0.25mm$ errors, resulted in a decrease in PSNR.
Interestingly, we observed that when the ratio of error in resolution and
overlap was similar, the performance was better than when the error was
present in only one parameter. Fig 6 provides examples where incorrect
modeling of resolution or overlap by $0.25mm$ led to artifacts, highlighting
the impact of inaccurate modeling.

Figure 7: Visual comparison of state-of-the-art CT images supervised super-
resolution methods SAINT [14], RPLHR [19], self-supervised method SMORE [20]
and our method SR4ZCT applied to a real-world sagittal image.

Fig 7 also demonstrates the importance of modeling the training images
accurately. We applied supervised methods SAINT [14], RPLHR [19], and self-
supervised SMORE [20], to Low Dose CT Grand Challenge [12]. We used the
provided pre-trained weights for SAINT and trained RPLHR as described on their
dataset. SAINT and RPLHR do not consider overlap in their method design. For
SMORE, we used convolution with a filter width of 4 and stride 3, which was
equivalent to simulating $2.96mm$ resolution and $0.74mm$ overlap images, the
closest possible combination to the actual $3mm$ resolution and $1mm$ overlap.
The blurriness caused by the overlap in the LR image was amplified in the
results of SAINT and RPLHR, while artifacts occurred in the output of SMORE.
In contrast, SR4ZCT improved the resolution and reduced blurriness by
accurately modeling the training images. This stresses the importance of
correctly modeling training images to ensure improvement for resolution
effectively.

##  4 Conclusion

In this work, we presented SR4ZCT, a self-supervised method designed to
improve the through-plane resolution of CT images with arbitrary resolution
and overlap. The method is based on the same assumption as [18, 20] that
images of a medical CT volume from different orientations often share similar
image features. By accurately simulating axial training images that match the
resolution and overlap of coronal and sagittal images, SR4ZCT effectively
enhances the through-plane resolution of CT images. Our experimental results
demonstrated that SR4ZCT outperformed existing methods for CT images with
complicated resolution and overlap combinations. It successfully enhanced the
resolution of real-world CT images without the need for reference images. We
also emphasized the crucial role of correctly modeling the simulated training
images for such off-axis training-based self-supervised method. In the future,
our method may have potential applications in a wide range of clinical
settings where differences in resolution across volume orientations currently
limit the 3D CT resolution.

## Acknowledgments

This research was financed by the European Union H2020-MSCA-ITN-2020 under
grant agreement no. 956172 (xCTing).

## References

  * [1] Angelopoulos, C., Scarfe, W.C., Farman, A.G.: A comparison of maxillofacial cbct and medical ct. Atlas of the oral and maxillofacial surgery clinics of North America 20(1), 1–17 (2012)
  * [2] Brink, J.A.: Technical aspects of helical (spiral) ct. Radiologic Clinics of North America 33(5), 825–841 (1995)
  * [3] Coward, J., Lawson, R., Kane, T., Elias, M., Howes, A., Birchall, J., Hogg, P.: Multi-centre analysis of incidental findings on low-resolution ct attenuation correction images. The British Journal of Radiology 87(1042), 20130701 (2014)
  * [4] Gavrielides, M.A., Zeng, R., Myers, K.J., Sahiner, B., Petrick, N.: Benefit of overlapping reconstruction for improving the quantitative assessment of ct lung nodule volume. Academic Radiology 20(2), 173–180 (2013)
  * [5] Hansen, P.C., Jørgensen, J., Lionheart, W.R.: Computed Tomography: Algorithms, Insight, and Just Enough Theory. SIAM (2021)
  * [6] He, L., Huang, Y., Ma, Z., Liang, C., Liang, C., Liu, Z.: Effects of contrast-enhancement, reconstruction slice thickness and convolution kernel on the diagnostic performance of radiomics signature in solitary pulmonary nodule. Scientific reports 6(1), 34921 (2016)
  * [7] Honda, O., Sumikawa, H., Johkoh, T., Tomiyama, N., Mihara, N., Inoue, A., Tsubamoto, M., Natsag, J., Hamada, S., Nakamura, H.: Computer-assisted lung nodule volumetry from multi-detector row ct: influence of image reconstruction parameters. European journal of radiology 62(1), 106–113 (2007)
  * [8] Iwano, S., Makino, N., Ikeda, M., Itoh, S., Tadokoro, M., Satake, H., Ishigaki, T.: Solitary pulmonary nodules: optimal slice thickness of high-resolution ct in differentiating malignant from benign. Clinical imaging 28(5), 322–328 (2004)
  * [9] Kasales, C., Hopper, K., Ariola, D., TenHave, T., Meilstrup, J., Mahraj, R., Van Hook, D., Westacott, S., Sefczek, R., Barr, J.: Reconstructed helical ct scans: improvement in z-axis resolution compared with overlapped and nonoverlapped conventional ct scans. AJR. American journal of roentgenology 164(5), 1281–1284 (1995)
  * [10] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. In: Bengio, Y., LeCun, Y. (eds.) 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings (2015), http://arxiv.org/abs/1412.6980
  * [11] Liu, Q., Zhou, Z., Liu, F., Fang, X., Yu, Y., Wang, Y.: Multi-stream progressive up-sampling network for dense ct image reconstruction. In: Medical Image Computing and Computer Assisted Intervention–MICCAI 2020: 23rd International Conference, Lima, Peru, October 4–8, 2020, Proceedings, Part VI 23. pp. 518–528. Springer (2020)
  * [12] McCollough, C.H., Bartley, A.C., Carter, R.E., Chen, B., Drees, T.A., Edwards, P., Holmes III, D.R., Huang, A.E., Khan, F., Leng, S., et al.: Low-dose ct for the detection and classification of metastatic liver lesions: results of the 2016 low dose ct grand challenge. Medical physics 44(10), e339–e352 (2017)
  * [13] Pelt, D.M., Sethian, J.A.: A mixed-scale dense convolutional neural network for image analysis. Proceedings of the National Academy of Sciences 115(2), 254–259 (2018)
  * [14] Peng, C., Lin, W.A., Liao, H., Chellappa, R., Zhou, S.K.: Saint: spatially aware interpolation network for medical slice synthesis. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 7750–7759 (2020)
  * [15] Ravenel, J.G., Leue, W.M., Nietert, P.J., Miller, J.V., Taylor, K.K., Silvestri, G.A.: Pulmonary nodule volume: effects of reconstruction parameters on automated measurements—a phantom study. Radiology 247(2), 400–408 (2008)
  * [16] Simpson, A.L., Antonelli, M., Bakas, S., Bilello, M., Farahani, K., Van Ginneken, B., Kopp-Schneider, A., Landman, B.A., Litjens, G., Menze, B., et al.: A large annotated medical image dataset for the development and evaluation of segmentation algorithms. arXiv preprint arXiv:1902.09063 (2019)
  * [17] Tsukagoshi, S., Ota, T., Fujii, M., Kazama, M., Okumura, M., Johkoh, T.: Improvement of spatial resolution in the longitudinal direction for isotropic imaging in helical ct. Physics in Medicine & Biology 52(3), 791 (2007)
  * [18] Xie, H., Lei, Y., Wang, T., Tian, Z., Roper, J., Bradley, J.D., Curran, W.J., Tang, X., Liu, T., Yang, X.: High through-plane resolution ct imaging with self-supervised deep learning. Physics in Medicine & Biology 66(14), 145013 (2021)
  * [19] Yu, P., Zhang, H., Kang, H., Tang, W., Arnold, C.W., Zhang, R.: Rplhr-ct dataset and transformer baseline for volumetric super-resolution from ct scans. In: Medical Image Computing and Computer Assisted Intervention–MICCAI 2022: 25th International Conference, Singapore, September 18–22, 2022, Proceedings, Part VI. pp. 344–353. Springer (2022)
  * [20] Zhao, C., Dewey, B.E., Pham, D.L., Calabresi, P.A., Reich, D.S., Prince, J.L.: Smore: a self-supervised anti-aliasing and super-resolution algorithm for mri using deep learning. IEEE transactions on medical imaging 40(3), 805–817 (2020)


The ALPACA experiment: The project of the first sub-PeV gamma-ray observation
in the southern sky

Teruyoshi Kawashima1$\star$ for The ALPACA Collaboration

1 Institute for Cosmic Ray Research, the University of Tokyo, 5-1-5,
Kashiwanoha, Kashiwa-shi, Chiba, 277-8582, Japan

*<EMAIL_ADDRESS>
21st International Symposium on Very High Energy Cosmic Ray Interactions
(ISVHE- CRI 2022)
Online, 23-28 May 2022 10.21468/SciPostPhysProc.?

## Abstract

The ALPACA experiment is a project aiming to observe sub-PeV gamma rays for
the first time in the southern hemisphere. The main goal of ALPACA is to
identify PeVatrons, the accelerators of Galactic PeV cosmic rays, by observing
sub-PeV pion-decay gamma rays generated in interactions between PeV cosmic
rays and the interstellar medium. This new air shower experiment is located at
an altitude of 4,740 m above sea level in the middle of Mt. Chakartaya in
Bolivia. The air shower array consists of 401 scintillation counters covering
an 83,000 m2 surface area. In addition, a water-Cherenkov-type muon detector
array with an area of 3,700 m2 is installed to discriminate gamma rays from
background cosmic rays. The prototype array ALPAQUITA will start data taking
in 2022 and will extend to ALPACA in 2024. We report on a general introduction
to ALPACA, the progress of the project, and the sensitivity to sub-PeV gamma
rays.

###### Contents

  1. 1 Introduction
  2. 2 Outline of the ALPACA project
    1. 2.1 ALPACA
    2. 2.2 The prototype array: ALPAQUITA
  3. 3 Construction status of ALPAQUITA
  4. 4 Conclusion

##  1 Introduction

The identity of the PeV-cosmic-ray accelerators in our galaxy, called
PeVatrons, is unknown, and identifying them is an important challenge in
astrophysics. When PeV cosmic rays interact with the interstellar medium,
neutral pions are produced, which immediately decay into two photons above 100
TeV (sub-PeV gamma rays). Therefore, observations of sub-PeV gamma rays are
essential to the search for PeVatrons and have been actively performed in
recent years. The Tibet AS$\gamma$ collaboration opened the door to sub-PeV
gamma-ray astronomy for the first time with the detection of sub-PeV gamma
rays from the Crab Nebula[1] followed by HAWC[2] and LHAASO[3]. Then, sub-PeV
gamma-ray observations from several sources were reported by these
collaborators, such as SNR G106.3+2.7, Cygnus Cocoon, and J1908+0621[4, 5, 6].
In addition, the detection of diffuse sub-PeV gamma rays by the Tibet
AS$\gamma$ collaboration has established the existence of PeVatrons in our
galaxy[7]. Thus, the sub-PeV gamma-ray astronomy in the Northern Hemisphere
has been greatly developed by these collaborators. However, the direct
identification of PeVatrons is still to be achieved.

As reported in the H.E.S.S. Galactic plane survey[8], more high-energy gamma-
ray sources exists in the southern sky including some promising PeVatron
candidates. The H.E.S.S. group reported gamma-ray emissions of up to 40 TeV
(spectral index: $-2.3$) from the Galactic Center region, strongly pointing
out the possibility of the Galactic Center being a PeVatron[9]. Furthermore,
the H.E.S.S. group reported the existence of an unidentified object HESS
J1702-420A, a PeVatron candidate having a hard gamma-ray spectrum (spectral
index: $-1.5$) extending above 100 TeV without cutoff[10]. However, in the
southern sky, sub-PeV gamma-ray observations are still almost untouched, so it
is important to start surveys.

The Andes Large area PArticle detector for Cosmic ray physics and Astronomy
(ALPACA) collaboration aims to demonstrate the first sub-PeV gamma ray
observations in the southern hemisphere using an air-shoer-array detector. The
array is going to be construct in Bolivia and plans to start observations in
2024. The construction of ALPAQUITA, a prototype array of ALPACA, is currently
underway and is scheduled to start observations in 2022. This paper presents a
general introduction to the ALPACA project, the sensitivites of ALPACA and
ALPAQUITA to celestial gamma rays, and the construction status of ALPAQUITA.

##  2 Outline of the ALPACA project

The ALPACA project is an international collaboration between Bolivia, Mexico,
and Japan. A new array based on Tibet AS$\gamma$ design[11] is constructed at
an altitude of 4,740 m above sea level in the middle of Mt. Chacaltaya,
Bolivia (16∘23’ S, 68∘08’ W and the atmospheric depth of
$572.4\,\mathrm{g/cm^{2}}$).

###  2.1 ALPACA

ALPACA array consists of two main components; surface air shower detectors
(SDs) and underground muon detectors (MDs) (see Figure 1). The array has a
geometrical area of $\sim 83,000\,\mathrm{m}^{2}$ and captures air showers
produced by the interactions of primary gamma rays and cosmic rays with the
atmosphere. Each SD consists of a $1\,\mathrm{m}\times 1\,\mathrm{m}\times
5\,\mathrm{cm}$ plastic scintillator and a 2-inch PMT covered by a steel box.
401 SDs are installed at 15 m intervals and can determine the energy and
direction of arrival of a primary particle. On the other hand, MDs comprises
water Cherenkov detectors constructed 2 m underground and its total area is
$\sim 3,700\,\mathrm{m}^{2}$. Each MD consists of 16 cells, each of which have
a area of $56\,\mathrm{m}^{2}$ and contains a water layer with a 1.5 m depth.
Most of the electromagnetic components in air showers are absorbed by the soil
before reaching MDs, but only muons with energies above 1 GeV penetrate into
the water layer of the MDs and produce Cherenkov radiations, which are
collected by a 20-inch PMT suspended downward at the cell’s ceiling after
several reflection on the wall and floor. Counting the number of muons in air
showers allow us to distinguish gamma-ray-induced air showers from (hadronic)
cosmic-ray-induced ones.

Figure 1: Schematic view of the ALPACA array. The small black and large grey
squares mean SDs and MDs, respectively. The gray rectangle in the center
indicates the central electronics hut. The solid octagon in the lower left
corner shows the location and area of the prototype array ALPAQUITA (see
Section 2.2).

ALPACA’s high duty cycle ($\sim 100\,\%$) and wide field of view ($\sim 2$
steradian) make it suitable for detecting sub-PeV gamma rays with low arrival
frequencies. For 100 TeV gamma rays, the angular and energy resolutions are
estimated at about 0.2 degrees and 25 %, respectively. The hadron rejection
capability of MDs is more than 99.9 % at 100 TeV, while keeping about 90 % of
the gamma-ray efficiency.

Figure 2 shows the sensitivity curve of ALPACA to gamma rays. The horizontal
and vertical axis indicates the energies of gamma rays and the integral gamma-
ray flux multiplied by energy. Blue points show the diffuse gamma-ray spectra
in the Galactic Center region measured by H.E.S.S. telescope. The blue dashed
and solid line indicates a power-law function (spectral index: -2.3) with an
exponential cutoff and the hadronic emission model, respectively, fitted to
the results of H.E.S.S., showing that ALPACA touches the Galactic Center flux
above 100 TeV in one year.

###  2.2 The prototype array: ALPAQUITA

The prototype-array ALPAQUITA consists of 97 SDs and one MD as shown in Figure
1. The ALPAQUITA construction is ongoing (see Section 3) and observations are
scheduled to begin in 2022. Figure 3 shows the sensitivity curve of ALPAQUITA
together with that of ALPACA sensitivity and the energy spectra of gamma-ray
sources that resides in the ALPACA field of view. The thick black curve
indicates that ALPAQTUIA can detect several gamma-ray sources including HESS
J1702-420A in one year.

Figure 2: Sensitivity curve of ALPACA (thick red curve) overlaid on the
integral flux of the diffuse gamma rays from the Galactic Center region. The
curve is estimated by MC simulation based on Tibet AS+MD within a zenith angle
of 60 degrees. Black points indicates the Crab Nebula spectrum measured by
Tibet AS$\gamma$[1] and the dashed line is a power-law function with an
exponential cutoff fitted to the spectrum. Figure 3: Sensitivity curves of
ALPAQUITA and ALPACA overlaid on the gamma-ray energy spectra of gamma-ray
sources in the ALPACA field of view[12]. For each energy spectrum, the solid
line represents the observed region by HESS[8][10] and HAWC[13] and the dashed
line shows the extrapolated region.

##  3 Construction status of ALPAQUITA

The central electronics hut has already been built and 20 SDs are already
located around the hut (see Figure 4). Installations of the perimeters, cables
trenches, lightining rods, and internet Connection are also completed. ALPACA
powerline is supplied by the branch of the Chacaltaya Observatory power line.
We have restarted the construction of ALPAQUITA in June 2022 as the COVID-19
situation is improved.

Figure 4: Photo of ALPACA (ALPAQUITA) site. Mt. Chakartaya can be seen in the
background. The building with the red roof is the central electronics hut. The
multiple white cubes around the hut are SDs.

##  4 Conclusion

Sub-PeV gamma-ray observaton is crutial to search PeV-cosmic-ray accelerators
called PeVatrons. Tibet AS$\gamma$, HAWC, LHAASO Collaboration has pioneered
sub-PeV energy gamma-ray astronomy in the Northern Hemisphere. The ALPACA
experiment is a project aiming to search for sub-PeV gamma-ray sources in the
Southern Hemisphere for the first time using a new air shower array in
Bolivia. ALPACA is scheduled to start its observatio in 2024 and is expected
to yield interesting results after one year of observations. ALPAQUITA, a
prototype array under construction, can detect several bright sources in one-
year observation such as HESS J1702-420A. The infrastructure for ALPAQUITA is
being built locally and its DAQ is going to begin in 2022. Ultimately, the
ALPACA project will identify Pevatrons in the southern sky and will help
understand their acceleration mechanism of PeV comis rays.

## Acknowledgements

The ALPACA project is supported by the Japan Society for the Promotion of
Science (JSPS) through Grants-in-Aid for Scientific Research (A) 19H00678,
Scientific Research (B) 19H01922, Scientific Research (B) 20H01920, Scientific
Research (S) 20H05640, Scientific Research (B) 20H01234, Scientific Research
(C) 22K03660, and Specially Promoted Research 22H04912, the LeoAtrox
supercomputer located at the facilities of the Centro de Análisis de Datos
(CADS), CGSAIT, Universidad de Guadalajara, México, and by the joint research
program of the Institute for Cosmic Ray Research (ICRR), The University of
Tokyo. K. Kawata is supported by the Toray Science Foundation. E. de la Fuente
thanks Coordinación General Académica y de Innovación (CGAI-UDG), cuerpo
académico PRODEP-UDG-CA-499, Carlos Iván Moreno, Cynthia Ruano, Rosario
Cedano, and Diana Naylleli, for financial and administrative support during
sabbatical year stay at the ICRR on 2021. I. Toledano-Juarez acknowledges
support from CONACyT, México; grant 754851.

## References

  * [1] M. Amenomori et al., First Detection of Photons with Energy beyond 100 TeV from an Astrophysical Source, Phys. Rev. Lett., 123, 051101 (2019), 10.1103/PhysRevLett.123.051101.
  * [2] A. U. Abeysekara et al., Measurement of the Crab Nebula Spectrum Past 100 TeV with HAWC, Astrophys. J., 881:134 (2019), 10.3847/1538-4357/ab2f7d.
  * [3] F. Aharonian et al., Observation of the Crab Nebula with LHAASO-KM2A $-$ a performance study, Chinese Phys. C, 45, 025002 (2021), 10.1088/1674-1137/abd01b.
  * [4] The Tibet AS$\gamma$ Collaboration, Potential PeVatron supernova remnant G106.3+2.7 seen in the highest-energy gamma rays, Nature Astronomy 5, 460-464 (2021), 10.1038/s41550-020-01294-9.
  * [5] A. U. Abeysekara et al., HAWC observations of the acceleration of very-high-energy cosmic rays in the Cygnus Cocoon, Nature Astronomy 5, 465-471 (2021), 10.1038/s41550-021-01318-y.
  * [6] Zhen Cao et al., Ultrahigh-energy photons up to 1.4 petaelectronvolts from 12 $\gamma$-ray Galactic sources, Nature, 594, 33-36 (2021), 10.1038/s41586-021-03498-z.
  * [7] M. Amenomori et al., First Detection of sub-PeV Diffuse Gamma Rays from the Galactic Disk: Evidence for Ubiquitous Galactic Cosmic Rays beyond PeV Energies, Phys. Rev. Lett., 126, 141101 (2021), 10.1103/PhysRevLett.126.141101.
  * [8] H.E.S.S. Collaboration, The H.E.S.S. Galactic plane survey, A&A 612, A1 (2018), 10.1051/0004-6361/201732098.
  * [9] H.E.S.S. Collaboration, Acceleration of petaelectronvolt protons in the Galactic Centre, Nature, 531, 476-479 (2016), 10.1038/nature17147.
  * [10] H. Abdalla et al., Evidence of 100 TeV $\gamma$-ray emission from HESS J1702-420: A new PeVatron candidate, A&A 653, A152 (2021), 10.1051/0004-6361/202140962.
  * [11] T. K. Sako et al., Exploration of a 100 TeV gamma-ray northern sky using the Tibet air-shower array combined with an underground water-Cherenkov muon-detector array Author links open overlay panel, Astroparticle Physics, 32, 177-184 (2009), 10.1016/j.astropartphys.2009.07.006.
  * [12] S.Kato et al., Detectability of southern gamma-ray sources beyond 100 TeV with ALPAQUITA, the prototype experiment of ALPACA, Experimental Astronomy 52, 85-107 (2021), 10.1007/s10686-021-09796-8.
  * [13] A. U. Abeysekara et al., Multiple Galactic Sources with Emission Above 56 TeV Detected by HAWC, Phys. Rev. Lett. 124, 021102 (2020), 10.1103/PhysRevLett.124.021102.