# Line-of-sight Elongation and Hydrostatic Mass Bias of the Frontier Fields
Galaxy Cluster Abell 370

Keiichi Umetsu Academia Sinica Institute of Astronomy and Astrophysics
(ASIAA), No. 1, Section 4, Roosevelt Road, Taipei 10617, Taiwan  Shutaro Ueda
Academia Sinica Institute of Astronomy and Astrophysics (ASIAA), No. 1,
Section 4, Roosevelt Road, Taipei 10617, Taiwan  Bau-Ching Hsieh Academia
Sinica Institute of Astronomy and Astrophysics (ASIAA), No. 1, Section 4,
Roosevelt Road, Taipei 10617, Taiwan  Mario Nonino INAF-Osservatorio
Astronomico di Trieste, Via G. B. Tiepolo 11, 34143, Trieste, Italy  I-Non
Chiu Tsung-Dao Lee Institute, and Key Laboratory for Particle Physics,
Astrophysics and Cosmology, Ministry of Education, Shanghai Jiao Tong
University, Shanghai 200240, China  Department of Astronomy, School of Physics
and Astronomy, and Shanghai Key Laboratory for Particle Physics and Cosmology,
Shanghai Jiao Tong University, Shanghai 200240, China  Academia Sinica
Institute of Astronomy and Astrophysics (ASIAA), No. 1, Section 4, Roosevelt
Road, Taipei 10617, Taiwan  Masamune Oguri Center for Frontier Science, Chiba
University, 1-33 Yayoi-cho, Inage-ku, Chiba 263-8522, Japan  Department of
Physics, Graduate School of Science, Chiba University, 1-33 Yayoi-Cho, Inage-
Ku, Chiba 263-8522, Japan  Sandor M. Molnar Academia Sinica Institute of
Astronomy and Astrophysics (ASIAA), No. 1, Section 4, Roosevelt Road, Taipei
10617, Taiwan  Anton M. Koekemoer Space Telescope Science Institute, 3700 San
Martin Drive, Baltimore, MD 21218, USA  Sut-Ieng Tam Academia Sinica Institute
of Astronomy and Astrophysics (ASIAA), No. 1, Section 4, Roosevelt Road,
Taipei 10617, Taiwan

###### Abstract

We present a detailed weak-lensing and X-ray study of the Frontier Fields
galaxy cluster Abell 370, one of the most massive known lenses on the sky,
using wide-field $BR_{\mathrm{C}}z^{\prime}$ Subaru/Sprime-Cam and Chandra
X-ray observations. By combining two-dimensional (2D) shear and azimuthally
averaged magnification constraints derived from Subaru data, we perform a
lensing mass reconstruction in a free-form manner, which allows us to
determine both radial structure and 2D morphology of the cluster mass
distribution. In a triaxial framework assuming a Navarro–Frenk–White density
profile, we constrain the intrinsic structure and geometry of the cluster halo
by forward modeling the reconstructed mass map. We obtain a halo mass
$M_{200}=(1.54\pm 0.29)\times 10^{15}\,h^{-1}\,M_{\odot}$, a halo
concentration $c_{200}=5.27\pm 1.28$, and a minor–major axis ratio
$q_{a}=0.62\pm 0.23$ with uninformative priors. Using a prior on the line-of-
sight alignment of the halo major axis derived from binary merger simulations
constrained by multi-probe observations, we find that the data favor a more
prolate geometry with lower mass and lower concentration. From triaxial lens
modeling with the line-of-sight prior, we find a spherically enclosed gas mass
fraction of $f_{\mathrm{gas}}=(8.4\pm 1.0)\%$ at
$0.7\,h^{-1}\,\mathrm{Mpc}\sim 0.7r_{500}$. When compared to the hydrostatic
mass estimate ($M_{\mathrm{HE}}$) from Chandra observations, our triaxial
weak-lensing analysis yields spherically enclosed mass ratios of $1-b\equiv
M_{\mathrm{HE}}/M_{\mathrm{WL}}=0.56\pm 0.09$ and $0.51\pm 0.09$ at
$0.7\,h^{-1}\,\mathrm{Mpc}$ with and without using the line-of-sight prior,
respectively. Since the cluster is in a highly disturbed dynamical state, this
represents the likely maximum level of hydrostatic bias in galaxy clusters.

cosmology: observations — dark matter — gravitational lensing: weak — X-rays:
galaxies: clusters — galaxies: clusters: individual (A370)

††software: IMCAT package (Kaiser et al., 1995), xspec (v12.11.1; Arnaud,
1996), Sherpa (Freeman et al., 2001; Doe et al., 2007; Burke et al., 2021),
CIAO (v4.13; Fruscione et al., 2006), PYATOMDB (Foster & Heuer, 2020), SCAMP
software (Bertin, 2006), SWARP (Bertin et al., 2002), EMCEE (Foreman-Mackey et
al., 2013), LEPHARE (Ilbert et al., 2006)

##  1 Introduction

Galaxy clusters can provide a range of valuable information from the physics
driving structure formation to the nature of dark matter and dark energy.
Their matter content reflects that of the universe: $\sim 85\%$ dark matter
and $\sim 15\%$ baryons, with $\sim 90\%$ of the baryons residing in the hot
intracluster medium (ICM). Determining the evolution of the abundance of rare
massive clusters provides powerful cosmological constraints, especially on the
matter density parameter, $\Omega_{\mathrm{m}}$, and the amplitude of linear
density fluctuations, $\sigma_{8}$ (e.g., see Mantz et al., 2015, and
references therein). Conversely, an accurate determination of the total mass
of galaxy clusters using direct mass probes, such as weak gravitational
lensing, is essential to harness the full potential of cluster cosmology
(e.g., Pratt et al., 2019; Chiu et al., 2021; Tam et al., 2022).

In the context of the standard $\Lambda$ cold dark matter ($\Lambda$CDM)
model, galaxy clusters are non-spherical in shape and better approximated as
triaxial halos (Jing & Suto, 2002), with a preference for prolateness over
oblateness and preferentially aligned with surrounding filaments (Bett et al.,
2007). Cluster-scale halos can be characterized as triaxial ellipsoids with a
typical minor-to-major axis ratio of $0.4$–$0.5$ (Bonamigo et al., 2015),
where more massive objects tend to be more prolate. Thus, while the intrinsic
shape and orientation of galaxy clusters contain unique cosmological
information (Okumura & Taruya, 2020), they can also introduce significant
scatter and bias in cluster mass estimates due to the unknown orientation of
cluster halos. In particular, gravitational lensing is sensitive to such
projection effects (e.g., Becker & Kravtsov, 2011).

According to cosmological $N$-body simulations, “superlens” clusters
characterized by large Einstein radii ($\theta_{\mathrm{Ein}}\geqslant
30\arcsec$ for a source redshift of $z_{s}=2$) represent the most lensing-
biased population of clusters, with their major axis preferentially aligned
with the observer’s line of sight (Hennawi et al., 2007; Oguri & Blandford,
2009; Meneghetti et al., 2010a, 2011). A statistical bias in favor of prolate
structure pointed close to the observer arises, because such a halo geometry
can boost the projected mass density and hence the lensing signal. In
particular, major mergers of two clusters colliding nearly along the line of
sight provide a possible mechanism for producing a powerful superlens (e.g.,
see the case of Cl0024+1654; Umetsu et al., 2010).

Abell 370 (hereafter A370; a.k.a. PSZ2 G172.98$-$53.55) at $z=0.375$ is known
as a prominent strong lens with an Einstein radius of
$\theta_{\mathrm{Ein}}\approx 34\arcsec$ (for $z_{s}=2$; see Table 1) and is
the first galaxy cluster where gravitational lensing has been observed in the
form of a giant luminous arc (Soucail et al., 1987, 1988). A370 is also among
the most massive clusters based on weak gravitational lensing, with an
estimated virial mass of $M_{\mathrm{vir}}\sim 2\times
10^{15}\,h^{-1}\,M_{\odot}$ (Umetsu et al., 2011; Hoekstra et al., 2015, all
relevant symbols are defined at the end of this section). Because of its large
projected mass and exceptional lensing strength, A370 was selected as one of
the six Hubble Frontier Fields (Lotz et al., 2017) and has recently been
targeted by the Beyond Ultra-deep Frontier Fields and Legacy Observations
(BUFFALO; Steinhardt et al., 2020) with the Hubble Space Telescope (HST),
which expands the area coverage of the Frontier Fields in optical and near-
infrared pass bands.

Lensing studies of A370 reveal a bimodal mass distribution in the core
elongated in the north–south direction (e.g., Kneib et al., 1993; Umetsu et
al., 1999; Richard et al., 2010; Diego et al., 2018; Lagattuta et al., 2017,
2019; Ghosh et al., 2021). Strait et al. (2018) combined strong and weak
lensing constraints from Frontier Fields imaging and spectroscopic
observations to reconstruct the central mass distribution of A370. Their mass
map shows two dominant peaks associated with the two brightest cluster
galaxies (BCGs), with the northern peak much less concentrated than the
southern one and slightly offset from the stellar mass distribution. These
bimodal and offset features are often an indication of recent major merger
activity (e.g., Bradač et al., 2008).

In contrast to its extreme mass and exceptional lensing properties, A370 is
intriguingly faint in both X-ray and Sunyaev–Zel’dovich effect (SZE) signals
and does not follow the X-ray/SZE observable–mass scaling relations (see
Czakon et al., 2015). The X-ray brightness distribution revealed from Chandra
observations is highly elongated in the north–south direction, showing a
disturbed morphology with the brightest X-ray peak located about halfway
between the two BCGs (Molnar et al., 2020). The irregular morphology in X-ray
emission with a large elongation similar to that of the mass distribution is a
strong indication that the cluster is far from hydrostatic equilibrium (Lee &
Suto, 2003).

Recently, $N$-body hydrodynamical simulations of binary cluster mergers
constrained by lensing, X-ray, SZE, and optical spectroscopic observations
suggest that A370 is a massive post-major merger viewed after the second core
passage in the infalling phase, just before the third core passage (Molnar et
al., 2020). In this post-collision phase, the gas has not settled into the
gravitational potential of the cluster, which explains the faintness of the
X-ray and SZE signals. These results also suggest that the mass distribution
of A370 is highly elongated along the current direction of the collision axis,
which is oriented close to the line of sight in their best-matching
simulation.

In this paper, we present a detailed weak-lensing and X-ray study of A370
using wide-field $BR_{\mathrm{C}}z^{\prime}$ imaging taken with Suprime-Cam on
the Subaru telescope and high-quality data from the Chandra X-ray Observatory.
The primary aims of this paper are to obtain an accurate inference of the
three-dimensional (3D) mass model of A370 from a full triaxial analysis of
two-dimensional (2D) weak-lensing data and to determine the level of
hydrostatic mass bias and the gas mass fraction as a function of cluster
radius. The key for this study is to perform a lensing mass reconstruction in
an unbiased manner, from which to constrain both radial structure and 2D
morphology of the cluster mass distribution. To this end, we perform an
improved joint shear and magnification analysis of 2D Subaru weak-lensing
data, revisiting our earlier one-dimensional (1D) work presented in Umetsu et
al. (2011). Since A370 is extremely massive and in a highly disturbed
dynamical state, this analysis will provide a constraint on the likely maximum
level of the hydrostatic bias expected in galaxy clusters.

This paper is organized as follows. Section 2 describes the basic theory of
cluster weak lensing and outlines the methodology used to reconstruct the
cluster mass distribution. Section 3 describes details of the Subaru
observations, reduction procedures, and weak-lensing analysis. Section 4
presents the results of our mass reconstruction, followed by our triaxial
modeling in Section 5. Section 6 describes the X-ray data analysis. Section 7
compares the weak-lensing and Chandra mass profiles. Finally, a summary is
given in Section 8.

Throughout this paper, we assume a spatially flat $\Lambda$CDM cosmology with
$\Omega_{\mathrm{m}}=0.3$, $\Omega_{\Lambda}=0.7$, and a Hubble constant of
$H_{0}=100~{}h$ km s-1 Mpc-1 with $h=0.7$. In this cosmology, $1\arcmin$
corresponds to $216.8\,h^{-1}\,\mathrm{kpc}$ at the cluster redshift of
$z=0.375$. The reference center of the cluster is taken to be the optical
cluster center defined by Lotz et al. (2017, see also ):
$\mathrm{R.A.}=\mathrm{02:39:52.9}$, $\mathrm{decl.}=-\mathrm{01:34:36.5}$
(see Table 1).

We denote the critical density of the universe at a particular redshift $z$ as
$\rho_{\mathrm{c}}(z)=3H^{2}(z)/(8\pi G)$, with $H(z)$ the Hubble function. We
generally denote spherical and projected radii from the cluster center as $r$
and $r_{\perp}$, respectively, and reserve the symbol $R$ for ellipsoidal
cluster radii. We adopt the standard notation $M_{\Delta}$ (or
$M_{\Delta\mathrm{m}}$) to denote the mass enclosed within a sphere of radius
$r_{\Delta}$ (or $r_{\Delta\mathrm{m}}$) within which the mean overdensity
equals $\Delta$ (or $\Delta_{\mathrm{m}}$) times $\rho_{\mathrm{c}}(z)$ (or
the mean background density $\rho_{\mathrm{m}}(z)$). We compute the virial
mass and radius, $M_{\mathrm{vir}}$ and $r_{\mathrm{vir}}$, using an
expression for $\Delta_{\mathrm{vir}}$ based on the spherical collapse model
(Bryan & Norman, 1998). For its ellipsoidal counterpart $R_{\Delta}$, see
Section 5.2. We use “$\log$” to denote the base-10 logarithm and “$\ln$” to
denote the natural logarithm. All quoted errors are at the $1\sigma$
confidence level unless otherwise stated. The AB magnitude system is used
throughout.

Table 1: Properties of the galaxy cluster A370

Parameter | Value
---|---
ID | A370
Reference center position (J2000.0) |
R.A. | 02:39:52.9
Decl. |  $-$01:34:36.5
X-ray emission centroid (J2000.0) |
R.A. | 02:39:53.2
Decl. |  $-$01:34:35.1
Redshift | $0.375$
Velocity dispersion (km s-1) | $1520\pm 93$
X-ray temperature (keV) | $8.77^{+0.33}_{-0.34}$
Einstein radius ($\arcsec$) |  $33.9\pm 1.1$ for $z_{s}=2$

Note. — The optical cluster center is at the midpoint of the two BCGs (Lotz et
al., 2017; Steinhardt et al., 2020). Units of right ascension are hours,
minutes, and seconds, and units of declination are degrees, arcminutes, and
arcseconds. The cluster velocity dispersion is derived from spectroscopic
observations of Lagattuta et al. (2022) in the core region of the cluster (see
also Lagattuta et al., 2019). The X-ray emission centroid is determined from a
2D $\beta$-model fit to Chandra X-ray observations (see Section 6). The
average temperature of the cluster is measured from the Chandra X-ray spectrum
in the radial range $\in[50,500]\,h^{-1}\,\mathrm{kpc}$ centered on the X-ray
centroid. The Einstein radius is constrained by detailed strong lens modeling
by Kawamata et al. (2018).

##  2 Weak-lensing Methodology

###  2.1 Basics of Galaxy–Cluster Weak Lensing

The effects of weak gravitational lensing on background galaxies are
characterized by the convergence, $\kappa$, and the shear with spin 2
rotational symmetry, $\gamma=|\gamma|e^{2i\phi_{\gamma}}$ (for reviews, see
Bartelmann & Schneider, 2001; Umetsu, 2020). In this work, we closely follow
the notation of Umetsu (2020).

The lensing convergence $\kappa$ alone causes an isotropic magnification of
galaxy images and it is defined as the surface mass density $\Sigma$ of a lens
in units of the critical surface density for gravitational lensing,
$\kappa=\Sigma/\Sigma_{\mathrm{cr}}$, where

$\displaystyle\Sigma_{\mathrm{cr}}(z_{l},z_{s})$
$\displaystyle=\frac{c^{2}}{4\pi GD_{l}(z_{l})}\frac{1}{\beta(z_{l},z_{s})},$
(1) $\displaystyle\beta(z_{l},z_{s})$
$\displaystyle=\begin{cases}D_{ls}(z_{l},z_{s})/D_{s}(z_{s})&\mathrm{for}~{}z_{s}>z_{l}\\\
0&\mathrm{for}~{}z_{s}\leqslant z_{l}\end{cases},$

with $c$ the speed of light, $G$ the gravitational constant, and
$D_{l}(z_{l})$, $D_{s}(z_{s})$, and $D_{ls}(z_{l},z_{s})$ the observer–lens,
observer–source, and lens–source angular diameter distances, respectively. The
dimensionless factor $\beta(z_{l},z_{s})$ describes the geometric lensing
efficiency as a function of lens redshift $z_{l}$ and source redshift $z_{s}$.
The shear and convergence thus depend on $(z_{l},z_{s})$ as well as on the
image position $\theta$.

The gravitational shear field $\gamma(\mbox{\boldmath$\theta$})$ is directly
observable from image ellipticities of background galaxies in the weak-lensing
regime, $|\kappa|\ll 1,|\gamma|\ll 1$. The shear and convergence fields are
related by

$\gamma(\mbox{\boldmath$\theta$})=\int\\!d^{2}\theta^{\prime}\,D(\mbox{\boldmath$\theta$}-\mbox{\boldmath$\theta$}^{\prime})\kappa(\mbox{\boldmath$\theta$}^{\prime})$
(2)

with $D(\mbox{\boldmath$\theta$})$ the complex kernel
$D(\mbox{\boldmath$\theta$})=(\theta_{2}^{2}-\theta_{1}^{2}-2i\theta_{1}\theta_{2})/(\pi|\mbox{\boldmath$\theta$}|^{4})$.
The key observable for weak shear lensing in the subcritical regime is the
complex reduced shear,

$g(\mbox{\boldmath$\theta$})=\frac{\gamma(\mbox{\boldmath$\theta$})}{1-\kappa(\mbox{\boldmath$\theta$})},$
(3)

which remains invariant under the global transformation
$\kappa(\mbox{\boldmath$\theta$})\to\lambda\kappa(\mbox{\boldmath$\theta$})+1-\lambda$
and
$\gamma(\mbox{\boldmath$\theta$})\to\lambda\gamma(\mbox{\boldmath$\theta$})$
with an arbitrary constant $\lambda\neq 0$ (for a fixed source redshift
$z_{s}$). This is referred to as the mass-sheet degeneracy (Schneider & Seitz,
1995). This degeneracy can be broken or alleviated, for example, by measuring
the magnification factor $\mu$ in the subcritical regime,

$\mu(\mbox{\boldmath$\theta$})=\frac{1}{[1-\kappa(\mbox{\boldmath$\theta$})]^{2}-|\gamma(\mbox{\boldmath$\theta$})|^{2}}\equiv\frac{1}{\Delta_{\mu}(\mbox{\boldmath$\theta$})}.$
(4)

We note that in practical applications to magnification bias measurements,
this degeneracy can be lifted only if the unlensed mean source background
density is known or can be estimated from the data (see Section 2.5). The
magnification factor $\mu$ transforms as
$\mu(\mbox{\boldmath$\theta$})\to\lambda^{-2}\mu(\mbox{\boldmath$\theta$})$.
For simplicity of notation, we often use the inverse magnification
$\Delta_{\mu}=\mu^{-1}$.

The reduced shear $g_{1,2}$ can be decomposed into the tangential component
$g_{+}=\gamma_{+}/(1-\kappa)$ and the $45^{\circ}$-rotated cross-shear
component $g_{\times}=\gamma_{\times}/(1-\kappa)$ with respect to a given
reference point. The tangential shear $\gamma_{+}(\theta)$ averaged around a
circle of projected radius $\theta$ is related to the excess surface mass
density $\Delta\Sigma(\theta)$ through the following identity:

$\displaystyle\Sigma_{\mathrm{cr}}\gamma_{+}(\theta)$
$\displaystyle=\Sigma(<\theta)-\Sigma(\theta)\equiv\Delta\Sigma(\theta),$ (5)

where $\Sigma(\theta)$ is the azimuthally averaged surface mass density at
radius $\theta$ and $\Sigma(<\theta)$ is the average surface mass density
interior to $\theta$. The azimuthally averaged cross-shear
$\gamma_{\times}(\theta)$ is expected to vanish if the signal is due to weak
lensing.

###  2.2 Source Redshift Distribution

We consider a population of source galaxies characterized by their mean
(unlensed) redshift distribution, $\overline{N}(z)$. In general, we use
different magnitude, color, size, and quality cuts in background selection for
measuring the shear and magnification effects. This results in different
$\overline{N}(z)$ for shear and magnification. The source-averaged mean
lensing depth $\langle\beta^{n}\rangle_{X}$ ($n=1,2,\dots$) for a given
population ($X=g,\mu$) is

$\langle\beta^{n}\rangle_{X}=\left[\int_{0}^{\infty}\\!dz\,\overline{N}_{X}(z)\beta^{n}(z_{l},z)\right]\left[\int_{0}^{\infty}\\!dz\,\overline{N}_{X}(z)\right]^{-1}.$
(6)

In general, $\overline{N}(z)$ for a given lens can include foreground
galaxies. The contribution from unlensed objects with $\beta=0$ is thus taken
into account in the calculation of $\langle\beta^{n}\rangle_{X}$.

We introduce the relative lensing strength of a given source population
$\langle W\rangle_{X}=\langle\beta\rangle_{X}/\beta_{\infty}$ with
$\beta_{\infty}\equiv\beta(z_{l},z_{s,\infty})$ defined relative to a
reference source in the far background at redshift $z_{s,\infty}$ (Bartelmann
& Schneider, 2001). We use a reference redshift of $z_{s,\infty}=20000$, which
was adopted in the CLASH program (Umetsu et al., 2014; Merten et al., 2015).
The associated critical surface density is
$\Sigma_{\mathrm{cr},\infty}(z_{l})=c^{2}/(4\pi GD_{l})\beta_{\infty}^{-1}$.
Hereafter, we use the far-background fields
$\kappa_{\infty}(\mbox{\boldmath$\theta$})$ and
$\gamma_{\infty}(\mbox{\boldmath$\theta$})$ to describe the projected mass
distribution of the cluster.

###  2.3 Pixelized Mass Distribution

We pixelize the convergence field,
$\kappa_{\infty}(\mbox{\boldmath$\theta$})=\Sigma_{\mathrm{cr},\infty}^{-1}\Sigma(\mbox{\boldmath$\theta$})$,
into a regular grid of pixels and describe
$\kappa_{\infty}(\mbox{\boldmath$\theta$})$ by a linear combination of basis
functions $B(\mbox{\boldmath$\theta$}-\mbox{\boldmath$\theta$}^{\prime})$ as

$\kappa_{\infty}(\mbox{\boldmath$\theta$})=\Sigma_{\mathrm{cr},\infty}^{-1}\sum_{n=1}^{N_{\mathrm{pix}}}B(\mbox{\boldmath$\theta$}-\mbox{\boldmath$\theta$}_{n})\,\Sigma_{n}.$
(7)

To avoid the loss of information due to oversmoothing, we choose the basis
function to be the Dirac delta function,
$B(\mbox{\boldmath$\theta$}-\mbox{\boldmath$\theta$}^{\prime})=(\Delta\theta)^{2}\delta^{2}_{\mathrm{D}}(\mbox{\boldmath$\theta$}-\mbox{\boldmath$\theta$}^{\prime})$,
with $\Delta\theta$ a constant grid spacing. The 2D cluster lensing signal is
specified by a vector of model parameters containing cell-averaged surface
mass densities (Umetsu et al., 2015, 2018),

$\mbox{\boldmath$m$}=\left\\{\Sigma_{n}\right\\}_{n=1}^{N_{\mathrm{pix}}}$ (8)

with $\Sigma_{n}=\Sigma(\mbox{\boldmath$\theta$}_{n})$
($n=1,2,\dots,N_{\mathrm{pix}}$). The complex shear field is then expressed as

$\gamma_{\infty}(\mbox{\boldmath$\theta$})=\Sigma_{\mathrm{cr},\infty}^{-1}\sum_{n=1}^{N_{\mathrm{pix}}}(D\otimes
B)(\mbox{\boldmath$\theta$}-\mbox{\boldmath$\theta$}_{n})\,\Sigma_{n}$ (9)

with $D\otimes B=\pi^{-1}(\Delta\theta)^{2}D$ an effective complex kernel
(Equation (2)). Hence, both $\kappa_{\infty}(\mbox{\boldmath$\theta$})$ and
$\gamma_{\infty}(\mbox{\boldmath$\theta$})$ can be expressed as linear
combinations of mass coefficients.111Because of the choice of the basis
function, an unbiased extraction of the mass coefficients
$\\{\Sigma_{n}\\}_{n=1}^{N_{\mathrm{pix}}}$ is possible by performing a
spatial integral of $\kappa_{\infty}(\mbox{\boldmath$\theta$})$ over a certain
area. Such operations include spatial smoothing, azimuthal averaging for the
radial profile extraction, and fitting with smooth parametric functions.

###  2.4 Reduced Shear Field

We use the reduced shear field as the primary constraint from our weak-lensing
observations. The source-averaged reduced shear
$g_{n}=g(\mbox{\boldmath$\theta$}_{n})$ is measured from shape measurements of
background galaxies onto a regular grid of $N_{\mathrm{pix}}$ pixels
($n=1,2,\dots,N_{\mathrm{pix}}$) as

$g_{n}=\left[\displaystyle\sum_{k}S(\mbox{\boldmath$\theta$}_{(k)},\mbox{\boldmath$\theta$}_{n})w_{g(k)}g_{(k)}\right]\left[\displaystyle\sum_{k}S(\mbox{\boldmath$\theta$}_{(k)},\mbox{\boldmath$\theta$}_{n})w_{g(k)}\right]^{-1}$
(10)

where $S(\mbox{\boldmath$\theta$},\mbox{\boldmath$\theta$}^{\prime})$ is a
spatial window function, $g_{(k)}$ is an estimate of
$g(\mbox{\boldmath$\theta$})$ for the $k$th galaxy at
$\mbox{\boldmath$\theta$}_{(k)}$, and $w_{g(k)}$ is its statistical weight,
$w_{g(k)}=1/(\sigma^{2}_{g(k)}+\alpha^{2}_{g})$, with $\sigma^{2}_{g(k)}$ the
error variance of $g_{(k)}$. The $\alpha_{g}$ parameter is set to a typical
value of the shear dispersion $\sigma_{g}=0.4$ found in Subaru weak-lensing
observations (e.g., Umetsu et al., 2009, 2014).

The source-averaged expectation (denoted by a hat symbol) for the observable
$g_{n}$ (Equation (10)) is given by (Seitz & Schneider, 1997; Umetsu et al.,
2015)

$\widehat{g}(\mbox{\boldmath$\theta$}_{n})\simeq\frac{\langle
W\rangle_{g}\gamma_{\infty}(\mbox{\boldmath$\theta$}_{n})}{1-f_{W,g}\langle
W\rangle_{g}\kappa_{\infty}(\mbox{\boldmath$\theta$}_{n})},$ (11)

where $\langle W\rangle_{g}$ is the source-averaged relative lensing strength
(see Section 2.2) and $f_{W,g}\equiv\langle W^{2}\rangle_{g}/\langle
W\rangle_{g}^{2}=\langle\beta^{2}\rangle_{g}/\langle\beta\rangle_{g}^{2}$ is a
dimensionless correction factor of the order unity. The error variance
$\sigma_{g,n}^{2}$ for $g_{n}$ is expressed as

$\displaystyle\sigma^{2}_{g,n}=\frac{\sum_{k}S^{2}(\mbox{\boldmath$\theta$}_{(k)},\mbox{\boldmath$\theta$}_{n})w_{g(k)}^{2}\sigma^{2}_{g(k)}}{\left[{\sum_{k}S(\mbox{\boldmath$\theta$}_{(k)},\mbox{\boldmath$\theta$}_{n})w_{g(k)}}\right]^{2}}.$
(12)

We adopt the top-hat window of radius $\theta_{\mathrm{f}}$ (Merten et al.,
2009; Umetsu et al., 2015, 2018),
$S(\mbox{\boldmath$\theta$},\mbox{\boldmath$\theta$}^{\prime})=H(\theta_{\mathrm{f}}-|\mbox{\boldmath$\theta$}-\mbox{\boldmath$\theta$}^{\prime}|)$,
with $H(x)$ the Heaviside function defined such that $H(x)=1$ if $x\geqslant
0$ and $H(x)=0$ otherwise. The shape-noise covariance matrix for
$g_{\alpha,n}=g_{\alpha}(\mbox{\boldmath$\theta$}_{n})$ is then given as
(Oguri et al., 2010)

$\left(C_{g}\right)_{\alpha\beta,mn}=\frac{1}{2}\delta_{\alpha\beta}\sigma_{g,m}\sigma_{g,n}\xi_{H}(|\mbox{\boldmath$\theta$}_{m}-\mbox{\boldmath$\theta$}_{n}|),$
(13)

where the indices $\alpha$ and $\beta$ run over the two components of the
reduced shear ($\alpha,\beta=1,2$), $\delta_{\alpha\beta}$ denotes the
Kronecker delta, and $\xi_{H}(x)$ is the autocorrelation of a pillbox of
radius $\theta_{\mathrm{f}}$ (White et al., 1999; Park et al., 2003; Umetsu et
al., 2015),

$\xi_{H}(x)=\frac{2}{\pi}\left[\cos^{-1}\left(\frac{x}{2\theta_{\rm
f}}\right)-\left(\frac{x}{2\theta_{\rm
f}}\right)\sqrt{1-\left(\frac{x}{2\theta_{\mathrm{f}}}\right)^{2}}\right]$
(14)

for $|x|\leqslant 2\theta_{\mathrm{f}}$ and $\xi_{H}(x)=0$ for
$|x|>2\theta_{\mathrm{f}}$.

###  2.5 Flux Magnification Bias

Lensing magnification influences the observed surface number density of
background sources behind lenses, enhancing the apparent source fluxes and
expanding the area of sky. The former effect increases the source counts above
the limiting flux, whereas the latter reduces the effective observing area in
the source plane, thus decreasing the observed source counts per unit solid
angle. The net effect, known as magnification bias (Broadhurst et al., 1995),
depends on the intrinsic slope of the source luminosity function.

Deep multi-band photometry can be used to sample the faint end of the
luminosity function for quiescent galaxies at $z\sim 1$ (Ilbert et al., 2010).
The effect of magnification bias for such a population is dominated by the
geometric area distortion, because there are relatively few fainter objects
that can be magnified into the flux-limited sample. This effect results in a
net depletion of source counts (e.g., Broadhurst et al., 2005; Ford et al.,
2012; Coe et al., 2012; Radovich et al., 2015; Ziparo et al., 2016). The key
advantage in the regime of density depletion, at the expense of deep multi-
band imaging, is that the effect is not sensitive to the exact form of the
source luminosity function (Umetsu et al., 2014).

In cluster–galaxy weak lensing, the change in magnitude $\delta
m=2.5\log_{10}\mu$ due to magnification is small compared to the range over
which the slope of the luminosity function varies. The source counts can thus
be locally approximated by a power law at a given cutoff magnitude
$m_{\mathrm{cut}}$ (Umetsu, 2020). Following Umetsu et al. (2014, 2016), we
interpret the source-averaged magnification bias as (see Appendix A)

$\displaystyle\widehat{b}_{\mu}(\mbox{\boldmath$\theta$})$
$\displaystyle\equiv\frac{\widehat{N}_{\mu}(\mbox{\boldmath$\theta$}|<m_{\mathrm{cut}})}{\overline{N}_{\mu}(<m_{\mathrm{cut}})}\simeq\Delta_{\mu}(\mbox{\boldmath$\theta$})^{1-2.5s},$
(15) $\displaystyle\Delta_{\mu}(\mbox{\boldmath$\theta$})$
$\displaystyle=\left[1-\langle
W\rangle_{\mu}\kappa_{\infty}(\mbox{\boldmath$\theta$})\right]^{2}-\langle
W\rangle_{\mu}^{2}|\gamma_{\infty}(\mbox{\boldmath$\theta$})|^{2},$

where the expected value of a weak-lensing observable is denoted by a hat
symbol,
$\overline{N}_{\mu}(<m_{\mathrm{cut}})=\int_{0}^{\infty}\\!dz\,\overline{N}_{\mu}(z|<m_{\mathrm{cut}})$
is the unlensed mean counts per cell, $\langle W\rangle_{\mu}$ is the source-
averaged relative lensing strength (Section 2.2), and $s$ is the logarithmic
count slope evaluated at the cutoff magnitude $m_{\mathrm{cut}}$,

$s(m_{\mathrm{cut}})=\frac{d\log_{10}\overline{N}(<m)}{dm}\Bigg{|}_{m_{\mathrm{cut}}}.$
(16)

Since a given magnitude cut corresponds to different luminosities at different
source redshifts, different source populations probe different regimes of
magnification bias (Umetsu, 2013). A net depletion (or enhancement) of source
counts results when $s<0.4$ (or $>0.4$). In this study, we measure the density
depletion signal using a source population with $s<0.4$. For simplicity, we
write
$N_{\mu}(\mbox{\boldmath$\theta$})=N_{\mu}(\mbox{\boldmath$\theta$}|<m_{\mathrm{cut}})$
and $\overline{N}_{\mu}=\overline{N}_{\mu}(<m_{\mathrm{cut}})$. In the weak-
lensing limit, $\widehat{b}_{\mu}-1\simeq(5s-2)\langle
W\rangle_{\mu}\kappa_{\infty}$.

The covariance matrix
$\mathrm{Cov}[N(\mbox{\boldmath$\theta$}_{m}),N(\mbox{\boldmath$\theta$}_{n})]\equiv(C_{N})_{mn}$
of the counts in cell includes the clustering and Poisson contributions,
$(C_{N})_{mn}=(\overline{N}_{\mu})^{2}\omega_{mn}+\delta_{mn}N_{\mu}(\mbox{\boldmath$\theta$}_{m})$
(Hu & Kravtsov, 2003) with $\omega_{mn}$ the cell-averaged angular correlation
function of source galaxies. As discussed in detail by Umetsu et al. (2015),
$C_{N}$ can be approximated as

$\left(C_{N}\right)_{mn}\simeq\left[\langle\delta
N_{\mu}^{2}(\mbox{\boldmath$\theta$}_{m})\rangle+N_{\mu}(\mbox{\boldmath$\theta$}_{m})\right]\delta_{mn},$
(17)

with $\langle\delta N_{\mu}^{2}(\mbox{\boldmath$\theta$}_{m})\rangle$ the
variance of the $m$th counts.

To overcome this noise, we azimuthally average the observed counts
$N_{\mu}(\mbox{\boldmath$\theta$})$ in a set of clustercentric annuli and
calculate the surface number density profile
$\\{n_{\mu,i}\\}_{i=1}^{N_{\mathrm{bin}}}$ of background galaxies as (Umetsu
et al., 2015, 2016)

$n_{\mu,i}=\frac{1}{(1-f_{\mathrm{mask},i})\Omega_{\mathrm{cell}}}\sum_{m}{\cal
P}_{im}N_{\mu}(\mbox{\boldmath$\theta$}_{m}),$ (18)

where $\Omega_{\mathrm{cell}}$ is the solid angle of each cell and ${\cal
P}_{im}=(\sum_{m}A_{mi})^{-1}A_{mi}$ is the projection matrix normalized by
$\sum_{m}{\cal P}_{im}=1$; $A_{mi}$ denotes the area fraction of the $m$th
cell lying within the $i$th radial bin and $f_{\mathrm{mask},i}$ is the mask
correction factor for the $i$th bin,
$(1-f_{\mathrm{mask},i})^{-1}\equiv\left[\sum_{m}(1-f_{m})A_{mi}\right]^{-1}\sum_{m}A_{mi}$,
with $f_{m}$ the masked area fraction in the $m$th cell due to saturated
objects, foreground galaxies, and cluster members. We use Monte Carlo
integration to calculate the area fractions $A_{mi}$ for individual cells
(Umetsu & Broadhurst, 2008). The Poisson and clustering contributions to the
uncertainty in $n_{\mu,i}$ are

$\sigma_{\mu,i}^{2}=\frac{1}{(1-f_{\mathrm{mask},i})^{2}\Omega_{\mathrm{cell}}^{2}}\sum_{m}{\cal
P}_{im}^{2}\left(C_{N}\right)_{mm}.$ (19)

Additionally, we account for systematic uncertainties in the magnification
analysis. In Appendix C, we describe the procedure used to estimate the
uncertainties $\sigma_{\mu,i}$ in $n_{\mu,i}$.

The expectation for the observable $n_{\mu,i}$ (Equation (18)) is

$\widehat{n}_{\mu,i}=\overline{n}_{\mu}\sum_{m}{\cal
P}_{im}\Delta_{\mu}(\mbox{\boldmath$\theta$}_{m})^{1-2.5s}$ (20)

with $\overline{n}_{\mu}=\overline{N}_{\mu}/\Omega_{\mathrm{cell}}$.

###  2.6 Mass Reconstruction Algorithm

A practical limitation of the shear-only lensing analysis is the inherent
mass-sheet degeneracy, which can be alleviated by using the complementary
combination of shear and magnification (Schneider et al., 2000; Umetsu &
Broadhurst, 2008; Rozo & Schmidt, 2010). Measuring the two complementary
effects also enables us to check the internal consistency of weak-lensing
measurements (Umetsu et al., 2014). Moreover, obtaining accurate mass maps has
the important advantage of being able to identify local mass structures and to
directly compare them with multiwavelength observations.

In this work, we use the mass inversion algorithm developed by Umetsu et al.
(2015), who generalized the cluster lensing mass inversion (clumi) code of
Umetsu (2013) into a 2D description of the pixelized mass distribution. This
free-form algorithm combines a 2D shear pattern
($g_{1}(\mbox{\boldmath$\theta$}),g_{2}(\mbox{\boldmath$\theta$})$) with
azimuthally averaged measurements of magnification bias
$\\{n_{\mu,i}\\}_{i=1}^{N_{\mathrm{bin}}}$. The latter imposes a set of
azimuthally averaged constraints on $\Sigma(\mbox{\boldmath$\theta$})$ to
effectively break the mass-sheet degeneracy. The clumi-2D algorithm takes full
account of the nonlinear subcritical regime of lensing.

Given a model $\lambda$ and observed data $d$, the Bayes’ theorem states that
the joint posterior probability
$P(\mbox{\boldmath$\lambda$}|\mbox{\boldmath$d$})$ is proportional to the
product of the likelihood ${\cal L}(\mbox{\boldmath$\lambda$})\equiv
P(\mbox{\boldmath$d$}|\mbox{\boldmath$\lambda$})$ and the prior probability
$P(\mbox{\boldmath$\lambda$})$. In our inversion problem, $\lambda$ is a
signal vector containing the pixelized mass coefficients
$\mbox{\boldmath$m$}=\\{\Sigma_{n}\\}_{n=1}^{N_{\mathrm{pix}}}$ (Section 2.3)
and calibration nuisance parameters $c$ (see Section 2.6.3), so that
$\mbox{\boldmath$\lambda$}\equiv(\mbox{\boldmath$m$},\mbox{\boldmath$c$})$.

We express the joint likelihood function ${\cal L}$ for combined weak-lensing
data $d$ as a product of the two separate likelihood functions, ${\cal
L}={\cal L}_{g}\times{\cal L}_{\mu}$ with ${\cal L}_{g}$ and ${\cal L}_{\mu}$
the likelihood functions for shear and magnification, respectively. We assume
that the observational errors follow a Gaussian distribution, so that ${\cal
L}\propto\exp(-\chi^{2}/2)$, with $\chi^{2}$ the standard misfit statistic.

####  2.6.1 Shear Log-likelihood Function

The log-likelihood function $l_{g}\equiv-\ln{\cal L}_{g}$ for 2D shear data is
written as (Oguri et al., 2010; Umetsu et al., 2015, 2018)

$\displaystyle l_{g}(\mbox{\boldmath$\lambda$})=$
$\displaystyle\frac{1}{2}\sum_{m,n=1}^{N_{\mathrm{pix}}}\sum_{\alpha,\beta=1}^{2}[g_{\alpha,m}-\widehat{g}_{\alpha,m}(\mbox{\boldmath$\lambda$})]\left({\cal
W}_{g}\right)_{\alpha\beta,mn}$ (21)
$\displaystyle\times[g_{\beta,n}-\widehat{g}_{\beta,n}(\mbox{\boldmath$\lambda$})]+\mathrm{const.},$

where $\widehat{g}_{\alpha,m}(\mbox{\boldmath$\lambda$})$ is the theoretical
expectation for $g_{\alpha,m}=g_{\alpha}(\mbox{\boldmath$\theta$}_{m})$ and
$({\cal W}_{g})_{\alpha\beta,mn}$ is the shear weight matrix,

$\left({\cal
W}_{g}\right)_{\alpha\beta,mn}=M_{m}M_{n}\left(C_{g}^{-1}\right)_{\alpha\beta,mn}.$
(22)

Here, $M_{m}$ is a mask weight, defined such that $M_{m}=0$ if the $m$th cell
is masked out and $M_{m}=1$ otherwise, and $C_{g}$ is the shear covariance
matrix given by Equation (13).

####  2.6.2 Magnification Log-likelihood Function

The log-likelihood function for magnification bias data
$l_{\mu}\equiv-\ln{\cal L}_{\mu}$ is written as (Umetsu et al., 2015, 2018)

$\displaystyle l_{\mu}(\mbox{\boldmath$\lambda$})=$
$\displaystyle\frac{1}{2}\sum_{i=1}^{N_{\mathrm{bin}}}[n_{\mu,i}-\widehat{n}_{\mu,i}(\mbox{\boldmath$\lambda$})]\left({\cal
W}_{\mu}\right)_{ij}$ (23)
$\displaystyle\times[n_{\mu,j}-\widehat{n}_{\mu,j}(\mbox{\boldmath$\lambda$})]+\mathrm{const.},$

where $\widehat{n}_{\mu,i}(\mbox{\boldmath$\lambda$})$ is the theoretical
expectation for $n_{\mu,i}$ and $({\cal W}_{\mu})_{ij}$ is the magnification
weight matrix, ${\cal W}_{\mu}=C^{-1}_{\mu}$, with $C_{\mu}$ the corresponding
covariance matrix,

$(C_{\mu})_{ij}=\sigma_{\mu,i}^{2}\delta_{ij},$ (24)

where the diagonal errors $\sigma_{\mu,i}$ ($i=1,2,\dots,N_{\mathrm{bin}}$)
are given by Equation (C1).

The $l_{\mu}$ function sets azimuthally integrated constraints on
$\Sigma(\mbox{\boldmath$\theta$})$, providing the otherwise unconstrained
normalization of $\Sigma(\mbox{\boldmath$\theta$})$ over a set of concentric
annuli where magnification measurements are obtained. No assumption is made
about the azimuthal symmetry of $\Sigma(\mbox{\boldmath$\theta$})$ in our
analysis. We use Monte Carlo integration to compute the projection matrix
${\cal P}_{im}$ (Equation (18)) of size $N_{\mathrm{bin}}\times
N_{\mathrm{pix}}$, which is necessary to predict
$\\{\widehat{n}_{\mu,i}(\mbox{\boldmath$\lambda$})\\}_{i=1}^{N_{\mathrm{bin}}}$
for a given model
$\mbox{\boldmath$\lambda$}=(\mbox{\boldmath$m$},\mbox{\boldmath$c$})$.

####  2.6.3 Calibration Parameters

In our joint likelihood analysis, we account for the uncertainty in the
observational calibration parameters,

$\mbox{\boldmath$c$}=\\{\langle W\rangle_{g},f_{W,g},\langle
W\rangle_{\mu},\overline{n}_{\mu},s\\},$ (25)

with $\langle W\rangle_{g}=\langle\beta\rangle_{g}/\beta_{\infty}$, $\langle
W\rangle_{\mu}=\langle\beta\rangle_{\mu}/\beta_{\infty}$, and
$f_{W,g}=\langle\beta^{2}\rangle_{g}/\langle\beta\rangle_{g}^{2}$ (Section
2.2). To this end, we include Gaussian priors on $c$ defined with mean values
and uncertainties directly estimated from data. Specifically, we use for each
parameter the mean and uncertainty estimated from the Suprime-Cam data (Tables
3 and 4) as the center and dispersion of the prior distribution, respectively.

###  2.7 Best-fit Solution and Covariance Matrix

The log-posterior function
$F(\mbox{\boldmath$\lambda$})=-\ln{P(\mbox{\boldmath$\lambda$}|\mbox{\boldmath$d$})}$
is written as a linear sum of the log-likelihood and log-prior (or quadratic
penalty) terms. The global maximum of the joint posterior probability
distribution function (PDF) over $\lambda$ is found by minimizing
$F(\mbox{\boldmath$\lambda$})$ with respect to $\lambda$. We use the
conjugate-gradient algorithm (see Press et al., 1992) to find the global
solution $\widehat{\mbox{\boldmath$\lambda$}}$. We employ an analytic
expression for the gradient function $\mbox{\boldmath$\nabla$}F$ obtained in
the nonlinear, subcritical regime (see Appendix B of Umetsu et al. 2018).

The reconstructed mass pixels are correlated primarily because the relation
between the shear and convergence is nonlocal (Equation (2)). Additionally,
the effects of spatial averaging (Equation (14)) and cosmic noise due to
projected uncorrelated large scale structure can produce a covariance between
different pixels. In our analysis, the effects of correlated errors are
modeled analytically. Specifically, we take into account the statistical and
cosmic-noise contributions to the total covariance matrix
$C_{mn}=\mathrm{Cov}(\Sigma_{m},\Sigma_{n})$
($m,n=1,2,\dots,N_{\mathrm{pix}}$) as

$C=C_{\mathrm{stat}}+C_{\mathrm{lss}},$ (26)

where $C_{\mathrm{stat}}$ is given by $(C_{\mathrm{stat}})_{mn}=\left({\cal
F}^{-1}\right)_{mn}$ with ${\cal F}$ the Fisher matrix evaluated at the best-
fit solution $\mbox{\boldmath$\lambda$}=\widehat{\mbox{\boldmath$\lambda$}}$
(see Appendix B of Umetsu et al. 2018),

${\cal
F}_{mn}=\left\langle\frac{\partial^{2}F(\mbox{\boldmath$\lambda$})}{\partial\lambda_{m}\partial\lambda_{n}}\right\rangle\Bigg{|}_{\widehat{\mbox{\boldmath$\lambda$}}},$
(27)

and
$(C_{\mathrm{lss}})_{mn}=\Sigma_{\mathrm{cr}}^{2}\xi_{\mathrm{lss}}(|\mbox{\boldmath$\theta$}_{m}-\mbox{\boldmath$\theta$}_{n}|)$
with $\xi_{\mathrm{lss}}(|\mbox{\boldmath$\theta$}|)$ the cell-averaged two-
point angular correlation function for the cosmic convergence field
$\kappa_{\mathrm{lss}}(\mbox{\boldmath$\theta$})$ (Kaiser, 1992). In this
work, we approximate the pixel window function (e.g., Hu & White, 2001) by a
Dirac delta function centered at each pixel and compute the elements of the
$C_{\mathrm{lss}}$ matrix for a given source population (see Section 3.3),
using the nonlinear matter power spectrum of Smith et al. (2003) for the
base-$\Lambda$CDM model from Planck 2018 cosmic microwave background (CMB)
anisotropy data in combination with CMB lensing (Planck Collaboration et al.,
2020, see their Table 2).

##  3 Subaru Data and Weak-lensing Analysis

In this section, we describe our new weak-lensing analysis of A370 based on
deep Suprime-Cam $BR_{\mathrm{C}}z^{\prime}$ imaging. In this study, we
analyze the Suprime-Cam data using our reduction and analysis pipelines
presented in Umetsu et al. (2014, see also ), who performed a homogeneous
weak-lensing analysis of 20 high-mass clusters targeted by the CLASH program.
As detailed in Section 3.1, the present analysis further implements an
improved astrometry based on the Gaia mission (Gaia Collaboration et al.,
2021).

###  3.1 Data and Photometry

Figure 1:  Subaru/Suprime-Cam $BR_{\mathrm{C}}z^{\prime}$ composite color image centered on A370, overlaid with mass contours from our weak-lensing analysis (Section 4). The image is $15\arcmin\times 15\arcmin$ in size. The mass map is smoothed with a Gaussian of $1.2\arcmin$ FWHM. The lowest contour level and the contour interval are $\Delta\kappa=0.08$. Also overlaid are logarithmically spaced X-ray brightness contours (cyan dashed) from Chandra observations in the $0.5$–$7$ keV energy band (Section 6). The horizontal bar represents $1\,h^{-1}\,\mathrm{Mpc}$ at the cluster redshift. North is up and east is to the left. Table 2: Subaru Suprime-Cam data Filter | Total exposure time | SeeingaaSeeing FWHM from the full stack of images. |  $m_{\mathrm{lim}}$bbLimiting magnitude for a $3\sigma$ detection in a $2\arcsec$ diameter aperture.
---|---|---|---
| (ks) | (arcsec) | (AB mag)
$B$ | $3.0$ | $0.74$ | $27.3$
$R_{\mathrm{C}}$ | $3.6$ | $0.63$ | $26.3$
$z^{\prime}$ | $8.4$ | $0.72$ | $25.8$

We analyze deep $BR_{\mathrm{C}}z^{\prime}$ images centered on A370 observed
with the wide-field camera Suprime-Cam ($34\arcmin\times 27\arcmin$; Miyazaki
et al., 2002) mounted at the prime focus of the 8.2 m Subaru Telescope.
Details of the Subaru/Suprime-Cam observations are summarized in Table 2. We
use existing archival data taken from SMOKA.222http://smoka.nao.ac.jp The
$R_{\mathrm{C}}$-band images used in this work were taken in excellent seeing
conditions on the night of 2005 December 4 (Proposal ID: o05319). The
$R_{\mathrm{C}}$ images were obtained at two different camera orientations
separated by 90 degrees. The $B$ images were taken on the night of 2010
October 12 (Proposal ID: o10314). For the $z^{\prime}$ band, we use data taken
on the nights of 2009 September 17 and 2010 October 12 (Proposal ID: o10314)
after the Suprime-Cam CCD upgrade in 2008. For the weak-lensing shape
measurements (Section 3.2), we use the $R_{\mathrm{C}}$-band data, which have
the best image quality in our data sets.

Figure 1 shows a Suprime-Cam $BR_{\mathrm{C}}z^{\prime}$ composite color image
of the cluster field, produced using the trilogy software (Coe et al., 2012).
The image is overlaid by mass contours from our weak-lensing analysis (Section
4) and X-ray brightness contours from our Chandra analysis (Section 6).

The image reduction pipeline used in this study derives from Nonino et al.
(2009). Several modifications and improvements have been applied to the
original pipeline (e.g., Umetsu et al., 2012, 2014, 2015; Medezinski et al.,
2013, 2016). In particular, it has been optimized separately for accurate
photometry and shape measurements. For multi-band photometry, standard
reduction steps include bias subtraction, super-flat-field correction, and
masking of saturated star trails and other artifacts. Photometric catalogs are
created using sextractor (Bertin & Arnouts, 1996) in the dual-image mode on
PSF-matched images, with the Suprime-Cam $z^{\prime}$ band image as the
detection image.

An accurate astrometric solution was derived with the scamp software (Bertin,
2006) using Gaia Data Release 2 (DR2; Gaia Collaboration et al., 2016, 2018)
as an external reference catalog. An astrometric solution has been obtained at
the camera level using Gaia DR2 sources extracted from individual exposures
for each CCD chip. This astrometric solution does not account for the proper
motions of Gaia DR2 sources since the epoch of the Suprime-Cam observations.
Comparing with the astrometric solution obtained from proper-motion-corrected
source positions based on Gaia Early Data Release 3 (EDR3; Gaia Collaboration
et al., 2021), we find a mean positional offset of $\approx 2$–$8$ mas and an
rms of $\approx 20$–$30$ mas for the $R_{\mathrm{C}}$-band astrometry.

The swarp software (Bertin et al., 2002) is used to stack individual exposures
on a common World Coordinate System (WCS) grid with pixel scale of
$0.2\arcsec$. No point spread function (PSF) matching is applied. For each
passband, we create a full stack of co-added images from which to measure
source photometry. For the weak-lensing band ($R_{\mathrm{C}}$), we
additionally create two separate co-added images, each from different camera
rotation angles (see Section 3.2).

Once the Subaru images had been combined into a full stack of co-added images,
a catalog was then produced and matched directly to the corresponding Gaia DR2
sources to validate the astrometric properties of the full stack. A total of
428 sources from the full stack catalog were matched directly to Gaia DR2
sources, and the positional offsets were measured between each of these
catalog sources and their matched Gaia DR2 counterparts. The resulting
distribution of positional offsets displays well-behaved symmetry, with an rms
uncertainty of 35 mas in R.A. and 33 mas in decl., demonstrating that the
positional accuracy of the full stack is in good agreement with the accuracy
of the Gaia DR2 alignment carried out on each of the individual single-
exposure frames that were used to construct the full stack, as previously
described.

Finally, after having verified that the positional uncertainties of the full
stack catalog sources were comparable to those of all the individual single-
exposure frames, this full stack catalog was then aligned to Gaia EDR3, to
ensure that the absolute astrometry could be as up-to-date as possible. The
astrometric difference between sources from the full-stack catalog (which was
still on Gaia DR2) and the matched sources from Gaia EDR3 show excellent
agreement, with only a small difference needing to be applied to place these
sources onto Gaia EDR3, namely 2.1 mas in R.A. and 1.4 mas in decl., perhaps
due to slight residual differences in proper motion corrections, and not
significant compared to the rms uncertainties of $33$–$35$ mas in the catalog
source positions.

The photometric zero point for the Suprime-Cam $z^{\prime}$ filter was
calibrated against stars from the Pan-STARRS data release 1 (DR1) catalog
(Flewelling et al., 2020). The zero points for the Suprime-Cam $B$ and
$R_{\mathrm{C}}$ filters were derived by matching the stellar locus in the
$B-R_{\mathrm{C}}$ vs. $R_{\mathrm{C}}-z^{\prime}$ diagram to the COSMOS2020
photometry (Weaver et al., 2022). These zero points were further refined by
matching the color distributions in the $B-R_{\mathrm{C}}$ vs.
$R_{\mathrm{C}}-z^{\prime}$ diagram between our Suprime-Cam data and the
COSMOS2020 data. The magnitudes for galaxies were corrected for foreground
Galactic extinction according to Schlegel et al. (1998). Full details of our
photometric calibration are described in Appendix B.

###  3.2 Shape Measurement

We use our shape measurement pipeline based in part on the imcat package
(Kaiser et al., 1995, KSB), with modifications incorporating several key
improvements developed by Umetsu et al. (2010, 2014). In this work, we perform
a weak shear analysis of A370 following the procedure of Umetsu et al. (2014).

Here we briefly summarize some of the main features and refer to Umetsu et al.
(2014) for details of the analysis pipeline. We select isolated galaxy images
for the shape measurement, reducing the impact of crowding and blending. After
the rejection of close pairs, objects detected with low significance
$\nu_{g}<10$ are excluded from our analysis. Here $\nu_{g}$ is the peak
detection significance given by imcat’s peak-finding algorithm. We select
galaxies detected with high significance $\nu_{g}\geqslant 30$ as a sample of
shape calibrators, which is a subset of the target galaxy sample with
$\nu_{g}\geqslant 10$. The key feature of our analysis pipeline is that only
those galaxies detected with sufficiently high significance, $\nu_{g}\geqslant
30$, are used to model the isotropic PSF correction as a function of object
size and magnitude (Umetsu et al., 2010). This calibration method is designed
to minimize the inherent noise bias and was employed by the CLASH and LoCuSS
collaborations in their cluster weak-lensing studies based on Subaru/Suprime-
Cam data (Umetsu et al., 2014; Okabe & Smith, 2016).

For the shape measurement, we separately stack $R_{\mathrm{C}}$-band images
collected at two different camera rotation angles (Section 3.1). In this way,
we do not smear individual exposures before stacking, so as not to degrade the
weak-lensing signal derived from the shapes of galaxies (Umetsu et al., 2014,
2015). A shape catalog is created for each camera rotation separately. The two
subcatalogs are combined by properly weighting and stacking the calibrated
distortion measurements for galaxies in the overlapping region (Umetsu et al.,
2014, see their Section 4.3). All galaxies with usable shape measurements are
matched to those in our $BR_{\mathrm{C}}z^{\prime}$-selected background
samples (see Section 3.3).

Our KSB+ implementation has been extensively tested and applied to ground-
based observations of a large number of massive clusters including 20 CLASH
clusters (Umetsu et al., 2014; Merten et al., 2015). Full details of our shear
recovery test based on simulated Subaru/Suprime-Cam images are found in Umetsu
et al. (2018, see their Appendix A). They found that the reduced shear signal
$g_{\alpha}$ ($\alpha=1,2$) can be recovered with $m_{\alpha}\approx-0.05$ of
the multiplicative calibration bias and $|c_{\alpha}|\sim 10^{-4}$ of the
additive shear bias. Here the observed and true values of the reduced shear
($g^{\mathrm{obs}},g^{\mathrm{true}}$) are related by (Heymans et al., 2006;
Massey et al., 2007),

$g_{\alpha}^{\mathrm{obs}}=(1+m_{\alpha})g_{\alpha}^{\mathrm{true}}+c_{\alpha}$
(28)

Accordingly, we include for each galaxy a shear calibration factor of $1/0.95$
($g\to g/0.95$) to account for the residual multiplicative bias.

###  3.3 Background Galaxy Selection

$\begin{array}[]{c@{\hspace{0.3in}}c}\includegraphics[width=195.12767pt,angle={0},clip]{Figs//f2a.pdf}\hfil\hskip
21.68121pt&\includegraphics[width=195.12767pt,angle={0},clip]{Figs//f2b.pdf}\\\
\end{array}$

Figure 2: Binned distribution of galaxies in color–color space for the COSMOS
field (left) and A370 (right). Color boundaries of the blue and red background
samples (left blue and lower-right red regions, respectively) selected on the
basis of Subaru $BR_{\mathrm{C}}z^{\prime}$ photometry are indicated in each
panel. In the right panel, the green polygon marks the boundaries of our green
sample dominated by red-sequence galaxies of A370 at $z=0.375$. The middle
peak with colors bluer than the cluster sequence shows the overdensity of
foreground galaxies (see also Figure 3). The plots in both panels are limited
to $z^{\prime}<26$ mag, which is close to our detection limit (Table 2).
Figure 3: Binned photometric redshift of COSMOS field galaxies displayed in
color–color space. Blue and red polygons mark the boundaries of the blue and
red background samples, respectively.

Contamination of background galaxy samples by unlensed objects, when not
accounted for, leads to a systematic underestimation of the true lensing
signal. Inclusion of foreground galaxies produces a dilution of the lensing
signal that is independent of the cluster radius. In contrast, the inclusion
of cluster members dilutes the lensing signal more strongly at smaller cluster
radii (Broadhurst et al., 2005). A secure selection of background galaxies is
thus essential for obtaining accurate cluster mass estimates from weak lensing
(e.g., Medezinski et al., 2010; Okabe et al., 2013; Gruen et al., 2014).

In this study, we employ the color–color (CC) selection method of Medezinski
et al. (2010) (see also Medezinski et al., 2018) to define background galaxy
samples for measuring both shear and magnification effects. We use
$BR_{\mathrm{C}}z^{\prime}$ photometry from Subaru/Suprime-Cam, which spans
the full optical wavelength range. The CC-cut selection method has been
calibrated with evolutionary color tracks of galaxies (Kotulla et al., 2009;
Medezinski et al., 2010, 2011) as well as with photometric-redshift
(photo-$z$) catalogs from deep multiwavelength surveys such as COSMOS (Ilbert
et al., 2009; Laigle et al., 2016; Weaver et al., 2022). For this purpose, we
use the photometric properties and redshifts derived from the COSMOS2020
catalog (Weaver et al., 2022) based on the farmer photometry using the lephare
code (Ilbert et al., 2006).

In Figure 2, we show the distribution of galaxies in the $B-R_{\mathrm{C}}$
vs. $R_{\mathrm{C}}-z^{\prime}$ plane obtained for the COSMOS field (left
panel) and A370 (right panel). Similarly, Figure 3 shows the binned average
photo-$z$ distribution of COSMOS field galaxies in CC space. As demonstrated
by Medezinski et al. (2010, 2011), the color region dominated by the
foreground population is well defined in CC space as a clear overdensity (at
$B-R_{\mathrm{C}}\sim 1$ and $R_{\mathrm{C}}-z^{\prime}\sim 0.3$ with $\langle
z\rangle\sim 0.5$). Following Medezinski et al. (2010, 2011), we select two
distinct populations that encompass the “red” and “blue” branches of
background galaxies in CC space, each with typical redshift distributions
$N(z)$ peaked around $z\sim 1$ and $\sim 2$, respectively (see Medezinski et
al., 2011; Lilly et al., 2007).

The color boundaries of our CC-cut samples are shown in Figure 2. The green
polygon shown in the right panel marks the boundaries of the “green” sample
comprising mostly the red-sequence galaxies of A370. We see in Figure 2 that
the foreground peak for A370 is more pronounced compared to the COSMOS field.
This enhancement is likely due to the contribution from bluer cluster members
and galaxies in the surrounding regions (see Umetsu et al., 2012, 2015).

Table 3: Background Galaxy Samples for Weak-lensing Shape Measurements

Sample | $N_{g}$ | $\overline{n}_{g}$ | $\langle\beta\rangle$ | $f_{W}$ | $z_{\mathrm{eff}}$ | $\mathrm{S/N}$
---|---|---|---|---|---|---
|  | (arcmin-2) |  |  |  |
Red | 9988 | $12.6$ | $0.56\pm 0.03$ | $1.05$ | $1.02$ | $13.5$
Blue | 6679 | $8.4$ | $0.64\pm 0.03$ | $1.10$ | $1.33$ | $10.9$
Blue+red | 16667 | $21.0$ | $0.59\pm 0.03$ | $1.08$ | $1.12$ | $16.8$

Note. — Subaru $BR_{\mathrm{C}}z^{\prime}$ selected samples of background
galaxies. We use the composite blue+red background sample for our weak-lensing
shear analysis. The mean lensing depth $\langle\beta\rangle$ and the spread
parameter $f_{W}=\langle\beta^{2}\rangle/\langle\beta\rangle^{2}$ for each
source population are estimated using photometric redshifts from the
COSMOS2020 farmer catalog. The quantity $z_{\mathrm{eff}}$ represents the
effective source redshift of each sample, defined as
$\beta(z_{\mathrm{eff}})=\langle\beta\rangle$. The $\mathrm{S/N}$ is the
detection significance for the tangential distortion profile $g_{+}(\theta)$.

Table 4: Background Galaxy Samples for Magnification-bias Measurements

Sample | Magnitude limits | $N_{\mu}$ | $\overline{n}_{\mu}$ | $s$ | $\langle\beta\rangle$ | $z_{\mathrm{eff}}$ | $\mathrm{S/N}$
---|---|---|---|---|---|---|---
| (AB mag) |  | (arcmin-2) |  |  |  |
Lensing cut | $21.0<z^{\prime}<25.6$ | 22142 | $21.1\pm 0.6$ | $0.168\pm 0.037$ | $0.57\pm 0.03$ | $1.04$ | $5.7$
Null-test | $21.0<z^{\prime}<23.6$ | 6344 | $5.6\pm 0.3$ | $0.404\pm 0.069$ | $0.57\pm 0.03$ | $1.04$ | $2.7$

Note. — Lensing-cut and null-test samples of CC-red background galaxies
selected for our weak-lensing magnification analysis. Apparent magnitude cuts
are applied in the reddest CC-selection band available ($z^{\prime}$) to avoid
incompleteness near the detection limit (Table 2). The mean lensing depth
$\langle\beta\rangle$ for each source population is estimated using
photometric redshifts from the COSMOS2020 farmer catalog. The quantity
$z_{\mathrm{eff}}$ represents the effective source redshift corresponding to
the mean lensing depth $\langle\beta\rangle$ of each sample, defined as
$\beta(z_{\mathrm{eff}})=\langle\beta\rangle$. The $\mathrm{S/N}$ is the
detection significance for the magnification bias profile
$b_{\mu}(\theta)=n_{\mu}(\theta_{)}/\overline{n}_{\mu}$.

To further reduce residual contamination by bright foreground objects, we
apply bright magnitude cuts of $z^{\prime}>21$ and $22$ mag for the red and
blue photometry samples, respectively (Medezinski et al., 2010, 2018). These
selection criteria yield a total of $31306$ and $14954$ galaxies for the red
and blue photometry samples, respectively. For our shear analysis, we use the
weak-lensing-matched, blue and red composite sample containing $16667$
galaxies with usable $R_{\mathrm{C}}$ shape measurements, corresponding to a
mean surface number density of $\overline{n}_{g}\approx 21$ galaxies arcmin-2
(Table 3).

To measure the magnification bias, we use magnitude-limited samples of CC-red
galaxies. For the measurement of density depletion (Section 2.5), we define a
“lensing-cut” sample by applying a faint magnitude cut of $z^{\prime}=25.6$
mag to the red photometry sample (Table 4).333Our CC-cut selection is not
expected to cause incompleteness at the faint end in the bluer filters (see
Hildebrandt et al. 2012 for a general discussion) because we have deeper
photometry in the bluer bands (Broadhurst et al. 1995) and our CC-red galaxies
are relatively blue in $B-R_{\mathrm{C}}$ (Figure 2). On the other hand, since
the net effect of magnification bias is expected to vanish at $s=0.4$ (Section
2.5), lensing magnification also provides a null test, which allows us to
assess the level of residual bias that could be present in the measurement for
the lensing-cut sample (see Chiu et al., 2020; Umetsu, 2020). To this end, we
define a “null-test” sample with a faint magnitude cut of $z^{\prime}=23.6$
mag, at which the count slope is found to be $s=0.404\pm 0.069$ (Table 4).

###  3.4 Lensing Depth Estimation

To assess the mean lensing depth
($\langle\beta\rangle,\langle\beta^{2}\rangle$; see Equation (6)) for our CC-
cut samples, we use the COSMOS2020 farmer catalog with robust photometry and
photo-$z$ measurements. For each background sample, we apply the same cuts to
the COSMOS multi-band photometry and obtain the redshift distribution $N(z)$
of the selected galaxies. The lensing weight $w_{g}$ (see Section 2.4) is not
taken into account in the depth estimation, because there are no photo-$z$
estimates available for our background sample in the A370 field.444The effect
of neglecting the lensing weight $w_{g}$ was checked using photo-$z$ and shape
catalogs based on Suprime-Cam 5-band imaging available for CLASH clusters at
similar redshifts, $z_{l}\in[0.35,0.40]$ (Umetsu et al., 2014). The fractional
differences in the estimated $\langle\beta\rangle$ values are found to be
$<1\%$, which is not significant compared to the total fractional uncertainty
of $5\%$ adopted in this study. The resulting depth estimates for our shear
and magnification analyses are summarized in Tables 3 and 4,
respectively.555The expected contribution of foreground galaxies with
$\beta(z_{l},z_{s})=0$ is accounted for in our lensing depth estimation (see
Equation (6)).

For a consistency check, we also make use of photo-$z$ estimates from
alternative aperture-based COSMOS2020 photometry, Classic (Weaver et al.,
2022). For each sample, we obtain consistent depth estimates (to within $1\%$)
from the farmer and classic catalogs. Taking into account the field-to-field
variance in $N(z)$ (Umetsu et al., 2014, see their Section 4.4), we assume a
fractional uncertainty of $5\%$ in the COSMOS-based estimates of
$\langle\beta\rangle$. We marginalize over this uncertainty in our mass
reconstruction (Section 4).

The level of residual cluster contamination for the CC-cut method has been
assessed by Umetsu et al. (2016) using large spectroscopic samples from the
CLASH-VLT program (Rosati et al., 2014). Combining VLT spectroscopic redshifts
and Subaru multi-band photometry available for 10 southern CLASH clusters with
a mean redshift of $z\approx 0.37$, Umetsu et al. (2016) found a mean
contamination fraction of $(2.4\pm 0.7)\%$ in the blue+red CC-cut sample. This
level of residual contamination is subdominant compared to other uncertainties
in our lensing analysis.

###  3.5 Null Tests

Figure 4: Azimuthally averaged radial profiles of the tangential and cross
shear components, $g_{+}$ (upper panel) and $g_{\times}$ (lower panel),
respectively, for our red (triangles), blue (circles), green (crosses), and
blue+red (squares) galaxy samples.

Figure 4 shows the azimuthally averaged tangential ($g_{+}$) and cross
($g_{\times}$) components of the reduced shear as a function of projected
cluster radius. We find a rising $g_{+}(\theta)$ profile toward the cluster
center from both blue and red background samples. In contrast, the
$g_{+}(\theta)$ signal for the green sample is suppressed by the inclusion of
cluster members and consistent with zero at $R\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}2\,h^{-1}\,\mathrm{Mpc}$, while it
becomes comparable to the pure background signal outside the cluster region.

In the absence of higher-order effects, weak lensing only produces tangential
shape distortions (Section 2.1). The presence of $\times$ distortions can thus
be used to check for systematic errors. Here we use a $\chi^{2}$ test to
assess the statistical significance of the measured $\times$-mode signal
against the null hypothesis. We find $\chi^{2}$ values of the null hypothesis
to be $10.6,16.4,9.8$, and $14.9$ for $N_{\mathrm{bin}}=12$ degrees of
freedom, for the red, blue, green, and blue+red samples, respectively. For all
the cases tested, the $\times$-component signal is statistically consistent
with a null detection.

Figure 5: Coverage- and mask-corrected surface number density profiles
$n_{\mu}$ of $BR_{\mathrm{C}}z^{\prime}$-selected red background samples. The
results are shown for our lensing-cut (circles) and null-test (crosses)
samples. The error bars include both Poisson and clustering contributions
estimated from the data. For the lensing-cut sample, a radial count depletion
due to magnification of the sky area is seen toward the cluster center. For
the null-test sample with $s\approx 0.4$, the net effect of magnification bias
is expected to vanish. The mean background levels estimated for the lensing-
cut and null-test samples are marked with solid and dashed horizontal lines,
respectively.

Figure 5 shows the coverage- and mask-corrected surface number density of
background galaxies as a function of projected cluster radius, for the
lensing-cut and null-test samples. In both cases, no clustering is observed
toward the center, demonstrating that there is no detectable contamination by
cluster members. The lensing-cut sample reveals a systematic decrease in their
counts toward the cluster center, caused by magnification of the sky area. In
contrast, the null-test sample shows no significant evidence for radial count
variations with $\chi^{2}=10.0$ for $12$ degrees of freedom, as expected by
their count slope. A more quantitative magnification analysis will be
discussed in Section 4.

##  4 Weak-lensing Mass Reconstruction

###  4.1 Mass Profile Reconstruction (WL-1D)

Figure 6:  Azimuthally averaged cluster lensing profiles of A370 obtained from
Subaru observations. The upper panel shows the reduced tangential shear
profile $g_{+}$ (blue squares) based on the blue+red background sample. The
lower panel shows the magnification bias profile $n_{\mu}$ measured from our
$z^{\prime}$-limited sample of red background galaxies, with (red circles) and
without (orange dots) the mask correction. For each observed profile, the
shaded area represents the $1\sigma$ confidence region of the joint
reconstruction from the shear and magnification profile measurements. The
horizontal bar (gray shaded region) shows the estimated mean background level.
Figure 7: Surface mass density profile $\Sigma(\theta)$ (upper panel, red
squares) derived from a joint analysis of azimuthally averaged shear and
magnification measurements shown in Figure 6. The lower panel shows the
corresponding cumulative mass profile $M_{\mathrm{2D}}(<\theta)$ (red
squares). The innermost bin is an integrated average inside
$\theta_{\mathrm{min}}=1.3\arcmin$. The gray shaded area in each panel
represents the marginalized $1\sigma$ confidence interval of the NFW fit to
the $\Sigma$ profile.

Before carrying out a 2D mass reconstruction, we first perform a weak-lensing
1D radial profile analysis (WL-1D) of our Subaru observations (Section 3).
A370 has two central BCGs separated by $\approx 37\arcsec$ (about
$140\,h^{-1}\,\mathrm{kpc}$ at $z=0.375$) along the north–south direction
(Figure 1). In this work, we adopt the optical center, or the midpoint of the
two BCGs (see Table 1), as the cluster center for our radial profile analysis.

We derive azimuthally averaged radial profiles of tangential reduced shear
($g_{+}$) and magnification bias ($n_{\mu}$) from our Subaru/Suprime-Cam data.
We calculate the binned lensing profiles,
$\\{g_{+,i}\\}_{i=1}^{N_{\mathrm{bin}}}$ and
$\\{n_{\mu,i}\\}_{i=1}^{N_{\mathrm{bin}}}$, in $N_{\mathrm{bin}}=12$
logarithmically spaced radial bins centered on the cluster, spanning the range
from $\theta_{\mathrm{min}}=1.3\arcmin$ to $\theta_{\mathrm{max}}=16\arcmin$,
with a logarithmic spacing of $\Delta\ln\theta\approx 0.21$. Our radial
profile analysis begins at $\theta_{\mathrm{min}}=1.3\arcmin$, which is
sufficiently large compared to twice the effective Einstein radius,
$2\theta_{\mathrm{Ein}}\sim 1.1\arcmin$ (for $z_{s}=2$; Table 1), determined
from strong-lens modeling by Kawamata et al. (2018). Hence, our analysis does
not include outer multiple images of strongly lensed galaxies lying at
$\theta\sim[\theta_{\mathrm{Ein}},2\theta_{\mathrm{Ein}}]$, and our data do
not resolve the central substructures. The outer boundary
$\theta_{\mathrm{max}}=16\arcmin$ ($\approx 3.5\,h^{-1}\,\mathrm{Mpc}$) is
large enough to encompass the entire cluster region with $r_{\mathrm{vir}}\sim
2\,h^{-1}\,\mathrm{Mpc}$ (Umetsu et al., 2011), but sufficiently small
compared to the size of the Suprime-Cam field of view so as to ensure accurate
PSF correction.

For the magnification analysis, the count normalization parameter
$\overline{n}_{\mu}$ is estimated in the reference background region at
$\theta\in[12\arcmin,16\arcmin]$.666The 2-halo term ($\kappa\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}10^{-2}$) does not cause bias in
the mass reconstruction, because the range of the prior on
$\overline{n}_{\mu}$ is sufficiently wide. See Umetsu et al. (2014, their
Section 7.4.2) for detailed discussion. The estimated values and errors for
$\overline{n}_{\mu}$ and $s$ are summarized in Table 4. Details of the error
analysis and the mask correction procedure are described in Appendix C.

We reconstruct the radial mass profile of A370 from a joint likelihood
analysis of azimuthally averaged shear and magnification constraints, using
the clumi code of Umetsu (2013). We have a total of 24 constraints
$\\{g_{+,i},n_{\mu,i}\\}_{i=1}^{N_{\mathrm{bin}}}$ in 12 radial bins. The
model is described by $N_{\mathrm{bin}}+1=13$ parameters,
$\mbox{\boldmath$m_{\mathrm{1D}}$}=\\{\Sigma_{\mathrm{min}},\Sigma_{i}\\}_{i=1}^{N_{\mathrm{bin}}}$,
where $\Sigma_{\mathrm{min}}\equiv\Sigma(<\theta_{\mathrm{min}})$ is the
average surface mass density interior to $\theta_{\mathrm{min}}$777The central
surface mass density $\Sigma(<\theta_{\mathrm{min}})$ can be fully determined
by the combination of tangential-shear and magnification measurements outside
$\theta_{\mathrm{min}}$ (see Umetsu, 2013). and $\Sigma_{i}$ is the surface
mass density averaged in the $i$th bin. In addition, we account for the
calibration uncertainty in the observational parameters $c$ (Equation (25);
see Tables 3 and 4).888The clumi algorithm employs uniform priors on the
calibration parameters $c$ and explores parameter space with a Markov chain
sampling method (Umetsu, 2013), whereas clumi-2D uses Gaussian priors on $c$
to find the global solution $\widehat{\mbox{\boldmath$\lambda$}}$ with the
conjugate-gradient method (Section 2.7). Following Umetsu et al. (2014), we
fix $f_{W,g}$ to the estimated value.

Figure 6 compares the observed lensing profiles
$\\{g_{+,i},n_{\mu,i}\\}_{i=1}^{N_{\mathrm{bin}}}$ with the respective joint
reconstructions. The joint solution has a $\chi^{2}$ value of $18.7$ for $11$
degrees of freedom, indicating a slight (but statistically not significant)
discrepancy between the two data sets. We see from the lower panel of Figure 6
that the measured $n_{\mu}$ value at $\theta\sim 5\arcmin$ is $\sim 2\sigma$
lower than expected from the joint reconstruction. This is consistent with the
result for the null-test sample, which exhibits a similar local deficit of the
galaxy counts in the same radial bin (see Figure 5). For the other bins, we
find a good agreement between the shear and magnification data. The
reconstructed $\Sigma(\theta)$ profile is shown in the upper panel of Figure
7, along with the $1\sigma$ confidence interval of the spherical
Navarro–Frenk–White (Navarro et al., 1996, 1997, hereafter NFW) model (see
Section 5 for details of the modeling). The corresponding cumulative mass
profile $M_{\mathrm{2D}}(<\theta)=\pi(D_{l}\theta)^{2}\Sigma(<\theta)$ is
shown in the lower panel of Figure 7.

###  4.2 Two-dimensional Map Making (WL-2D)

Figure 8:  Spatial distribution of reduced shear constraints ($g_{1},g_{2}$)
averaged onto a mass grid of $48\times 48$ pixels, covering a field of
$24\arcmin\times 24\arcmin$ centered on A370. Each point represents a single
pixel with measured ($g_{1},g_{2}$) averaged within a top-hat region with
radius $\theta_{\mathrm{f}}=0.4\arcmin$. We exclude from our analysis those
pixels lying within the central $1\arcmin$ region (red circle) and those
having no background galaxies with usable shape measurements. Azimuthally
averaged magnification constraints are obtained in $12$ logarithmically spaced
annuli centered on the cluster spanning the range
$\theta\in[1.3\arcmin,16\arcmin]$.  Figure 9: Projected mass distribution
$\Sigma(\mbox{\boldmath$\theta$})$ of A370 reconstructed from a joint analysis
of 2D shear and azimuthally averaged magnification data (Figure 8). The map is
$24\arcmin\times 24\arcmin$ in size ($5.2\,h^{-1}\,\mathrm{Mpc}$ on a side at
$z=0.375$) and centered on the cluster. The color bar indicates
$\kappa=\Sigma_{\mathrm{cr}}^{-1}\Sigma$ scaled to the mean lensing depth,
$1/\langle\Sigma_{\mathrm{cr}}^{-1}\rangle=3.77\times
10^{15}~{}h\,M_{\odot}\,\mathrm{Mpc}^{-2}$. For visualization purposes, the
mass map is smoothed with a $3\times 3$ pixel boxcar kernel. The dashed
contours show the surface density distribution of red-sequence cluster
galaxies, smoothed with a Gaussian of $1.2\arcmin$ FWHM. The lowest contour
level and the contour interval are both $10\%$ of the peak density $n\approx
47$ galaxies arcmin-2. The red circle indicates the cluster radius of
$r_{200}\approx 1.7\,h^{-1}\,\mathrm{Mpc}$. North is to the top, east to the
left.

We apply our clumi-2D method (Section 2.6) to our Subaru/Suprime-Cam data
(Section 3) for obtaining an unbiased recovery of the projected mass
distribution $\Sigma(\mbox{\boldmath$\theta$})$ in A370. In this analysis
(WL-2D), we combine the observed shear field
$(g_{1}(\mbox{\boldmath$\theta$}),g_{2}(\mbox{\boldmath$\theta$}))$ with the
azimuthally averaged magnification data
$\\{n_{\mu,i}\\}_{i=1}^{N_{\mathrm{bin}}}$ (Section 4.1), which impose a set
of azimuthally integrated constraints on the underlying
$\Sigma(\mbox{\boldmath$\theta$})$ field. clumi-2D takes into account the
nonlinear subcritical regime of the lensing properties.

For mapmaking, we pixelize the $\gamma_{\infty}(\mbox{\boldmath$\theta$})$ and
$\kappa_{\infty}(\mbox{\boldmath$\theta$})$ fields into a uniform grid of
$48\times 48$ pixels with $\Delta\theta=0.5\arcmin$ spacing, covering the
central $24\arcmin\times 24\arcmin$ field. The model
$\mbox{\boldmath$\lambda$}=(\mbox{\boldmath$m$},\mbox{\boldmath$c$})$ is
specified by $N_{\mathrm{pix}}=48^{2}$ parameters,
$\mbox{\boldmath$m$}=\\{\Sigma(\mbox{\boldmath$\theta$}_{n})\\}_{n=1}^{N_{\mathrm{pix}}}$,
and a set of calibration parameters $c$ to marginalize over (Tables 3 and 4).
We utilize the fftw implementation of fast Fourier transforms to compute
$\gamma_{\infty}(\mbox{\boldmath$\theta$})$ from
$\kappa_{\infty}(\mbox{\boldmath$\theta$})$ using Equation (9). To avoid
spurious aliasing effects from the periodic boundary condition, the maps are
zero padded to twice the data length in each spatial dimension (e.g., Seljak,
1998; Umetsu & Broadhurst, 2008).

We use a top-hat window of $\theta_{\mathrm{f}}=0.4\arcmin$ (Section 2.4) to
average over a local ensemble of galaxy image ellipticities at each grid
point, accounting for the intrinsic ellipticity distribution of background
galaxies. To avoid potential systematic errors (see Section 2.6.1), we exclude
from our analysis 12 pixels lying within the central
$\theta_{\mathrm{cut}}=1\arcmin$ and one pixel containing no background
galaxies. For distortion measurements
($g_{1}(\mbox{\boldmath$\theta$}),g_{2}(\mbox{\boldmath$\theta$})$), this
leaves us with a total of $2291$ usable measurement pixels (blue points in
Figure 8), corresponding to $4582$ constraints. For magnification
measurements, we have $12$ azimuthally averaged constraints
$\\{n_{\mu,i}\\}_{i=1}^{N_{\mathrm{bin}}}$ (Figure 8). The total number of
constraints is thus $N_{\mathrm{data}}=4594$, yielding
$N_{\mathrm{data}}-N_{\mathrm{pix}}=2290$ degrees of freedom.

In Figure 9, we show the $\Sigma(\mbox{\boldmath$\theta$})$ field
reconstructed from the joint analysis of the 2D shear and azimuthally averaged
magnification data. The $\chi^{2}$ value for the global maximum posterior
solution is $\chi^{2}(\widehat{\mbox{\boldmath$\lambda$}})=2871$ for $2290$
degrees of freedom. For comparison, we plot in Figure 9 the surface density
distribution of the green sample (dashed contours) composed mostly of cluster
members. The projected mass distribution is elongated in the north–south
direction and similar to that of cluster member galaxies (Figure 9). Our mass
reconstruction barely resolves substructure features (e.g., a north–south mass
extension located about $1\arcmin$ north and south of the cluster center)
revealed by the free-form mass inversion of Ghosh et al. (2021) based on
BUFFALO strong-lensing data. We defer a more detailed investigation of weak-
lensing substructures in the A370 field to a forthcoming paper (S.-I. Tam et
al. 2022, in preparation).

Table 5: Projected total mass estimates for A370 from the WL-2D analysis Aperture radius, $\theta$ |  $r_{\perp}$aaClustercentric radius in physical units, $r_{\perp}=D_{l}\theta$. | $M_{\mathrm{2D}}(<\theta)$
---|---|---
(arcmin) | ($\,h^{-1}\,\mathrm{Mpc}$) | ($10^{15}\,h^{-1}\,M_{\odot}$)
$1.3$ | $0.28$ | $0.51\pm 0.05$
$1.6$ | $0.35$ | $0.52\pm 0.05$
$2.0$ | $0.43$ | $0.58\pm 0.06$
$2.4$ | $0.53$ | $0.82\pm 0.07$
$3.0$ | $0.65$ | $0.99\pm 0.09$
$3.7$ | $0.80$ | $1.16\pm 0.11$
$4.6$ | $0.99$ | $1.37\pm 0.14$
$5.6$ | $1.22$ | $1.61\pm 0.18$
$6.9$ | $1.50$ | $2.09\pm 0.24$
$8.5$ | $1.85$ | $2.61\pm 0.33$
$10.5$ | $2.28$ | $3.11\pm 0.47$
$13.0$ | $2.81$ | $3.69\pm 0.68$
$16.0$ | $3.47$ | $4.30\pm 1.02$

We construct the binned radial profiles $\Sigma(\theta)$ and $\Sigma(<\theta)$
and their associated covariance matrices from an optimally weighted projection
of the $\Sigma(\mbox{\boldmath$\theta$})$ map using the method described in
Appendix D. We thus obtain model-independent constraints on the projected
total mass $M_{\mathrm{2D}}(<\theta)=\pi(D_{l}\theta)^{2}\Sigma(<\theta)$ from
our WL-2D analysis. The resulting projected mass estimates are listed in Table
5.

###  4.3 Radial Mass Profiles

Figure 10:  Comparison of reconstructed surface mass density profiles
$\Sigma(\theta)$ of A370. The black circles (red squares) with error bars show
the results from our 2D (1D) weak-lensing reconstruction. The blue shaded area
represents the $\Sigma$ profile with $2\sigma$ uncertainty derived from strong
lens modeling of Hubble Frontier Fields data (Kawamata et al., 2018). The
orange diamonds with error bars show the results from our earlier Subaru weak-
lensing analysis (Umetsu et al., 2011). For each case, the innermost central
bin $\Sigma(<\theta_{\mathrm{min}})$ is marked with a horizontal bar. The gray
shaded area represents the marginalized $1\sigma$ confidence interval of the
spherical NFW fit to the $\Sigma$ profile from our 1D reconstruction.

In Figure 10, we compare the surface mass density profiles $\Sigma(\theta)$ of
A370 obtained from our WL-1D (Section 4.1) and WL-2D (Section 4.2) analyses.
Our 1D- and 2D-based $\Sigma$ profiles are consistent within the errors in
each radial bin. The gray shaded area in the figure represents the $1\sigma$
confidence region of the spherical NFW fit to the 1D-based $\Sigma$ profile
(see Section 5 for details of the modeling).

For comparison, we overplot in Figure 10 the azimuthally averaged $\Sigma$
profile (shown out to $2\theta_{\mathrm{Ein}}\approx 1.1\arcmin$ for
$z_{s}=2$) based on strong lens modeling of Hubble Frontier Fields data
performed by Kawamata et al. (2018, see also ), obtained using the technique
detailed in Oguri (2021) to speed up lensing calculations. The inner $\Sigma$
profile derived from HST strong lensing is in excellent agreement with our
WL-1D constraints on the NFW profile.

In Figure 10, our $\Sigma$ profiles are also compared with the 1D results of
Umetsu et al. (2011) based on their joint shear and magnification analysis of
Suprime-Cam data. In the Umetsu et al. (2011) analysis, the innermost
measurement radius was taken to be $\theta_{\mathrm{min}}=0.7\arcmin$
($\approx 1.2\theta_{\mathrm{Ein}}$ for $z_{s}=2$), in contrast to the
conservative choice adopted in this work ($\theta_{\mathrm{min}}=1.3\arcmin$).
The shear and magnification measurements of Umetsu et al. (2011) are based on
Suprime-Cam images reduced by Medezinski et al. (2010), who used the sdfred
package (Yagi et al., 2002; Ouchi et al., 2004) for flat-fielding,
instrumental distortion correction, PSF matching, sky subtraction, and image
stacking. Moreover, implementation details of the 1D reconstruction of Umetsu
et al. (2011) (e.g., the choice of summary statistics) are different from
those of the clumi code that has been tested and calibrated with simulations
(see Umetsu, 2013). Nevertheless, our 1D and 2D results are both in agreement
with the $\Sigma$ profile of Umetsu et al. (2011) in the region of overlap.

##  5 Mass Modeling of A370

In this section, we present mass modeling of A370. With ground-based Subaru
weak-lensing observations alone, we cannot spatially resolve the bimodal
structure of the cluster in the supercritical region (Sections 4.1 and 4.2).
In this study, we thus restrict ourselves to single-component mass models of a
spherical or ellipsoidal halo. We forward model projected cluster lensing
observations by projecting a triaxial or spherical NFW halo model along the
line of sight (e.g., Corless et al., 2009; Sereno & Umetsu, 2011; Umetsu et
al., 2015; Chiu et al., 2018b).

###  5.1 Dynamical State of A370

Molnar et al. (2020) studied the dynamics of A370 using dedicated $N$-body
hydrodynamical simulations of binary cluster mergers constrained by multi-
probe observations. In their simulations, the initial virial masses of two
progenitors were fixed to $1.2\times 10^{15}\,h^{-1}\,M_{\odot}$ and
$1.1\times 10^{15}\,h^{-1}\,M_{\odot}$ with a sum of $2.3\times
10^{15}\,h^{-1}\,M_{\odot}$, according to the result of Umetsu et al. (2011)
who conducted a weak-lensing shear and magnification analysis (see Section
4.3) of five superlens clusters with Suprime-Cam data. In combination with
mass profile constraints from HST strong-lensing data, Umetsu et al. (2011)
obtained $M_{\mathrm{vir}}=2.28^{+0.26}_{-0.22}\times
10^{15}\,h^{-1}\,M_{\odot}$ for A370 assuming a generalized form of the NFW
profile.

Molnar et al. (2020) found that initial conditions of the two progenitors with
an infall velocity of $3500$ km s-1 and an impact parameter of
$70\,h^{-1}\,\mathrm{kpc}$ can reproduce the positions and the offsets between
the peaks of the X-ray emission and the total mass surface density, the
amplitude of the integrated SZE signal (Czakon et al., 2015), and the relative
line-of-sight velocity between the two BCGs ($V\approx 1024$ km s-1).
Moreover, the best-matching simulation reproduces well the velocity dispersion
and the line-of-sight velocity distribution of cluster member galaxies
(Lagattuta et al., 2019; Molnar et al., 2020). These simulation results
support the large total mass of A370 derived from lensing (Umetsu et al.,
2011).

The binary merger simulations of Molnar et al. (2020) suggest that A370 is a
post-major merger of two similar-mass clusters, viewed after the second core
passage in the infalling phase, just before the third core passage. These
results also suggest that the mass distribution of A370 is highly elongated
along the current direction of the collision axis, which is oriented close to
the line of sight in their best simulation, with a viewing angle of
$\vartheta=17.6^{\circ}\pm 3.5^{\circ}$, or $\cos\vartheta=0.95\pm 0.02$.

###  5.2 Triaxial NFW Model

Triaxial modeling of density profiles gives an improved description of
simulated $\Lambda$CDM halos over the conventional spherical model (Jing &
Suto, 2002; Kasun & Evrard, 2005). In this work, we model the cluster mass
distribution with a triaxial NFW density profile. The radial dependence of the
spherical NFW profile is given by (Navarro et al., 1996, 1997)

$\rho(r)=\frac{\rho_{\mathrm{s}}}{(r/r_{\mathrm{s}})(1+r/r_{\mathrm{s}})^{2}},$
(29)

where $\rho_{\mathrm{s}}$ is the scale density and $r_{\mathrm{s}}$ is the
characteristic scale radius at which the logarithmic slope of the density
profile equals $-2$. We generalize the NFW density profile $\rho(r)$ to obtain
its triaxial expression by replacing $r$ and $r_{\mathrm{s}}$ with their
respective ellipsoidal radii $R$ and $R_{\mathrm{s}}$ as

$R^{2}=\frac{X^{2}}{q_{a}^{2}}+\frac{Y^{2}}{q_{b}^{2}}+Z^{2},$ (30)

with $q_{a}$ and $q_{b}$ the minor–major and intermediate–major axis ratios,
respectively. By definition, we have $0<q_{a}\leqslant q_{b}\leqslant 1$.

The degree of triaxiality is defined as (Sereno et al., 2013)

${\cal T}=\frac{1-q_{b}^{2}}{1-q_{a}^{2}},$ (31)

where $0\leqslant{\cal T}\leqslant 1$ by construction. The value of ${\cal T}$
approaches unity at $q_{a}=q_{b}$ (or zero at $q_{b}=1$), if the halo shape is
maximally prolate (or oblate). For $q_{a}=q_{b}=1$, Equation (29) reduces to
the spherical NFW profile $\rho(r)$ with $r=\sqrt{X^{2}+Y^{2}+Z^{2}}$.

We define an ellipsoidal overdensity radius $R_{\Delta}$ (Corless et al.,
2009; Sereno & Umetsu, 2011; Buote & Humphrey, 2012) such that the mean
interior density contained within an ellipsoidal volume of semimajor axis
$R_{\Delta}$ is $\Delta\times\rho_{\mathrm{c}}(z_{l})$. The total mass
enclosed within $R_{\Delta}$ is expressed as

$M_{\Delta}=\frac{4\pi\Delta}{3}q_{a}q_{b}\rho_{\mathrm{c}}(z_{l})R_{\Delta}^{3}.$
(32)

Spherical-equivalent overdensity radii $r_{\Delta}$ are defined by

$r_{\Delta}=(q_{a}q_{b})^{1/3}R_{\Delta}.$ (33)

Similarly, we define $r_{\mathrm{s}}=(q_{a}q_{b})^{1/3}R_{\mathrm{s}}$.

The triaxial concentration parameter is defined as the ratio of the
ellipsoidal overdensity radius $R_{\Delta}$ to the scale radius
$R_{\mathrm{s}}$ along the major axis,

$c_{\Delta}:=\frac{R_{\Delta}}{R_{\mathrm{s}}}=\frac{r_{\Delta}}{r_{\mathrm{s}}}.$
(34)

The characteristic density is then expressed as
$\rho_{\mathrm{s}}=(\Delta\rho_{\mathrm{c}}/3)\times
c_{\Delta}^{3}/[\ln(1+c_{\Delta})-c_{\Delta}/(1+c_{\Delta})]$. In this study,
we use $\Delta=200$ to define the halo mass, $M_{200}$, and the concentration
parameter, $c_{200}$.

A triaxial halo is projected onto the lens plane as elliptical isodensity
contours, which can be expressed as a function of the intrinsic halo axis
ratios ($q_{a},q_{b}$) and orientation angles ($\vartheta,\phi,\psi$) with
respect to the observer’s line of sight. Following Umetsu et al. (2015), we
adopt the $z$-$x$-$z$ convention of Euler angles ($\vartheta,\phi,\psi$) to be
consistent with Stark (1977, see also ). The angle $\vartheta$ represents the
inclination of the major axis ($Z$) with respect to the line of sight.

After a rotation by the first two Euler angles ($\vartheta,\phi$), elliptical
isodensity contours of the projected ellipsoid can be described as a function
of the elliptical radius $\zeta$, expressed in terms of projected Cartesian
coordinates $(x^{\prime},y^{\prime})$ as

$\zeta^{2}=\frac{1}{f}\left(jx^{\prime 2}+2kx^{\prime}y^{\prime}+ly^{\prime
2}\right),$ (35)

where

$\displaystyle j$
$\displaystyle=\cos^{2}\vartheta\left(\frac{1}{q_{a}^{2}}\cos^{2}\phi+\frac{1}{q_{b}^{2}}\sin^{2}\phi\right)+\frac{1}{q_{a}^{2}q_{b}^{2}}\sin^{2}\vartheta,$
(36) $\displaystyle k$
$\displaystyle=\sin\phi\cos\phi\cos\vartheta\left(\frac{1}{q_{a}^{2}}-\frac{1}{q_{b}^{2}}\right),$
$\displaystyle l$
$\displaystyle=\frac{1}{q_{a}^{2}}\sin^{2}\phi+\frac{1}{q_{b}^{2}}\cos^{2}\phi,$
$\displaystyle f$
$\displaystyle=\sin^{2}\vartheta\left(\frac{1}{q_{a}^{2}}\sin^{2}\phi+\frac{1}{q_{b}^{2}}\cos^{2}\phi\right)+\cos^{2}\vartheta.$

The minor–major axis ratio ($\leqslant 1$) of the elliptical isodensities is
expressed as
$q_{\perp}(q_{a},q_{b},\vartheta,\phi)=\sqrt{\frac{j+l-\sqrt{(j-l)^{2}+4k^{2}}}{j+l+\sqrt{(j-l)^{2}+4k^{2}}}}.$
Finally, the third Euler angle $\psi$ describes the additional rotational
degree of freedom in the sky plane to specify the observer’s coordinate system
$(x,y)$, defined such that $x^{\prime}=x\cos\psi-y\sin\psi$ and
$y^{\prime}=x\sin\psi+y\cos\psi$.

For a self-similar mass model expressed as
$\rho(R)=\rho_{\mathrm{s}}f_{\mathrm{3D}}(R/R_{\mathrm{s}})$, the projected
mass density $\Sigma(\zeta)$ is related to $\rho(R)$ (see Equation (29)) as
(Umetsu et al., 2015)

$\Sigma(\zeta)=\frac{2R_{\mathrm{s}}\rho_{\mathrm{s}}}{\sqrt{f}}\int_{\zeta/R_{\mathrm{s}}}^{\infty}\frac{f_{\mathrm{3D}}(u)udu}{\sqrt{u^{2}-(\xi/\xi_{\mathrm{s}})^{2}}}\equiv\Sigma_{\mathrm{s}}f_{\mathrm{2D}}(\xi/\xi_{\mathrm{s}}),$
(37)

where $\Sigma_{\mathrm{s}}$ is the scale surface mass density defined by

$\Sigma_{\mathrm{s}}=2R_{\mathrm{s}}\rho_{\mathrm{s}}/\sqrt{f},$ (38)

$\xi=\sqrt{x^{\prime\prime 2}+y^{\prime\prime 2}/q_{\perp}^{2}}$, and
$\xi_{\mathrm{s}}$ is the semi-major scale length of the projected halo. Here
we have chosen the new coordinate system $(x^{\prime\prime},y^{\prime\prime})$
such that the $x^{\prime\prime}$ axis is aligned with the major axis of the
projected ellipse. In this study, we employ the radial dependence of the
projected NFW profile $f_{\mathrm{2D}}(u)$ as given by Wright & Brainerd
(2000).

To summarize, our mass model is specified by a total of seven parameters
describing the total matter ellipsoid, namely, halo mass and concentration
($M_{200},c_{200}$), intrinsic axis ratios ($q_{a},q_{b}$), and three Euler
angles ($\vartheta,\phi,\psi$):

$\mbox{\boldmath$p$}=\\{M_{200},c_{200},q_{a},q_{b},\vartheta,\phi,\psi\\}.$
(39)

In this way, for a given set of the model parameters, we can project a
triaxial (or spherical) NFW halo onto the lens plane and compute the surface
mass density $\Sigma(x,y)$ at each angular position. As discussed in Umetsu et
al. (2015) (see also Sereno & Umetsu, 2011), however, it should be noted that
2D lensing observations can effectively constrain only four observationally
accessible parameters, namely, $\Sigma_{\mathrm{s}}$, $\xi_{\mathrm{s}}$,
$q_{\perp}$, and the position angle of the projected major axis (Gavazzi,
2005). That is, the deprojection of triaxial systems is intrinsically
underconstrained (Limousin et al., 2013). On the other hand, the spherical NFW
model ($q_{a}=q_{b}=1$) is specified by two parameters, $(M_{200},c_{200})$,
which can be constrained by data in principle.

###  5.3 Bayesian Inference Procedure

The likelihood function ${\cal L}$ of the 2D mass distribution data
$\mbox{\boldmath$m$}=\\{\Sigma(\mbox{\boldmath$\theta$}_{n})\\}_{n=1}^{N_{\mathrm{pix}}}$
given a set of model parameters $p$ is expressed as (Oguri et al., 2005;
Umetsu et al., 2015)

$\displaystyle-2\ln{\cal L}(\mbox{\boldmath$p$})=$
$\displaystyle\sum_{m,n=1}^{N_{\mathrm{pix}}}\left[\Sigma-\widehat{\Sigma}(\mbox{\boldmath$p$})\right]_{m}\left(C^{-1}\right)_{mn}\left[\Sigma-\widehat{\Sigma}(\mbox{\boldmath$p$})\right]_{n}$
(40) $\displaystyle+\ln\left[(2\pi)^{N_{\mathrm{pix}}}\mathrm{det(C)}\right]$

where
$\widehat{\Sigma}_{n}(\mbox{\boldmath$p$})=\Sigma(\mbox{\boldmath$\theta$}_{n}|\mbox{\boldmath$p$})$
is the surface mass density at the grid position
$\mbox{\boldmath$\theta$}_{n}$ predicted by the model $p$ and
$C=C_{\mathrm{stat}}+C_{\mathrm{lss}}$ is the total covariance matrix
(Equation (26)).

We use a Bayesian Markov Chain Monte Carlo (MCMC) algorithm to obtain a well-
characterized inference of the model $p$. We consider the following three
different modeling approaches: (1) spherical modeling with uninformative
uniform priors on $\log{M_{200}}$ and $\log{c_{200}}$, (2) fiducial triaxial
modeling with uninformative uniform priors on all parameters, and (3) triaxial
modeling incorporating an informative line-of-sight (LOS) prior from Molnar et
al. (2020). For simplicity, we refer to these three modeling approaches as
Spherical, Triaxial, and Triaxial+LOS modeling, respectively.

Here we briefly summarize the assumed priors for each case.

  1. 1.

Spherical modeling: We float only two parameters ($M_{200},c_{200}$) and fix
the remaining parameters ($q_{a}=q_{b}=1$ and $\vartheta=\phi=\psi=0$). We
employ uninformative log-uniform priors for $M_{200}$ and $c_{200}$ in the
range $M_{200}/(\,h^{-1}\,M_{\odot})\in[10^{14},10^{16}]$ and
$c_{200}\in[1,10]$.999It is appropriate to assume a log-uniform prior, instead
of a uniform prior, for a positive-definite quantity such as $M_{200}$ and
$c_{200}$, especially if the quantity spans a wide dynamic range (Umetsu,
2020). Since the corresponding prior distributions in $M_{200}$ and $c_{200}$
scale as $1/M_{200}$ and $1/c_{200}$, the choice of their lower bounds is
relatively important.

  2. 2.

Fiducial triaxial modeling: We use uniform priors on the intrinsic shapes
($q_{a},q_{b}$) and orientation angles ($\cos\vartheta,\phi,\psi$), while
keeping the same log-uniform priors on $M_{200}$ and $c_{200}$ as in the
spherical case. We assume the following form of the prior PDF for the
intrinsic axis ratios:

$P(q_{a},q_{b})=P(q_{b}|q_{a})P(q_{a}),$ (41)

where

$\displaystyle
P(q_{a})=\begin{cases}1/(1-q_{\mathrm{min}})&\mathrm{for}~{}q_{\mathrm{min}}<q_{a}\leqslant
1\\\ 0&\mathrm{for}~{}q_{a}\leqslant q_{\mathrm{min}}\end{cases},$ (42)
$\displaystyle
P(q_{b}|q_{a})=\begin{cases}1/(1-q_{a})&\mathrm{for}~{}q_{b}\geqslant q_{a}\\\
0&\mathrm{for}~{}q_{b}<q_{a}\end{cases},$

and $q_{\mathrm{min}}=0.1$ is the lower bound of the minor-to-major axis ratio
$q_{a}$ (e.g., Oguri et al., 2005; Chiu et al., 2018b), which is introduced to
exclude unstable configurations that are not expected for cluster halos. For
the orientation angles, we consider a population of randomly oriented halos
with $P(\cos\vartheta)=1$ for $0\leqslant\cos\vartheta\leqslant 1$,
$P(\phi)=1/\pi$ for $-\pi/2\leqslant\phi\leqslant\pi/2$, and $P(\psi)=1/\pi$
for $-\pi/2\leqslant\psi\leqslant\pi/2$.

  3. 3.

Triaxial+LOS modeling: We adopt an informative prior on $\cos\vartheta$ based
on the binary merger simulations of Molnar et al. (2020). For the other
parameters, we use the same priors as for the fiducial triaxial modeling.
Specifically, we employ a Gaussian prior on $\cos\vartheta$ of $0.95\pm 0.02$
(Section 5.1) truncated in the range $0\leqslant\cos\vartheta\leqslant 1$.

For comparison purposes, we perform spherical NFW modeling with the surface
mass density profile
$\mbox{\boldmath$m_{\mathrm{1D}}$}=\\{\Sigma_{\mathrm{min}},\Sigma_{i}\\}_{i=1}^{N_{\mathrm{bin}}}$
derived from the WL-1D analysis (Section 4.1). The likelihood function ${\cal
L}(\mbox{\boldmath$p$})$ for the WL-1D analysis is defined as in Equation (26)
of Umetsu et al. (2014). In the covariance matrix of WL-1D, we account for
systematic effects due to the residual mass-sheet degeneracy, in addition to
the measurement error and cosmic noise contributions (see Section 2.7). This
residual uncertainty is estimated in each $\Sigma$ bin as a difference between
the joint and marginal posterior solutions (see Umetsu et al.,
2014).101010Because of the large number of parameters involved, we do not
explore the whole likelihood surface in the clumi-2D code, and thus we are not
able to include the systematic term in the WL-2D analysis. As we have seen in
Figure 10, our WL-1D and WL-2D results are consistent with each other, with no
significant evidence for a systematic offset.

Similarly, we also perform spherical NFW modeling with the reduced tangential
shear profile $\\{g_{+,i}\\}_{i=1}^{N_{\mathrm{bin}}}$ obtained in our WL-1D
analysis, because this tangential shear fitting is the standard approach to
infer cluster masses from weak-lensing data (e.g., Okabe et al., 2013;
Applegate et al., 2014; Hoekstra et al., 2015; Schrabback et al., 2018). Here
we account for the measurement error and cosmic noise contributions in the
covariance matrix (see Section 4.4 of Umetsu, 2020).

###  5.4 Posterior Parameter Constraints

Figure 11: Marginalized 1D and 2D (68$\%$, 95$\%$, and 99.7$\%$ confidence
level contours) posterior distributions for the NFW model parameters
($\log_{10}{M_{\mathrm{200}}},\log_{10}{c_{\mathrm{200}}}$) obtained using
uniform priors assuming spherical symmetry. Blue shaded contours show the
constraints obtained from the 1D mass reconstruction (WL-1D; see Figure 7).
Orange contours show the constraints from the 2D reconstruction (WL-2D; see
Figure 9). Figure 12: Marginalized 1D and 2D (68$\%$, 95$\%$, and 99.7$\%$
confidence level contours) posterior distributions for the triaxial NFW
parameters of interest
$(\log_{10}{M_{200}},\log_{10}{c_{200}},q_{a},q_{b},\cos\vartheta)$. Orange
and blue shaded contours show the results obtained with and without the line-
of-sight prior on $\cos\vartheta$ based on binary merger simulations of Molnar
et al. (2020). For each parameter, the dashed line shows the biweight central
location ($C_{\mathrm{BI}}$) of the marginalized 1D distribution. Table 6:
Posterior summary of the NFW halo parameters for A370

Method | Data | Prior | $M_{200}$ | $c_{200}$ | $q_{a}$ | $q_{b}$ | $\cos\vartheta$ | ${\cal T}$
---|---|---|---|---|---|---|---|---
|  |  | $(10^{15}\,h^{-1}\,M_{\odot})$ |  |  |  |  |
Tangential shear |  $g_{+}$ profile | Spherical | $1.55\pm 0.26$ | $5.61\pm 1.68$ | — | — | — |
WL-1D |  $\kappa$ profile | Spherical | $1.72\pm 0.32$ | $4.90\pm 1.58$ | — | — | — |
WL-2D |  $\kappa$ map | Spherical | $1.51\pm 0.22$ | $4.91\pm 0.93$ | — | — | — | —
WL-2D |  $\kappa$ map | Triaxial (fiducial) | $1.54\pm 0.29$ | $5.27\pm 1.28$ | $0.62\pm 0.23$ | $0.86\pm 0.13$ | $0.58\pm 0.29$ | $0.51\pm 0.29$
WL-2D |  $\kappa$ map | Triaxial+LOS | $1.38\pm 0.20$ | $4.45\pm 0.93$ | $0.59\pm 0.23$ | $0.80\pm 0.21$ | $0.95\pm 0.02$ | $0.58\pm 0.30$

Note. — Cluster halo parameters derived from a spherical or triaxial NFW fit
to Subaru weak-lensing data. We adopt a concordance cosmology of $h=0.7$,
$\Omega_{\mathrm{m}}=0.3$, and $\Omega_{\Lambda}=0.7$. We note that the degree
of triaxiality ${\cal T}$ is a derived parameter that depends on $q_{a}$ and
$q_{b}$ (see Equation (31)). As posterior summary statistics, we use the
biweight estimator of Beers et al. (1990) to represent the center location
($C_{\mathrm{BI}}$) and the spread ($S_{\mathrm{BI}}$) of marginalized 1D
posterior distributions. For each parameter, symmetrized biweight statistics
$C_{\mathrm{BI}}\pm S_{\mathrm{BI}}$ are shown. The $\kappa(\theta)$ profile
is reconstructed from the WL-1D analysis of the
$\\{g_{+}(\theta),n_{\mu}(\theta)\\}$ data set (Section 4.1), while the
$\kappa(\mbox{\boldmath$\theta$})$ map from the WL-2D analysis of the
$\\{g_{1}(\mbox{\boldmath$\theta$}),g_{2}(\mbox{\boldmath$\theta$}),n_{\mu}(\theta)\\}$
data set (Section 4.2).

Table 7: Weak-lensing mass estimates for A370 Method | Data | Prior | $M_{2500}$ | $M_{1000}$ | $M_{500}$ | $M_{\mathrm{vir}}$ | $M_{\mathrm{200m}}$ |  $M_{\mathrm{s}}$aaTotal mass enclosed within the NFW scale radius $r_{\mathrm{s}}$, $M(<r_{\mathrm{s}})$
---|---|---|---|---|---|---|---|---
|  |  | $(10^{15}\,h^{-1}\,M_{\odot})$ | $(10^{15}\,h^{-1}\,M_{\odot})$ | $(10^{15}\,h^{-1}\,M_{\odot})$ | $(10^{15}h^{-1}\,h^{-1}\,M_{\odot})$ | $(10^{15}\,h^{-1}\,M_{\odot})$ | $(10^{15}\,h^{-1}\,M_{\odot})$
Tangential shear |  $g_{+}$ profile | Spherical | $0.53\pm 0.07$ | $0.85\pm 0.10$ | $1.13\pm 0.15$ | $1.76\pm 0.32$ | $1.86\pm 0.35$ | $0.29\pm 0.11$
WL-1D |  $\kappa$ profile | Spherical | $0.54\pm 0.08$ | $0.90\pm 0.11$ | $1.23\pm 0.18$ | $1.97\pm 0.40$ | $2.09\pm 0.44$ | $0.36\pm 0.14$
WL-2D |  $\kappa$ map | Spherical | $0.48\pm 0.05$ | $0.80\pm 0.09$ | $1.09\pm 0.13$ | $1.73\pm 0.27$ | $1.84\pm 0.29$ | $0.31\pm 0.08$
WL-2D |  $\kappa$ map | Triaxial (fiducial) | $0.50\pm 0.11$ | $0.82\pm 0.16$ | $1.11\pm 0.20$ | $1.76\pm 0.33$ | $1.86\pm 0.35$ | $0.31\pm 0.08$
WL-2D |  $\kappa$ map | Triaxial+LOS | $0.42\pm 0.07$ | $0.71\pm 0.09$ | $0.97\pm 0.13$ | $1.58\pm 0.24$ | $1.68\pm 0.26$ | $0.31\pm 0.08$

The main results from our Bayesian inference of the spherical and triaxial NFW
models are summarized in Table 6. As summary statistics, we employ the
biweight estimator of Beers et al. (1990) to represent the center location
($C_{\mathrm{BI}}$) and the scale or spread ($S_{\mathrm{BI}}$) of
marginalized 1D posterior PDFs (e.g., Umetsu et al., 2020). For a lognormally
distributed quantity, $C_{\mathrm{BI}}$ approximates the median of the
distribution.

Triaxial modeling allows for a more general description of the intrinsic shape
of cluster halos, leading to broader posterior distributions than the
spherical case (Oguri et al., 2005; Sereno & Umetsu, 2011). The parameter
constraints become more degenerate because of the lack of information of the
halo elongation along the line of sight. These trends are found in the
posterior distributions from our data.

Our spherical modeling of the WL-2D data yields $M_{200}=(1.51\pm 0.22)\times
10^{15}\,h^{-1}\,M_{\odot}$ (or $r_{200}=(1.64\pm
0.08)\,h^{-1}\,\mathrm{Mpc}$) and $c_{200}=4.91\pm 0.93$, which are consistent
with results from the tangential-shear and WL-1D methods. From triaxial
modeling, we obtain $M_{200}=(1.54\pm 0.29)\times 10^{15}\,h^{-1}\,M_{\odot}$
(or $r_{200}=(1.65\pm 0.10)\,h^{-1}\,\mathrm{Mpc}$), $c_{200}=5.27\pm 1.28$,
and a degree of triaxiality ${\cal T}=0.51\pm 0.29$ with the fiducial priors
and $M_{200}=(1.38\pm 0.20)\times 10^{15}\,h^{-1}\,M_{\odot}$ (or
$r_{200}=(1.59\pm 0.08)\,h^{-1}\,\mathrm{Mpc}$), $c_{200}=4.45\pm 0.93$, and
${\cal T}=0.58\pm 0.30$ when the LOS prior is employed. Thus, the level of
mass bias due to the LOS elongation is found to be $\sim 10\%$ for A370.

We now turn to the impact of the prior on the inferred uncertainty in the mass
determination. Spherical (or Triaxial+LOS) modeling of A370 yields a
fractional uncertainty in $M_{200}$ of
$S_{\mathrm{BI}}(M_{200})/C_{\mathrm{BI}}(M_{200})\approx 15\%$. In contrast,
the uncertainty in $M_{200}$ from Triaxial modeling is $\approx 19\%$, which
is a factor of $\sim 1.3$ larger than that of Spherical modeling. It is
insightful to compare our results with those of Umetsu et al. (2015), who
performed a WL-2D analysis of the superlens cluster Abell 1689 based on deeper
Suprime-Cam observations. Analyzing their WL-2D data, Umetsu et al. (2015)
obtained fractional uncertainties in $M_{200}$ of $\approx 8\%$ and $20\%$ for
their spherical and full-triaxial NFW models, respectively (see the first and
second rows of Table 7 in Umetsu et al. 2015). For both clusters, the
fractional mass uncertainty in full triaxial modeling is $\sim 20\%$,
suggesting that the mass accuracy in deep weak-lensing observations is
essentially limited by the uncertainty in the intrinsic shape and orientation
of the cluster. Similar trends are also found for the concentration parameter.
To accurately infer the cluster mass and concentration from lensing, it is
thus necessary to directly model or marginalize over the 3D shape of clusters;
when spherical symmetry is assumed, the effect of the intrinsic shape of the
cluster should be accounted for in the error analysis (e.g., Gruen et al.,
2015; Umetsu et al., 2016).

We also derive summary statistics on the total mass $M_{\Delta}$ evaluated at
several characteristic interior overdensities $\Delta$. Table 7 lists the
results of our cluster mass estimates. Our estimates of $M_{\mathrm{vir}}$
obtained without the LOS information are consistent within the errors with
$M_{\mathrm{vir}}=2.28^{+0.26}_{-0.22}\times 10^{15}\,h^{-1}\,M_{\odot}$ from
the combined weak- and strong-lensing analysis of Umetsu et al. (2011, see
Section 5.1). In particular, our WL-1D analysis yields
$M_{\mathrm{vir}}=(1.97\pm 0.40)\times 10^{15}\,h^{-1}\,M_{\odot}$ and
$c_{\mathrm{vir}}=5.91\pm 1.87$, in agreement with the results of Umetsu et
al. (2011). We find that spherical mass estimates from the WL-2D analysis are
slightly lower than but consistent within the errors with the WL-1D results.

It should be noted that the halo mass $M_{\Delta}$ constrained using the LOS
prior is likely to be considerably lower than the sum of the initial bound
masses of the two progenitors, because A370 is expected to be in a highly
disturbed dynamical state (see Section 5.1). Our estimates of
$M_{\mathrm{vir}}$ obtained without the LOS prior are consistent to better
than $2\sigma$ with the total mass of the system $M_{\mathrm{vir}}=2.3\times
10^{15}\,h^{-1}\,M_{\odot}$ adopted in the binary merger simulations of Molnar
et al. (2020).

Overall, our results agree well with weak-lensing mass estimates of Hoekstra
et al. (2015), who obtained $M_{\mathrm{vir}}=2.13^{+0.39}_{-0.37}\times
10^{15}\,h^{-1}\,M_{\odot}$ and $M_{500}=(1.23\pm 0.22)\times
10^{15}\,h^{-1}\,M_{\odot}$ for A370. Our estimates of $M_{500}$ are also in
agreement with the recent caustic mass estimate by Lagattuta et al. (2022),
$M_{500}=(0.92\pm 0.11)\times 10^{15}\,h^{-1}\,M_{\odot}$, obtained from a
detailed phase-space analysis of HST BUFFALO imaging and VLT Multi-Unit
Spectroscopic Explorer (MUSE) spectroscopic observations.

In Figure 11, we show posterior constraints on the NFW parameters
$(\log{M_{200}},\log{c_{200}})$ inferred from spherical modeling of both WL-1D
and WL-2D data. The blue and orange contours in the lower-left panel represent
the joint posterior PDFs for WL-1D and WL-2D, respectively, showing good
agreement between the two methods. In both cases, the marginalized posterior
PDFs for $(\log{M_{200}},\log{c_{200}})$ are unimodal and symmetric.

Figure 12 displays the marginalized posterior PDFs for the triaxial NFW
parameters of interest
($\log{M_{200}},\log{c_{200}},q_{a},q_{b},\cos{\vartheta}$) based on our WL-2D
analysis. The results with and without employing the LOS prior are compared in
the figure. The posterior PDFs for $\log{M_{200}}$ and $\log{c_{200}}$ are
clearly unimodal and fairly symmetric. For both parameters, there is no
significant shift in the PDF with respect to the spherical case (Figure 11).
The posterior PDFs for the shape and orientation parameters
($q_{a},q_{b},\cos{\vartheta}$) from our fiducial modeling are very broad,
reflecting the fact that the deprojection of triaxial halos is intrinsically
underconstrained. In contrast, the axis ratio of the projected mass
distribution, $q_{\perp}(q_{a},q_{b},\vartheta,\phi)$ (Section 5.2), can be
directly constrained by the WL-2D data. Our posterior inference of the
projected axis ratio is $q_{\perp}=0.78\pm 0.13$ and $0.79\pm 0.13$ with and
without using the LOS prior, respectively. This is slightly larger than, but
consistent with, the median axis ratio $\overline{q}_{\perp}\sim 0.6$ expected
for randomly oriented cluster-scale CDM halos (Umetsu et al., 2018, see also
Bonamigo et al. 2015; Suto et al. 2016).

Compared to the fiducial results obtained with uniform priors, our inference
with the informative Gaussian prior on $\cos{\vartheta}$ (Triaxial+LOS)
prefers a more prolate geometry with lower mass and lower concentration. In
fact, there is a slight increase in the posterior probability for a prolate
configuration ($q_{a}\sim q_{b}\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}0.6$) with lower mass and lower
concentration. This can be understood as a consequence of the boosted surface
mass density of the cluster lens due to the strongly aligned configuration.

##  6 X-ray Data and Analysis

Figure 13: Logarithmically scaled, exposure-corrected, and background-
subtracted Chandra ACIS image of A370 in the $0.5$–$7$ keV band, smoothed with
a Gaussian of $\mathrm{FWHM}=4.6\arcsec$. The image is $5.2\arcmin\times
5.2\arcmin$ ($1.1\,h^{-1}\,\mathrm{Mpc}$ on a side at $z=0.375$) in size and
centered on the optical cluster center (open diamond symbol). The positions of
the two BCGs are marked with $+$ symbols. The X-ray centroid position is
marked with a $\times$ symbol. A bright X-ray source in the north
corresponding to a foreground elliptical galaxy is masked by a black solid
circle. North is up and east is to the left.

Here we describe our analysis of archival Chandra X-ray data (Section 6.1). We
use two complementary approaches to determining the 3D gas density and
temperature profiles of A370 under the assumption of spherical symmetry.
First, we derive gas densities and temperatures of the cluster in concentric
spherical shells from a spectral deprojection analysis (Section 6.2). Second,
we perform forward modeling to simultaneously fit X-ray surface brightness
profiles binned in multiple energy bands to infer the 3D gas density and
temperature profiles in a parametric form (Section 6.3). With the forward-
fitting method, we will also derive the total mass profile assuming
hydrostatic equilibrium.

###  6.1 Chandra Data Reduction

We analyze archival X-ray data of A370 taken with the Advanced CCD Imaging
Spectrometer (ACIS; Garmire et al., 2003) on board the Chandra X-ray
Observatory. The observation identification (ObsID) numbers of Chandra
observations analyzed in this study are 515 and 7715. Our analysis uses the
Chandra Interactive Analysis of Observations software (ciao, version 4.13;
Fruscione et al., 2006) and the Chandra Calibration Database (caldb, version
4.9.5). We checked the light curve of each data set using the lc_clean task in
ciao, filtering flare data. The net exposure time of each data set is 62.9 ks
and 7.1 ks. Point sources were identified using the wavdetect task in ciao and
excluded from the analysis.

In our spectral analysis, we use the X-ray Spectral Fitting Package (Arnaud,
1996, xspec version 12.11.1) and the atomdb code (Heuer et al., 2021, version
3.0.9) for plasma emission modeling, assuming that the ICM is in collisional
ionization equilibrium (Smith et al., 2001). The abundance table of Anders &
Grevesse (1989) is used in xspec. Here, the abundance of a given element is
defined as
$Z_{i}=(n_{i,\mathrm{obs}}/n_{\mathrm{H,obs}})/(n_{i,\odot}/n_{\mathrm{H,\odot}})$,
where $n_{i}$ and $n_{\mathrm{H}}$ are the number densities of the $i$th
element and hydrogen, respectively. We use the iron abundance to represent the
ICM metal abundance, such that the abundance of other elements is tied to the
iron abundance as $Z_{i}=Z_{\mathrm{Fe}}$ (Ueda et al., 2021). The Galactic
absorption column density is estimated at $N_{\mathrm{H}}=2.89\times 10^{20}$
cm-2 according to HI4PI Collaboration et al. (2016) and fixed in our X-ray
spectral analysis. The blank-sky data included in caldb are used to determine
the background contribution.

To determine the centroid of X-ray emission in A370, we fit the surface
brightness distribution with a 2D $\beta$-model using the sherpa fitting
package in ciao (Freeman et al., 2001; Doe et al., 2007; Burke et al., 2021).
The surface brightness map was extracted from the ACIS S3 chip in the data set
of ObsID 515 to reduce the uncertainty in the background determination. A
bright foreground galaxy lying about $2\arcmin$ north of the cluster center
was masked with a circle of radius $30\arcsec$ from its X-ray peak. From the
best-fit model, we find the X-ray centroid of
$\mathrm{R.A.}=\mathrm{2:39:53.2}$ and $\mathrm{decl.}=-\mathrm{1:34:35.1}$
(Table 1), with positional uncertainties of
$(\Delta\mathrm{R.A.},\Delta\mathrm{decl.})=(0.33\arcsec,0.39\arcsec)$.

Figure 13 shows the exposure-corrected and background-subtracted Chandra ACIS
image of A370 in the $0.5$–$7$ keV band, smoothed with $4.6\arcsec$ FWHM
Gaussian. The X-ray emission centroid determined from Chandra observations is
$4.9\arcsec$ ($\approx 18\,h^{-1}\,\mathrm{kpc}$) away from the optical center
defined as the midpoint of the two BCGs.

###  6.2 Spectral Deprojection Analysis

Figure 14:  Three-dimensional gas temperature $T_{\mathrm{gas}}(r)$ (upper panel) and enclosed gas mass $M_{\mathrm{gas}}(<r)$ (lower panel) profiles of A370 derived from Chandra X-ray observations. The red shaded region in each panel shows the marginalized $1\sigma$ confidence region of the respective profile obtained from forward modeling of Chandra X-ray data. Blue open boxes show the gas temperatures with $1\sigma$ uncertainties from spectral deprojection in concentric spherical shells. Blue filled circles with error bars in the lower panel show the results obtained from spectral deprojection.  Table 8: Best-fit parameters from the Chandra spectral deprojection analysis $\theta_{\mathrm{in}}$ | $\theta_{\mathrm{out}}$ | $T_{\mathrm{gas}}$ | $Z$ | $n_{\mathrm{e}}$
---|---|---|---|---
($\arcsec$) | ($\arcsec$) | (keV) | ($Z_{\odot}$) | ($10^{-3}$ cm-3)
$0$ | $25$ | $10.16_{-2.06}^{+3.11}$ | $0.58_{-0.41}^{+0.45}$ | $4.754\pm 0.191$
$25$ | $50$ | $9.77_{-0.97}^{+1.20}$ | $0.71_{-0.21}^{+0.23}$ | $3.241\pm 0.067$
$50$ | $100$ | $8.40_{-0.61}^{+0.77}$ | $0.35\pm 0.11$ | $1.595\pm 0.020$
$100$ | $200$ | $9.29_{-1.25}^{+1.64}$ |  $0.3$ (fixed) | $0.502\pm 0.005$

We jointly fit background-subtracted Chandra spectra in the $0.4$–$7.0$ keV
band extracted from four concentric annular regions
($\theta_{\mathrm{in}},\theta_{\mathrm{out}}$) of $0\arcsec$–$25\arcsec$,
$25\arcsec$–$50\arcsec$, $50\arcsec$–$100\arcsec$, and
$100\arcsec$–$200\arcsec$ centered on the X-ray centroid using the projct
routine implemented in xspec (Smith et al., 2002). Here projct allows us to
fit spectra extracted from a series of concentric annuli simultaneously,
assuming spherical symmetry to calculate suitable geometric factors (e.g.,
Fabian et al., 1980, 1981; Kriss et al., 1983; Arabadjis et al., 2002). In
this analysis, the cluster redshift and $N_{\mathrm{H}}$ are fixed at $0.375$
and $2.89\times 10^{20}$ cm-2, respectively. The metal abundance of the ICM in
the $100\arcsec$–$200\arcsec$ region is assumed to be $0.3Z_{\odot}$ (Fujita
et al., 2008; Werner et al., 2013; Urban et al., 2017; Ghizzardi et al.,
2021).

To set the outer boundary conditions, we fit the background-subtracted X-ray
spectrum in the $0.4$–$2.0$ keV band extracted from the outermost annular
region of $200\arcsec$–$400\arcsec$, ignoring the emission from gas outside
the outermost shell (see Humphrey et al., 2006) and fixing the metal abundance
to $0.3Z_{\odot}$. The best-fit parameters for the outermost shell,
$T_{\mathrm{3D}}=9.7^{+10.2}_{-3.3}$ keV and $n_{\mathrm{e}}=(1.76\pm
0.08)\times 10^{-4}$ cm-3, are included and fixed in our deprojection analysis
of the inner concentric regions.

The resulting best-fit parameters for each spherical shell are summarized in
Table 8 (see also Figure 14). The gas density $\rho_{\mathrm{gas}}$ is related
to the electron number density $n_{\mathrm{e}}$ as
$\rho_{\mathrm{gas}}=\mu_{\mathrm{e}}m_{\mathrm{p}}n_{\mathrm{e}}$, with
$\mu_{\mathrm{e}}\approx 1.11$ the mean mass per electron in units of proton
mass $m_{\mathrm{p}}$.

###  6.3 Parametric Forward Fitting

We perform a forward model fitting of the Chandra observations for A370. The
mbproj2 algorithm developed by Sanders (2017, see also ) is capable of
modeling radial X-ray surface brightness profiles in multiple energy bins,
with or without assuming hydrostatic equilibrium. Motivated by their work, we
have implemented a forward-modeling algorithm to simultaneously fit the X-ray
brightness profiles binned in multiple energy bands to infer the 3D gas
density and temperature profiles in a parametric form, without assuming
hydrostatic equilibrium. Both algorithms assume spherical symmetry.

We model the 3D gas density profile $n_{\mathrm{e}}(r)$ as a $\beta$-profile
and the 3D temperature profile $T_{\mathrm{3D}}(r)$ as a universal temperature
profile of Vikhlinin et al. (2006):

$\displaystyle n_{\mathrm{e}}(r)$
$\displaystyle=n_{\mathrm{e0}}\left[1+(r/r_{\mathrm{c}})^{2}\right]^{-3\beta/2},$
(43) $\displaystyle T_{\mathrm{3D}}(r)$
$\displaystyle=T_{0}\frac{(r/r_{\mathrm{t}})^{a}}{\left[1+(r/r_{\mathrm{t}})^{2}\right]^{c/2}},$

where $n_{\mathrm{e0}}$ is the central electron number density, $\beta$ is the
slope parameter, and $r_{\mathrm{c}}$ is the core radius of the $\beta$
profile; $T_{0}$ is the central gas temperature, $a$ and $c$ are the
temperature slope parameters, and $r_{\mathrm{t}}$ is the temperature scale
radius.

In this analysis, we fix the cluster redshift to $z=0.375$, the Galactic
absorption column density to $N_{\mathrm{H}}=2.89\times 10^{20}$ cm-2, and the
metal abundance of the ICM to $Z=0.3Z_{\odot}$. We use the spectroscopic-like
temperature $T_{\mathrm{2D}}(r)$ of Mazzotta et al. (2004) to approximate
spectroscopic temperatures extracted from Chandra X-ray observations:

$T_{\mathrm{2D}}=\frac{\int\\!wT_{\mathrm{3D}}dV}{\int\\!wdV},$ (44)

with $w=n_{\mathrm{e}}^{2}(r)T_{\mathrm{3D}}^{-3/4}(r)$. The X-ray surface
brightness $S_{\mathrm{X}}(r_{\perp})$ as a function of projected cluster
radius $r_{\perp}=D_{l}\theta$ is modeled by the following equation (Ettori,
2000):

$S_{\mathrm{X}}(r_{\perp})=n_{\mathrm{e0}}n_{\mathrm{p0}}r_{\mathrm{c}}\Lambda_{\mathrm{X}}B\left(3\beta-\frac{1}{2},\frac{1}{2}\right)\left[1+\left(\frac{r_{\perp}}{r_{\mathrm{c}}}\right)^{2}\right]^{\frac{1}{2}-3\beta},$
(45)

where $\Lambda_{\mathrm{X}}(T_{\mathrm{2D}},Z)$ is the cooling function,
$n_{\mathrm{p0}}\approx n_{\mathrm{e0}}/1.17$ is the central proton number
density, and $B(x,y)$ is the beta function. We use the pyatomdb python package
(Foster & Heuer, 2020) to evaluate the cooling function
$\Lambda_{\mathrm{X}}(T_{\mathrm{2D}},Z)$ in each energy band for a given
value of the spectroscopic-like temperature $T_{\mathrm{2D}}(r_{\perp})$.

We have extracted the radial profiles of X-ray surface brightness in
$N_{\mathrm{spec}}=10$ energy bands between neighbouring energies of 0.5,
0.75, 1, 1.25, 1.5, 2, 3, 4, 5, 6, and 7 keV. In each energy band, the X-ray
surface brightness is sampled in $49$ linearly spaced radial bins in the range
$\theta\in[10\arcsec,200\arcsec]$ centered on the X-ray centroid. Following
Sanders et al. (2018), we have chosen these bands so as to capture most of the
spectral information without overly increasing the computational time. We
estimate in each radial bin the pixel-to-pixel variance of X-ray brightness
over the $0.4$–$7$ keV energy band, finding that the standard errors of the
mean based on the estimated variance are highly consistent with the errors
determined based on the photon counts. In this work, we use the standard error
based on the estimated variance to characterize the uncertainty in the mean
X-ray surface brightness in each bin.

The background contribution in each energy band is determined from the blank-
sky data included in caldb (Section 6.1). We estimate the count rate of the
blank-sky data in the spectral range of $9$–$12$ keV dominated by the particle
background (Hickox & Markevitch, 2006). Using the ratio between the count rate
observed in A370 and the background one in the $9$–$12$ keV band, we rescale
the background contribution in each energy band to match the observations of
A370, accounting for the difference in exposure times. We then construct the
azimuthally averaged radial profile of the background map in each energy band.
Similarly, we create azimuthally averaged radial profiles of exposure maps in
the 10 energy bands.

We simultaneously fit the observed X-ray surface brightness profiles in the 10
energy bands with our model using affine-invariant MCMC sampling (Goodman &
Weare, 2010) implemented by the emcee python package (Foreman-Mackey et al.,
2013). The log-likelihood function for the data is defined by (up to a
normalization constant)

$-2\ln\mathcal{L}=\sum_{i,j}\frac{\left[d_{ij}-\left(w_{ij}T_{i}\widehat{S}_{\mathrm{X},ij}+{\cal
N}_{i}\times\mathrm{BGD}_{ij}\right)\right]^{2}}{\sigma_{ij}^{2}},$ (46)

where $i$ and $j$ run over all energy bands and all radial bins, respectively,
$d_{ij}$ is the binned X-ray brightness measured in units of counts per pixel,
$\sigma_{ij}$ is the statistical uncertainty of the measurement in each bin,
$T_{i}$ represents the Galactic transmission in the $i$th energy band
calculated by xspec using the photoionization cross sections of Verner et al.
(1996), $\widehat{S}_{\mathrm{X},ij}$ is the model prediction in each bin for
the X-ray surface brightness given by Equation (45), $w_{ij}$ is the
conversion factor proportional to the product of the effective area and the
net exposure time in each bin, $\mathrm{BGD}_{ij}$ denotes the background
contribution in each bin given in units of counts per pixel, and ${\cal
N}_{i}$ is a dimensionless calibration factor of the background in the $i$th
energy band.

Table 9: X-ray model parameter constraints derived from forward modeling of Chandra observations Parameter | Posterior summary
---|---
$n_{\mathrm{e0}}$ ($10^{-3}$ cm-3) | $5.03\pm 0.09$
$r_{\mathrm{c}}$ ($h^{-1}$ kpc) | $190\pm 7$
$\beta$ | $0.693\pm 0.019$
$T_{0}$ (keV) | $9.63\pm 2.38$
$r_{\mathrm{t}}$ ($h^{-1}$ kpc) | $176\pm 100$
$a$ | $-0.042\pm 0.248$
$c$ | $0.20\pm 0.44$
${\cal N}(0.50\mathrm{-}0.75~{}\mathrm{keV})$ | $0.94\pm 0.06$
${\cal N}(0.75\mathrm{-}1.00~{}\mathrm{keV})$ | $0.99\pm 0.11$
${\cal N}(1.00\mathrm{-}1.25~{}\mathrm{keV})$ | $0.74\pm 0.10$
${\cal N}(1.25\mathrm{-}1.50~{}\mathrm{keV})$ | $0.76\pm 0.10$
${\cal N}(1.50\mathrm{-}2.00~{}\mathrm{keV})$ | $0.90\pm 0.06$
${\cal N}(2.00\mathrm{-}3.00~{}\mathrm{keV})$ | $1.01\pm 0.04$
${\cal N}(3.00\mathrm{-}4.00~{}\mathrm{keV})$ | $0.98\pm 0.05$
${\cal N}(4.00\mathrm{-}5.00~{}\mathrm{keV})$ | $1.02\pm 0.04$
${\cal N}(5.00\mathrm{-}6.00~{}\mathrm{keV})$ | $1.01\pm 0.04$
${\cal N}(6.00\mathrm{-}7.00~{}\mathrm{keV})$ | $1.00\pm 0.03$

Our model has a total of 17 parameters, of which seven parameters describe the
cluster X-ray emission (see Equation (43)), namely
$(n_{\mathrm{e0}},r_{\mathrm{c}},\beta,T_{0},r_{\mathrm{t}},a,c)$, and the
rest are calibration nuisance parameters, $\\{{\cal
N}_{i}\\}_{i=1}^{N_{\mathrm{spec}}}$. For the parameters describing the
cluster X-ray emission, we use uninformative uniform priors of
$n_{\mathrm{e0}}\in[0,1]$ cm-3,
$r_{\mathrm{c}}\in[0,350]\,h^{-1}\,\mathrm{kpc}$, $\beta\in[0,2]$,
$T_{0}\in[1,20]$ keV, $r_{\mathrm{t}}\in[0,350]\,h^{-1}\,\mathrm{kpc}$,
$a\in[-0.5,0.5]$, and $c\in[-1,1]$. For each calibration parameter $\\{{\cal
N}_{i}\\}_{i=1}^{N_{\mathrm{spec}}}$, we adopt a Gaussian prior of $1.0\pm
0.2$. We sample the posterior PDFs of all model parameters over the full
parameter space allowed by the priors. Posterior summaries of the model
parameters are listed in Table 9. In Appendix E, we show the Chandra X-ray
brightness profiles along with the best-fit model.

Finally, we use the posterior samples obtained with the MCMC algorithm to
derive constraints on the gas mass $M_{\mathrm{gas}}(<r)$ enclosed within the
spherical radius $r$ and the hydrostatic equilibrium mass
$M_{\mathrm{HE}}(<r)$ of A370. The hydrostatic mass $M_{\mathrm{HE}}(<r)$ is
given by

$M_{\mathrm{HE}}(<r)=-\frac{k_{\mathrm{B}}T_{\mathrm{gas}}(r)r}{G\mu_{\mathrm{g}}m_{\mathrm{p}}}\left[\frac{d\ln\rho_{\mathrm{gas}}(r)}{d\ln
r}+\frac{d\ln T_{\mathrm{gas}}(r)}{d\ln r}\right],$ (47)

where $k_{\mathrm{B}}$ is the Boltzmann constant and $\mu_{\mathrm{g}}\approx
0.60$ is the mean molecular weight. We will compare the resulting hydrostatic
mass profile, $M_{\mathrm{HE}}(<r)$, with our weak-lensing results in Section
7.1.

In the top panel of Figure 14, we show the marginalized $1\sigma$ confidence
region of $T_{\mathrm{3D}}(r)$ obtained from our forward modeling, along with
the deprojected temperatures inferred from our spectral deprojection analysis
(Section 6.2). Similarly, we compare in the bottom panel of Figure 14 our
determinations of $M_{\mathrm{gas}}(<r)$ from both methods. For both
comparisons, we find that the two complementary approaches yield highly
consistent results.

##  7 Discussion

###  7.1 Hydrostatic Mass Bias

Figure 15: Comparison of spherically enclosed total mass profiles
$M_{\mathrm{tot}}(<r)$ of A370 as a function of spherical radius $r$. The
orange and blue shaded areas represent the marginalized $1\sigma$ confidence
regions from triaxial NFW modeling of the $\Sigma(\mbox{\boldmath$\theta$})$
map (Figure 9) with and without using the LOS prior, respectively. The dashed
line shows the posterior mean from spherical NFW modeling of the
$\Sigma(\mbox{\boldmath$\theta$})$ map. The red hatched area represents the
marginalized $1\sigma$ confidence region of the total hydrostatic mass
$M_{\mathrm{HE}}(<r)$ obtained from forward modeling of Chandra X-ray data
that cover the radial range $\theta\in[10\arcsec,200\arcsec]$.  Figure 16:
X-ray hydrostatic mass bias $b(r)=1-M_{\mathrm{HE}}(<r)/M_{\mathrm{WL}}(<r)$
in A370 as a function of spherical radius $r$ derived from combined Chandra
X-ray and Subaru weak-lensing observations. The orange and blue shaded areas
represent the marginalized $1\sigma$ confidence regions based on triaxial NFW
modeling of the $\Sigma(\mbox{\boldmath$\theta$})$ map (Figure 9) with and
without using the LOS prior, respectively. The dashed line shows the posterior
mean based on spherical NFW modeling of the $\Sigma(\mbox{\boldmath$\theta$})$
map.

Hydrostatic mass estimates $M_{\mathrm{HE}}$ are expected to be biased low,
depending on details of nongravitational processes and the level of residual
gas motions in the ICM. Determining the level of hydrostatic mass bias for a
representative sample of galaxy clusters has important implications for both
cluster cosmology and astrophysics (Planck Collaboration et al., 2014, 2016;
Pratt et al., 2019). Cosmological hydrodynamical simulations suggest a modest
level of hydrostatic mass bias for an ensemble of galaxy clusters,
$b_{\mathrm{HE}}\equiv 1-\langle M_{\mathrm{HE}}/M_{\mathrm{true}}\rangle\sim
5\%$–$20\%$ at $r\leqslant r_{500}$ (Nagai et al., 2007; Lau et al., 2009;
Meneghetti et al., 2010b; Nelson et al., 2012; Angelinelli et al., 2020;
Ansarifard et al., 2020), defined with respect to the true enclosed mass
$M_{\mathrm{true}}$. Since A370 is a highly disturbed system, the cluster is
likely to exhibit a higher than typical value of mass bias, which could serve
as an extreme limit expected for galaxy clusters.

With the aim of characterizing the level of hydrostatic mass bias in A370, we
compare our lensing-based determinations of the cluster mass profile (Section
5) to the hydrostatic mass profile $M_{\mathrm{HE}}(<r)$ derived from Chandra
X-ray data (Section 6.3). For this purpose, we compute the total mass
$M_{\mathrm{tot}}(<r)$ of a triaxial halo enclosed within a sphere of radius
$r$:

$M_{\mathrm{tot}}(<r)=\iiint\limits_{\cal
V}\\!\rho(X,Y,Z)dV=\int\limits_{4\pi}\\!d\Omega\int_{0}^{r}\\!r^{\prime 2}\rho
dr^{\prime},$ (48)

where $\rho(X,Y,Z)$ is the density function (Equations (29) and (30)), the
region of integration ${\cal V}$ is $\sqrt{X^{2}+Y^{2}+Y^{2}}\leqslant r$, and
$d\Omega=\sin\vartheta d\vartheta d\phi$ is the solid angle element in
spherical coordinates.

In Figure 15, we compare the spherically enclosed total mass profiles
$M_{\mathrm{tot}}(<r)$ obtained from our WL-2D and X-ray analyses. Here we
have extrapolated the X-ray forward model beyond the range of the fitted data
($\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}720\,h^{-1}\,\mathrm{kpc}$) to
compute $M_{\mathrm{HE}}(<r)$ out to larger cluster radii. We note that in
contrast to the triaxial lensing constraints on $M_{\mathrm{tot}}(<r)$, the
hydrostatic mass $M_{\mathrm{HE}}(<r)$ obtained assuming spherical symmetry is
not corrected for the projection effect due to the LOS elongation of the gas
distribution. Since the shape of the collisional gas is rounder than the
underlying matter (e.g., Suto et al., 2017), the level of projection bias in
the gas distribution is expected to be less than $\sim 10\%$ found in the
total mass distribution (see Section 5.4).

In Figure 16, we show the hydrostatic mass bias as a function of spherical
radius $r$, defined with respect to the total mass $M_{\mathrm{WL}}(<r)$
determined from weak lensing:

$b(r)=1-\frac{M_{\mathrm{HE}}(<r)}{M_{\mathrm{WL}}(<r)}.$ (49)

The results are shown for the three different priors on the halo shape
employed in our mass modeling of the WL-2D data (see Table 6). We find no
significant evidence for a strong variation of $b(r)$ both within and beyond
the radial range probed by the Chandra data ($\lower
2.15277pt\hbox{$\;\buildrel<\over{\sim}\;$}720\,h^{-1}\,\mathrm{kpc}$). At
each radius $r$, we find similar central values of the distributions from the
spherical and the fiducial triaxial cases (see Table 6).

From the triaxial lens modeling, we obtain mass ratios of $1-b(r)=0.56\pm
0.09$ and $0.51\pm 0.09$ at $r=0.7\,h^{-1}\,\mathrm{Mpc}\sim 0.7r_{500}$, with
and without using the LOS prior, respectively. When the X-ray forward model is
extrapolated out to $r_{500}\sim 1\,h^{-1}\,\mathrm{Mpc}$, we find
$1-b(r_{500})=0.54\pm 0.12$ and $0.50\pm 0.11$ with and without the LOS prior,
respectively. The range of mass bias inferred for A370, $b\in[0.34,0.61]$ at
the $1\sigma$ level, is on the high side of the distribution expected from
cosmological cluster simulations (see Nagai et al., 2007; Lau et al., 2009;
Ansarifard et al., 2020) and is in better agreement with the value of
$1-b_{\mathrm{HE}}=\langle M_{\mathrm{HE}}/M_{\mathrm{true}}\rangle=0.58\pm
0.04$ required to bring the Planck CMB and cluster constraints into full
agreement in the base $\Lambda$CDM cosmology of Planck Collaboration et al.
(2016). However, it should be noted again that the mass bias found for this
cluster should be considered as an extreme value expected for galaxy clusters.

###  7.2 Gas Mass Fraction