© The Authors 2025

1. Introduction

The formation of large-scale structures in the Universe, predicted by N-body simulations based on the Λ cold dark matter (ΛCDM) cosmological model, follows a “bottom-up” process. In this scenario, smaller structures form first and gradually merge into larger ones under the influence of gravity. These hierarchical structures make up the cosmic web (Zel’dovich 1970), where significant concentrations of matter accumulate at the intersections, leading to the formation of massive dark matter (DM) halos. They host the formation of galaxy clusters, with lower-mass halos mainly merging along filamentary structures (Peebles 1980). Galaxy clusters, the largest self-gravitationally bound systems observed, play a critical role in shaping the formation and evolution of the Universe. Therefore, studying the properties of these clusters would provide crucial insights into the growth of cosmic structures and the nature of DM.

The key properties that characterize galaxy clusters, such as the halo mass function (e.g., Press & Schechter 1974; Sheth & Tormen 1999; Jenkins et al. 2001; Tinker et al. 2008; Watson et al. 2013) and the density profile (e.g., Navarro et al. 1997; Bullock et al. 2001; Macciò et al. 2007; Duffy et al. 2008; Mandelbaum et al. 2008; Ludlow et al. 2016; Shan et al. 2017; Child et al. 2018), can be used to test and refine our understanding of structure formation and the evolution of large-scale structures, and are thus essential for constraining cosmological models (Ghirardini et al. 2024). A variety of methods have been developed to measure the properties of halos, each with its strengths and limitations. X-ray observations, for example, allow for the measurement of the total mass of clusters based on the temperature and flux profiles of the hot gas assuming spherical symmetry and hydrostatic equilibrium (Nagai et al. 2007), which can introduce biases in clusters that deviate from spherical shapes and equilibrium. The Sunyaev-Zel’dovich (SZ) effect links the ionized gas in clusters with cosmic microwave background (CMB) observations, and provides crucial, redshift-independent observables that directly probe the integrated electron pressure of the intracluster medium, enabling us to derive fundamental cluster properties such as the mass and thermodynamic state of galaxy clusters across cosmic time (Planck Collaboration XXIV 2016). However, point sources and foreground signals could contaminate the SZ effect and complicate the measurements (Carlstrom et al. 2002; Vale & White 2006). Other methods, such as dynamics-based estimates, focus on galaxy velocity dispersion and distribution (Brainerd 2005; Sohn et al. 2020), but are also sensitive to the galaxy population within clusters (Nagai & Kravtsov 2005).

Weak gravitational lensing (WL) provides an independent and powerful method of studying galaxy clusters. Unlike other techniques, WL does not rely on assumptions about the hydrodynamic state or symmetry of the clusters (e.g., Shan et al. 2018; Li et al. 2023a). The lensing effect arises when the light from background galaxies is distorted due to the gravitational field of massive foreground structures, including galaxy clusters (Bartelmann & Schneider 2001; Hoekstra & Jain 2008; Fu & Fan 2014; Kilbinger 2015; Schrabback et al. 2025; Diego et al. 2025). These distortions thus offer a direct way to probe the distribution of DM in the clusters and provide valuable information about their mass, shape, and orientation. Recent WL surveys, including the Canada–France Hawaii Telescope Lensing Survey (CFHTLenS¹; Heymans et al. 2012), the Dark Energy Survey (DES²; Dark Energy Survey Collaboration 2016), the Kilo-Degree survey (KiDS³; de Jong et al. 2013), the Hyper Suprime-Cam Subaru Strategic Program survey (hereafter the HSC-SSP⁴ survey; Aihara et al. 2018), and the Euclid space telescope of the European Space Agency (ESA⁵; Laureijs et al. 2011), have significantly improved our ability to measure cosmological parameters and the properties of DM halos. These surveys make WL a promising tool for addressing long-standing questions about DM and the growth of large-scale structures in the Universe.

One specific technique, galaxy–galaxy lensing (GGL), measures the correlation between the positions of foreground galaxies and the shear of background galaxies (Mandelbaum et al. 2006a). This technique provides valuable clues about the relationship between galaxies and their surrounding DM halos (Mandelbaum et al. 2005, 2006b, 2013; Cacciato et al. 2009; Zu & Mandelbaum 2015; Shan et al. 2017; Cui et al. 2018; Xu et al. 2021; Luo et al. 2022, 2024; Wang et al. 2024; Li et al. 2024). However, systematic errors such as miscentering and nonsphericity still affect the accuracy of mass estimates in weak-lensing studies. Miscentering occurs when the observed center of the cluster differs from the center of the lensing potential (hereafter the lensing center), leading to biases in the measurement of mass and shape (George et al. 2012). Nonsphericity, which reflects the fact that halos are often triaxial rather than spherical (Jing & Suto 2002), can also introduce projection effects, which affect mass estimates (Osato et al. 2018; Gonzalez et al. 2021) and offer important insights into the nature of DM particles (Peter et al. 2013).

To address these challenges, it is essential to employ two-dimensional (2D) WL analyses (Cypriano et al. 2004; Deb et al. 2010; Oguri et al. 2010), by considering the full 2D shear map. Incorporating realistic cluster halo features like ellipticity and miscentering, 2D WL can provide more robust estimates of individual halo properties compared to traditional stacking methods. This study presents one of the first large-sample applications of 2D WL analysis to investigate the properties of 299 massive redMaPPer clusters with the shear data from the first year of the HSC-SSP survey, which significantly enhances statistical capabilities. We individually fit the observed shear map of each cluster to the elliptical Navarro-Frenk-White (NFW) (Navarro et al. 1996, 1997) models. This facilitates measuring key cluster parameters such as concentration, ellipticity, and mass centroid position, revealing new perspectives on the structure and evolution of massive galaxy clusters.

The structure of this paper is as follows. We describe the observational data, including the lens and source samples in Section 2, and introduce the methodology and fitting procedure in Section 3. In Section 4, we show our results. Finally, conclusions are presented in Section 5. Throughout the paper, we assume a flat Universe with the matter density Ω_m = 1 − Ω_Λ = 0.315, the dimensionless Hubble constant h = 0.674, the baryon density Ω_bh² = 0.0224, the spectral index n_s = 0.965, and the normalization of the matter power spectrum σ₈ = 0.811 (Planck Collaboration VI 2020). Unless otherwise stated, we quote parameter constraints using the median of the posterior distribution as the best-fit value, with uncertainties corresponding to the 68% credible interval.

2. Data

2.1. Cluster sample

In this paper, we studied the properties of clusters from the red-sequence Matched-filter Probabilistic Percolation (redMaPPer) cluster finding algorithm (version 6.3, Rykoff et al. 2014), found in the Sloan Digital Sky Survey (SDSS) Data Release 8 (DR8, Aihara et al. 2011) catalog. This catalog covers a 10 000 deg² sky area and contains 26 111 galaxy clusters over the redshift range z ∈ [0.08, 0.55], each assigned a richness parameter λ. The parameter λ is the sum of the membership probabilities of galaxies within a richness-scaling radius R_λ = (1.0 h⁻¹ Mpc)(λ/100.0)^0.2, which is proportional to halo mass.

We selected clusters with spectroscopic redshift measurements to ensure that we could precisely determine their positions along the lines of sight and require that the richness λ > 20, yielding a sample of 15 657 massive clusters. For each redMaPPer cluster, five galaxy candidates were assigned probabilities, P_cen, indicating their likelihood of being the optical center galaxy. We then designated the galaxy with the highest probability as the observed center. The displacement between this observed center and the center determined by weak lensing was referred to as the brightest cluster galaxy (BCG) offset. For a full description of the redMaPPer v6.3 SDSS DR8 catalog, refer to Rykoff et al. (2016).

2.2. Weak-lensing data

The HSC-SSP is an extensive, wide-field multiband imaging survey designed to tackle a broad range of scientific inquiries, spanning from cosmology to solar system bodies (Aihara et al. 2018). A key objective of this survey is the weak lensing study of wide-layer observations, targeting a sky coverage of approximately 1400 deg². The wide-layer data in the latest third data release has achieved a full depth of ∼26 mag at 5σ across all five bands (grizy) for an area of about 670 deg² (Aihara et al. 2022). The three-year shear catalog has also been completed (Li et al. 2022), but was only made publicly available very recently.

Therefore, in this paper, we utilized the first-year HSC-SSP shear catalog (S16A) (Mandelbaum et al. 2018). The shapes of galaxies were obtained using the re-Gaussianization point spread function (PSF) correction method on the i-band coadded images (Hirata & Seljak 2003). To ensure the robustness of shear measurements and exclude spurious detections as well as blends that may introduce contamination, we retained only those galaxies that meet the selection criteria specified in Table 4 of Mandelbaum et al. (2018). The photometric redshift (photo-z) for each galaxy was determined on the basis of HSC five-band photometry, and the specific catalog derived from the direct empirical photometric code (DEMP; Hsieh & Yee 2014) was used in our study (Tanaka et al. 2018). It has been demonstrated to exhibit high accuracy, with a scatter of σ[Δz_p/(1 + z_p)] ∼ 0.05 and an outlier rate of ∼15% for galaxies in the redshift range of 0.2 ≤ z_p ≤ 1.5, where z_p is the optimal photo-z value derived from its probability distribution for a given galaxy. Here, we adopted z_p as the redshift estimate z_s for each source galaxy, and then selected the lens (cluster)-source pairs that satisfy the criterion Δz = z_s − z_l > 0.2, where z_l = z_cl denotes the spectroscopic redshift of the cluster studied. This criterion ensures a clear separation between the foreground cluster and background galaxies. Furthermore, to ensure an adequate signal-to-noise ratio, we required the average number density of source galaxies to be ≥10/arcmin², within a circular region centered on the observed center of the foreground cluster and extending to a radius of 0.5 deg.

After cross-matching the redMaPPer clusters with the sky coverage of the HSC-SSP survey and applying all the aforementioned selection criteria, we ultimately obtained a sample consisting of 299 lens (cluster)-source pairs. The normalized redshift distributions of the entire lens and source samples are presented in Fig. 1.

Fig. 1. Normalized distribution of spectroscopic redshifts for the lens sample (dashed purple) and photometric redshifts for the source sample (solid blue).

3. Methods

3.1. Lensing signal

The weak-lensing effect can be characterized by the second derivatives of the lensing potential, ψ, via the Jacobian matrix, 𝒜 (Bartelmann & Schneider 2001):

$$\begin{array}{c}\hfill \mathcal{A}=({\delta }_{\mathit{ij}}-\frac{{\partial }^2\psi (\mathit{\theta })}{\partial {\theta }_i\partial {\theta }_j})=\left(\begin{array}{cc}1-\kappa -{\gamma }_1& -{\gamma }_2\\ -{\gamma }_2& 1-\kappa +{\gamma }_1\end{array}\right),\end{array}$$(1)

where the convergence, κ, and lensing shear, γ, in the complex form of γ₁ + iγ₂, respectively, induce an isotropic size change and elliptical shape distortion on the observed image of the background galaxy. Both of them are related to the second derivatives of the lensing potential, and are thus not independent of each other.

The convergence, specifically, is directly related to the projection of line-of-sight density fluctuations weighted by the lensing efficiency factor, which is expressed as the ratio Σ/Σ_crit, with Σ indicating the projected mass density of the lens and Σ_crit denoting the critical mass density.

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{crit}}=\frac{c^2}{4\pi G}\frac{D_{\mathrm{s}}}{D_{\mathrm{l}}{D}_{\mathrm{ls}}},\end{array}$$(2)

where D_s, D_l, and D_ls are the angular diameter distances from the observer to the source, from the observer to the lens, and from the lens to the source, respectively. c denotes the speed of light, and G stands for the Newtonian gravitational constant. As an indicator of lensing efficiency, the critical mass density, Σ_crit, for each galaxy cluster was calculated from the redshift distribution of its background sources. On the other hand, the observables for a weak-lensing survey are often the ellipticities of galaxies, which encompass both the intrinsic shape noise of the galaxies and the weak lensing shears γ, or more precisely, the reduced shears g = γ/(1 − κ) (Seitz & Schneider 1997).

The first-year HSC-SSP shear catalog was produced on the coadded i-band images using the re-Gaussianization PSF correction method (Hirata & Seljak 2003; Mandelbaum et al. 2018). The core idea of estimating galaxy shape with this method is to apply a Gaussian profile featuring elliptical isophotes to fit the image, and the components of the distortion are defined as follows:

$$\begin{array}{c}\hfill (e_1,e_2)=\frac{1-{(b/a)}^2}{1+{(b/a)}^2}(cos2\varphi ,sin2\varphi ),\end{array}$$(3)

where b/a is the axis ratio and ϕ is the position angle of the major axis with respect to the equatorial coordinate system (sky coordinates). The ensemble average distortion is then an estimator for the reduced shear, g,

$$\begin{array}{c}\hfill ({\widehat{g}}_1,{\widehat{g}}_2)=\frac{1}{2\mathcal{R}}⟨(e_1,e_2)⟩,\end{array}$$(4)

where ℛ is known as the “shear responsivity” and represents the response of the distortion to a small shear (Kaiser et al. 1995; Bernstein & Jarvis 2002; Jarvis et al. 2003). $\mathcal{R}\approx 1-e_{\mathrm{rms}}^2$, where e_rms is the rms intrinsic distortion per component, typically e_rms ≃ 0.4.

To carry out comprehensive analyses of our cluster sample, we worked directly on the 2D distortion maps. For all the clusters, we adopted a pixel size of 1 × 1 arcmin² centered on the observed center (the position of the BCG), and the lensing distortion was measured in a statistical sense, i.e., by averaging the galaxy shapes over an arcminute scale to mitigate the impact of dominant intrinsic shape noise (Oguri et al. 2010). Due to their varying masses and redshifts, different clusters occupy different angular sizes on the sky. Therefore, for each cluster, we gave a first estimate of its mass M_est using the richness-mass relation provided in Rykoff et al. (2012), calculated the cosmic matter density ρ_m(z_cl) at the cluster’s redshift z_cl, and derived an estimate of halo radius $r_{200\mathrm{m}}=\sqrt[3]{3M_{\mathrm{est}}/(4\pi ·200{\rho }_{\mathrm{m}}(z_{\mathrm{cl}}))}$, which is defined as the radius within which the mean halo density is 200 times the cosmic matter density ρ_m(z_cl). The total number of pixels on each side could be subsequently obtained by rounding up the angular size θ_cl = 2 × r_200m/D_A(z_cl) of each cluster in units of arcminute, where D_A(z_cl) is the angular diameter distance of the cluster. Pixels without any background galaxies were treated as masked and were excluded from the fitting process; these regions naturally appear blank in the shear map (e.g., Fig. 2).

Fig. 2. Upper panels: The left panel shows the observed 2D shear map for an illustrative redMaPPer cluster (M ≃ 1014.73 h−1 M⊙), while the right panel displays the corresponding best-fitting shear field predicted by the elliptical NFW model. Each map covers an area of 18′×18′. In both maps, the orientation and length of the sticks within each 1′×′ pixel indicate the local tangential direction and amplitude of the distortion, respectively. For visualization purposes, both shear fields have been smoothed with a Gaussian kernel with a full width at half maximum (FWHM) of ≈2′ (Oguri et al. 2010). The blue dot marks the weak-lensing-derived center, and the red pentagram denotes the observed center. Lower panel: Posterior constraints on structural parameters of the cluster, including the halo concentration, ellipticity, and projected offsets of the mass centroid. The contours represent the 1σ and 2σ confidence levels.

We followed the approach outlined in Mandelbaum et al. (2018) to compute the weighted shear responsivity factor, $\widehat{\mathcal{R}}$, utilizing the per-object estimates of rms distortion e_rms, i within each 1 × 1 arcmin² pixel and the inverse variance weights w_i which encode the shape-measurement uncertainty and intrinsic ellipticity dispersion for each galaxy. Specifically, the shear responsivity factor $\widehat{\mathcal{R}}({\mathit{\theta }}_l)$ in the l-th pixel (its angular position θ_l) is given by

$$\begin{array}{c}\hfill \widehat{\mathcal{R}}({\mathit{\theta }}_l)=1-{\displaystyle \frac{{\sum }_{i\in {\mathit{\theta }}_l}{w}_i{e}_{\mathrm{rms},i}^2}{{\sum }_{i\in {\mathit{\theta }}_l}{w}_i}},\end{array}$$(5)

where i ∈ θ_l indicates that the summation runs over source galaxies contained within the l-th pixel.

Similarly, the ensemble estimates for calibration biases in the l-th pixel, consisting of both the multiplicative term $\widehat{m}({\mathit{\theta }}_l)$ and additive term ${\widehat{c}}_{\alpha }({\mathit{\theta }}_l)$, are also derived as a weighted sum of their respective catalog estimates m_i and c_α, i within the l-th pixel,

$$\begin{array}{cc}& \widehat{m}({\mathit{\theta }}_l)={\displaystyle \frac{{\sum }_{i\in {\mathit{\theta }}_l}{w}_i{m}_i}{{\sum }_{i\in {\mathit{\theta }}_l}{w}_i}},\hfill \end{array}$$(6)

$$\begin{array}{cc}& {\widehat{c}}_{\alpha }({\mathit{\theta }}_l)={\displaystyle \frac{{\sum }_{i\in {\mathit{\theta }}_l}{w}_i{c}_{\alpha ,i}}{{\sum }_{i\in {\mathit{\theta }}_l}{w}_i}},\hfill \end{array}$$(7)

where α = 1, 2, denoting the two components of distortion throughout the paper.

Finally, the reduced shear in the l-th pixel is estimated as (Oguri et al. 2010; Mandelbaum et al. 2018)

$$\begin{array}{c}\hfill g_{\alpha }({\mathit{\theta }}_l)={\displaystyle \frac{{\sum }_{i\in {\mathit{\theta }}_l}{w}_i{e}_{\alpha ,i}}{2\mathcal{R}({\mathit{\theta }}_l)(1+\widehat{m}({\mathit{\theta }}_l)){\sum }_{i\in {\mathit{\theta }}_l}{w}_i}}-\frac{{\widehat{c}}_{\alpha }({\mathit{\theta }}_l)}{1+\widehat{m}({\mathit{\theta }}_l)}·\end{array}$$(8)

For the angular scales of interest, the weak-lensing approximation holds, with g_α ≈ γ_α. In subsequent analyses, this approximation was used only when constructing the large-scale-structure covariance (see Eq. (22)); both the model and the data vectors were expressed in terms of the reduced shear g_α.

To estimate shape noise per pixel, we randomly rotated source galaxies to eliminate correlated lensing signals and applied the same procedures described above to generate 2D noise distortion maps. For each cluster, we performed 1000 random rotations and consequently obtained 1000 noise distortion maps. The shape noise σ_γ(θ_l; α), also comprising two components, was then determined by calculating the standard deviation for each pixel across these 1000 noise distortion maps.

3.2. Cluster mass model

Clusters in the collisionless cold dark matter (CDM) Universe typically exhibit a highly nonspherical structure that was well fitted by a triaxial density profile (e.g., Jing & Suto 2002). We thus constructed an elliptical lens model (Oguri et al. 2010) by incorporating an ellipticity into the isodensity contour, under the assumption that the cluster mass distribution can be described by a single halo component whose radial distribution follows the NFW density profile. Such an elliptical model provided a better description of halos compared to the spherical one (Jing & Suto 2002; Kasun & Evrard 2005; Allgood et al. 2006), since the projection of a triaxial halo along any direction results in elliptical isodensity contours on the convergence map (Oguri et al. 2003; Oguri & Keeton 2004; Heyrovský & Karamazov 2024). Specifically, the following mass model was adopted in the analyses:

$$\begin{array}{cc}& \kappa (x,y)={\kappa }_{\mathrm{sph}}(\eta ),\hfill \end{array}$$(9)

$$\begin{array}{cc}& {\eta }^2=x^{\prime 2}+\frac{y^{\prime 2}}{{(1-e)}^2},\hfill \end{array}$$(10)

$$\begin{array}{cc}& x\prime =xcos{\theta }_e+ysin{\theta }_e,\hfill \end{array}$$(11)

$$\begin{array}{cc}& y\prime =-xsin{\theta }_e+ycos{\theta }_e,\hfill \end{array}$$(12)

where κ_sph(η) denotes the radial convergence for the spherical NFW model. The halo ellipticity e is defined as e = 1 − b/a, where a and b are the major and minor axis lengths of the isodensity contour. Here we adopted a coordinate system where the coordinate origin is set at the observed center of the cluster (x_c, y_c) and the position angle θ_e measured counterclockwise from the positive x axis.

To calculate κ_sph(η), we adopted the NFW density profile (Navarro et al. 1996, 1997), which is given by

$$\begin{array}{c}\hfill \rho (r)=\frac{{\rho }_s}{(r/r_s){(1+r/r_s)}^2},\end{array}$$(13)

where ρ_s and r_s are the characteristic mass density and scale radius, respectively. It is completely characterized by two parameters, the mass and the halo concentration. The corresponding convergence κ_sph(η) = (2ρ_sr_s/Σ_crit)f(x = η/r_s) with (Wright & Brainerd 2000)

$$\begin{array}{c}\hfill f(x)=\{\begin{array}{cc}\frac{1}{x^2-1}[1-\frac{2}{\sqrt{1-x^2}}\text{arctanh}\sqrt{\frac{1-x}{1+x}}]\hfill & x<1\hfill \\ \frac{1}{3}\hfill & x=1\hfill \\ \frac{1}{x^2-1}[1-\frac{2}{\sqrt{x^2-1}}\text{arctan}\sqrt{\frac{x-1}{1+x}}]\hfill & x>1\hfill \end{array}.\end{array}$$(14)

To compute the lensing efficiency, we used the spectroscopic redshift z_cl of each cluster and the point-estimate photometric redshifts z_p, i of background galaxies to calculate the critical surface mass density Σ_crit, i for each source galaxy, according to Eq. (2). The effective critical surface density was then obtained by averaging over all selected background sources within the cluster region using lensing weights w_i:

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathrm{crit},\mathrm{eff}}={\displaystyle \frac{{\sum }_{i\in {\mathit{\theta }}_{\mathrm{cl}}}{w}_i{\mathrm{\Sigma }}_{\mathrm{crit},i}}{{\sum }_{i\in {\mathit{\theta }}_{\mathrm{cl}}}{w}_i}},\end{array}$$(15)

where the summation is over all background galaxies within a θ_cl × θ_cl region of the cluster center.

The lensing shear for the elliptical mass model was calculated from κ(x, y) using the relationship between the convergence and the shear (Kaiser & Squires 1993; Seitz & Schneider 1995; Squires & Kaiser 1996). In particular, we used their relation in the Fourier space with

$$\begin{array}{c}\hfill \widehat{\mathit{\gamma }}(\mathit{k})={\pi }^{-1}\widehat{D}(\mathit{k})\widehat{\kappa }(\mathit{k}),\end{array}$$(16)

and

$$\begin{array}{c}\hfill \widehat{D}(\mathit{k})=\pi \frac{k_1^2-k_2^2+2i{k}_1{k}_2}{k_1^2+k_2^2}·\end{array}$$(17)

The reduced shear was then obtained as g = γ/(1 − κ).

In summary, the elliptical mass model was specified by six parameters:

$$\begin{array}{c}\hfill \mathit{p}\equiv \{M_{200\mathrm{m}},c_{200\mathrm{m}},e,\mathrm{\Delta }x,\mathrm{\Delta }y,{\theta }_e\},\end{array}$$(18)

where M_200m is the mass within the halo radius r_200m, c_200m is the concentration parameter defined as the ratio r_200m/r_s, e is the halo ellipticity, Δx and Δy represent the offsets along the x and y axes, respectively, from the observed cluster center (x_c, y_c), and θ_e is the position angle.

3.3. 2D weak-lensing fitting

To constrain cluster properties, we compared the pixelized 2D distortion field g_α(θ_l), described in Section 3.1, with the reduced shear g_α^m(θ_l; p) predicted by the elliptical mass model, where p denotes the parameter set defined in Section 3.2. Specifically, we calculated the χ² as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\chi }^2=& {\displaystyle \sum _{\alpha ,\beta =1}^2\sum _{l,n=1}^{N_{\mathrm{pix}}}}[g_{\alpha }({\mathit{\theta }}_l)-g_{\alpha }^{\mathrm{m}}({\mathit{\theta }}_l;\mathit{p})]{\left[{\mathit{C}}^{-1}\right]}_{\alpha \beta ,ln}\hfill \\ \hfill & \times [g_{\beta }({\mathit{\theta }}_n)-g_{\beta }^{\mathrm{m}}({\mathit{\theta }}_n;\mathit{p})],\hfill \end{array}\end{array}$$(19)

where the indices α and β run over the two components of reduced shear and the indices l and n denote the position of the pixels. C is the error covariance matrix and C⁻¹ is its inverse. It is expressed as

$$\begin{array}{c}\hfill \mathit{C}={\mathit{C}}^{\mathrm{shape}}+{\mathit{C}}^{\mathrm{lss}}.\end{array}$$(20)

The term C^shape is dominated by intrinsic shape noise and measurement error. This noise term is expected to be uncorrelated between different pixels; therefore, it is diagonal and is given by

$$\begin{array}{c}\hfill {\left[{\mathit{C}}^{\mathrm{shape}}\right]}_{\alpha \beta ,ln}={\delta }_{\alpha \beta }^K{\delta }_{\mathit{ln}}^K{\sigma }_{\gamma }^2({\mathit{\theta }}_l;\alpha ),\end{array}$$(21)

where δ_αβ^K and δ_ln^K represent the Kronecker delta function, and σ_γ(θ_l; α) denotes the shape noise estimated from the noise distortion maps, as described in Section 3.1.

The term C^lss presents the covariance due to the large-scale structure, which is given in Hoekstra (2003).

$$\begin{array}{c}\hfill {\left[{\mathit{C}}^{\mathrm{lss}}\right]}_{\alpha \beta ,ln}={\xi }_{\alpha \beta }(|{\theta }_l-{\theta }_n|),\end{array}$$(22)

the cosmic shear correlation function, ξ_αβ, is assumed to depend only on the length of the vector connecting the two points θ_l and θ_n, due to the statistical isotropy of the Universe. Specifically, the shear correlation functions ξ_± are defined through combinations of the two shear components as follows:

$$\begin{array}{cc}\hfill {\xi }_{11}(r)& ={(cos2\varphi )}^2{\xi }_{\mathit{tt}}(r)+{(sin2\varphi )}^2{\xi }_{\times \times }(r),\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill {\xi }_{22}(r)& ={(sin2\varphi )}^2{\xi }_{\mathit{tt}}(r)+{(cos2\varphi )}^2{\xi }_{\times \times }(r),\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill {\xi }_{12}(r)& =cos2\varphi sin2\varphi [{\xi }_{\mathit{tt}}(r)-{\xi }_{\times \times }(r)],\hfill \end{array}$$(25)

here, r = |θ_l − θ_n|, and ϕ represents the position angle between the coordinate x axis and the vector θ_l − θ_n. The functions ξ_tt and ξ_×× correspond to the tangential and cross-component shear correlation functions, respectively (Bartelmann & Schneider 2001). Both of the functions were calculated by the pyccl package⁶ (Chisari et al. 2019).

We utilized the Markov chain Monte Carlo (MCMC) sampler, emcee package⁷ (Foreman-Mackey et al. 2013), to explore the χ² surface, employing a standard Metropolis-Hastings sampling method with a multivariate Gaussian proposal distribution.

Among the six parameters in our model, we imposed a narrow Gaussian prior on the halo mass, with the mean value derived from the mass–richness relation presented in Rykoff et al. (2012),

$$\begin{array}{c}\hfill ln\left(\frac{M_{200\mathrm{m}}}{h_{70}^{-1}{10}^{14}{M}_{\odot }}\right)=A+B\times ln\left(\frac{\lambda }{60}\right),\end{array}$$(26)

where A = 1.72 and B = 1.08. The cluster richness, λ, and its measurement uncertainty, σ_λ, were taken from the redMaPPer catalog, and the corresponding mass uncertainty, σ_M, was obtained by propagating errors through the above relation. This allowed us to construct a cluster-specific Gaussian prior on halo mass. Our primary interest lay in constraining the structural properties of galaxy clusters, including halo concentration, ellipticity, and the offset between the WL mass centroid and the BCG position. However, uncertainties in halo mass could propagate into these structural parameters and bias their estimation–for example, a massive halo with a downward fluctuation in mass may be incorrectly assigned to a lower mass bin, and vice versa (Du & Fan 2014). Incorporating a mass prior thus helps reduce such biases and improves the robustness of our inference on the connection between halo structure and mass.

For the other five parameters, flat priors in the following ranges were adopted:

• 0 < c_200m < 20

• 0 < e < 0.9

• −θ_cl/2 < Δ_x < θ_cl/2

• −θ_cl/2 < Δ_y < θ_cl/2

• 0 < θ_e < π.

The best-fitting model parameters were determined by minimizing the χ² statistic, while constraints on individual parameters were derived by projecting the MCMC-sampled likelihood distributions onto the parameter space, with marginalization performed over the uncertainties of the other parameters. We calculated the dispersion (standard deviation) of parameters (concentration, ellipticity, and offset) from their posterior distributions of each cluster, and considered the inverse of the square of this dispersion as the weight to calculate the weighted average best-fit parameters in our subsequent analyses.

4. Results

In this section, we presented our main results. The results for individual clusters are described in Section 4.1, with the observed and fitted 2D shear map and the corresponding fitting results for one redMaPPer cluster shown as an example. The concentration-mass relation, the BCG offset distribution, and the cluster ellipticity distribution of our sample are discussed in Sections 4.2, 4.3, and 4.4, respectively.

4.1. Results for individual clusters

For our sample of 299 redMaPPer-selected clusters, we carefully constructed their 2D shear maps using the first-year HSC shear data, following the methodology described in Section 3.1. We then constrained the structural parameters of each cluster by directly fitting these shear maps with elliptical NFW models introduced in Section 3.2 via MCMC (Section 3.3). The analyses primarily focused on the halo concentration, ellipticity, and the distribution of miscentering offsets. The best-fit shear maps predicted by the elliptical NFW model closely reproduce the observed shear patterns, particularly in regions where the tangential shear signal has a high signal-to-noise ratio.

The observed 2D shear map (left) and its corresponding best-fit model (right) for a redMaPPer cluster are presented in the upper panels of Fig. 2 as a representative example. For visualization purposes only, both observed and fitted shear fields have been smoothed with a Gaussian kernel of FWHM ≈ 2′ (Oguri et al. 2010); please note that the unsmoothed shear fields were used in model fitting. In this example, the size of the map is 18′×18′. The orientation and length of the sticks in each 1′×1′ pixel represent the tangential direction and amplitude of the distortion, respectively. The blue dot and the red pentagram represent the center determined by weak lensing and the observed center (position of the BCG), respectively.

The lower panel of Fig. 2 presents posterior constraints on structural parameters, including concentration, ellipticity, and mass centroid offsets of the specific cluster. The contours are for 1σ and 2σ confidence levels. For each parameter, the best-fitting value and the corresponding 68% confidence interval are indicated above the respective one-dimensional marginalized distribution. The projected offsets Δ_x and Δ_y are expressed in arcseconds, while the position angle θ_e is measured in radians, defined counterclockwise from the positive x axis.

4.2. Concentration-mass relation

The concentration parameter characterizes the degree of central mass concentration within a DM halo or cluster. It depends on both the halo’s mass and redshift, and is closely linked to its formation and assembly history. A common observational approach to infer the average concentration is to stack lens–source pairs in bins of similar cluster mass and redshift, then fit the stacked shear profiles with parametric models (Xu et al. 2021).

Here, we measured the concentration of individual dark-matter halos and combined the results from all 299 clusters to statistically constrain the concentration–mass (c–M) relation. In Fig. 3, each gray point represents the best-fitting concentration and mass for a single cluster. We divided the range of log₁₀(M_200m) into five equal-width bins. Within each bin, the weighted average of the best-fit concentration and mass was computed, and the corresponding standard deviations are shown as black points with error bars.

Fig. 3. Concentration–mass (c–M) relation. Light gray points represent the best-fitting concentration and mass for individual clusters, while black points show the weighted average values in each mass bin, with error bars indicating the standard deviation. For comparison, the results for strong-lensing selected halos from Meneghetti et al. (2014) are shown by a solid orange line. Weak-lensing-based results from Xu et al. (2021) (upturn model) are plotted as a dashed-dotted blue line, and the power-law model from Duffy et al. (2008) derived from DM-only numerical simulations, is plotted as a dashed red line, respectively.

To place our results in context, we compared the derived c–M relation with several previous studies, as shown in Fig. 3. The dashed red line corresponds to the power-law model from Duffy et al. (2008), derived from DM-only numerical simulations. The dashed blue line corresponds to the upturn model proposed by Xu et al. (2021), based on weak-lensing measurements using the DECaLS DR8 shear catalog, with a particular emphasis on the high-mass end. We also included results based on strong-lensing selected halos. Meneghetti et al. (2014) estimated the concentrations of CLASH strong-lensing clusters through numerical simulations, using the projected mass distributions to fit NFW profiles, shown as the solid orange line. For a fair comparison, all reference models were evaluated at the mean redshift of our sample, z = 0.339, and we applied appropriate conversions to ensure that their definitions of halo mass and concentration are consistent with those adopted in this work.

As is shown in Fig. 3, our measured c–M relation lies closer to the strong-lensing-selected results of Meneghetti et al. (2014), appears steeper than the power-law relation of Duffy et al. (2008), and shows no indication of an upturn at the high-mass end. This may suggest that our sample selection and fitting methodology are more sensitive to halo concentration at the high-mass end. However, neither the upturn model derived from weak-lensing studies nor the power-law prediction obtained from DM-only numerical simulations can be ruled out within the current level of observational uncertainty.

We further investigated the dependence of concentration in mass-redshift space. In Fig. 4, gray points denote the mass and redshift of individual clusters in our sample. We divided the sample into bins with a redshift interval of 0.15 and a log₁₀(M_200m) interval of 0.2. Bins containing fewer than ten clusters were excluded from the analysis. Within each bin, the weighted mean concentration was computed and is shown as diamond symbols. The color of each diamond reflects the average concentration value in that bin, providing a visual representation of the concentration gradient across the mass–redshift plane.

Fig. 4. Concentration–mass–redshift (c–M–z) relation. Gray points indicate the positions of individual clusters in the mass–redshift plane. Colored diamonds represent the weighted mean concentration within each bin. Bins containing fewer than ten clusters are excluded from the analysis. Colored contour lines trace surfaces of constant concentration across the plane, and the color bar (labeled “concentration”) gives their numerical values.

By following the direction of gradient descent along the concentration contours, we find a clear trend: halos with lower redshifts and higher masses exhibit systematically lower concentrations. This trend is consistent with the expectation that more massive clusters at lower redshift are more dynamically relaxed and have had more time to evolve toward equilibrium.

To quantitatively characterize the joint dependence of concentration on both mass and redshift, we adopted a fitting function from Meneghetti et al. (2014) and applied it to the binned, weighted-mean concentrations in our sample:

$$\begin{array}{c}\hfill c(M,z)=A{\left(\frac{1+z_{\mathrm{mean}}}{1+z}\right)}^B{\left(\frac{M}{M_{\mathrm{mean}}}\right)}^C,\end{array}$$(27)

here, z_mean = 0.339 and M_mean = 3.058 × 10¹⁴ h⁻¹ M_⊙ represent the mean redshift and the mean mass of the sample, respectively. The best-fitting parameters and their associated uncertainties corresponding to the 68% confidence interval are summarized in Table 1.

Due to limitations in sample size and observational conditions, the precision of our quantitative characterization of the c–M–z relation is somewhat constrained. Nevertheless, the overall trend reveals a negative correlation of concentration with both halo mass and redshift.

4.3. Offset between lensing center and BCG

In observational studies, halo centers are typically identified with either the position of a massive galaxy (e.g., the BCG) or the peak of the hot gas distribution. This is because the lensing center of mass, dominated by invisible DM, is not directly observable. Miscentering between the true mass center and these observational tracers is a significant source of systematic uncertainty in stacked weak-lensing analyses (Johnston et al. 2007; Sommer et al. 2025). Correcting for this effect not only improves the accuracy of mass measurements but also enhances our understanding of the nature of DM (Harvey et al. 2015).

Weak lensing provides a unique and independent means of locating the lensing center of a DM halo, as it directly probes the projected gravitational potential and is sensitive to the full matter distribution. It thus offers a complementary perspective to optical indicators (e.g., BCGs) and thermal signals from the Sunyaev–Zel’dovich (SZ) effect (Ding et al. 2025).

In modeling the miscentering effect in weak lensing analyses, a 2D Gaussian distribution is widely adopted as a realistic representation of offset distributions. Following the approach of Oguri et al. (2018), we modeled the BCG offset distribution using a two-component 2D Gaussian function, as shown in Fig. 5, which enables a flexible description of both well-centered and miscentered populations:

Fig. 5. Upper panel: Probability distribution function (PDF) of the projected offset between the BCG position and the weak-lensing-determined halo center. Red data points show the observed BCG offsets in our sample, with error bars estimated via 1000 bootstrap resamplings. The solid blue and black lines represent the best-fitting single- and two-component 2D Gaussian models, respectively. The dashed purple and orange lines correspond to the well-centered and miscentered components of the two-component model. All offset scale parameters, including σ, σ1, and σ2, are expressed in units of h−1 kpc. Lower panels: Fitting results for the two-component (left) and single-component (right) 2D Gaussian models. Contours indicate 1σ and 2σ confidence levels, with best-fitting values and 68% confidence intervals shown above each 1D marginalized distribution.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill P(R_{\mathrm{off}})=& f_{\mathrm{cen}}\frac{R_{\mathrm{off}}}{{\sigma }_1^2}exp(-R_{\mathrm{off}}^2/2{\sigma }_1^2)\hfill \\ \hfill & +(1-f_{\mathrm{cen}})\frac{R_{\mathrm{off}}}{{\sigma }_2^2}exp(-R_{\mathrm{off}}^2/2{\sigma }_2^2).\hfill \end{array}\end{array}$$(28)

Here, R_off denotes the projected separation between the weak-lensing-determined halo center and the BCG position. The parameter f_cen represents the fraction of well-centered clusters, while σ₁ and σ₂ describe the characteristic offset scales of the well-centered and miscentered components, respectively. Parameter uncertainties were estimated via 1000 bootstrap resamplings. The reduced χ² values for the single- and two-component models are 8.58 and 0.87, respectively, confirming the superior fit of the two-component model.

The upper panel of Fig. 5 displays the weighted distribution of BCG–WL centroid offsets. The solid blue line represents a single 2D Gaussian fit with a characteristic scale of $\sigma ={203}_{-26}^{+16}{h}^{-1}$ kpc, which evidently fails to capture the bimodality. The best-fitting parameters for the two-component model are f_cen = 0.43 ± 0.05, ${\sigma }_1=115.2{\pm }_{-12}^{+9.9}{h}^{-1}$ kpc for the well-centered population, and ${\sigma }_2={409}_{-38}^{+27}{h}^{-1}$ kpc for the miscentered population.

We notice that, for a few clusters, the inferred offsets were close to the edge of the default fitting range [ − θ_cl/2, θ_cl/2]. To test the robustness of our offset measurements and to evaluate whether this boundary could bias the results, we repeated the fitting for these clusters using an expanded data range of [ − θ_cl, θ_cl], while keeping the offset prior unchanged. The resulting parameter constraints remain consistent within uncertainties, confirming the stability and reliability of the fitted offset distribution.

We further divided the sample into high- and low-mass subsamples based on the median halo mass. Offsets in the low-mass bin exhibit larger values of σ₁ and σ₂, suggesting that lower-mass halos in our sample may be more affected by surrounding large-scale structure and projection effects, leading to greater displacements between the WL center and the BCG.

Similarly, when dividing the sample into high- and low-redshift bins, we find that the low-redshift clusters tend to have higher values of f_cen and smaller offset scales. This trend is consistent with the expectation that halos become more dynamically relaxed over time, resulting in better alignment between the mass centroid and the BCG. These comparisons are illustrated in Fig. 6, which shows the projected offset distributions and the corresponding best-fit parameters for different redshift and mass bins.

Fig. 6. Probability distribution functions (PDFs) of the projected offsets between the BCG and the weak-lensing-determined halo center for different redshift and mass bins. All error bars are estimated from 1000 bootstrap resamplings. Red and blue points represent the low- and high-redshift subsamples, respectively, with solid lines showing the best-fitting two-component 2D Gaussian models. Corresponding vertical lines indicate the values of σ1 and σ2. Pink and light blue points show the low- and high-mass subsamples, respectively, with dashed lines and vertical markers indicating the best fits and offset scales.

We note that the miscentering parameters inferred from our BCG-WL offsets differ from those reported in previous studies. For example, the f_cen value shown in Fig. 5 is lower than that obtained by Oguri et al. (2018), while our offset scales (σ₁ and σ₂) are larger than their reported values of 46 ± 9 h⁻¹ kpc and 260 ± 40 h⁻¹ kpc, respectively. This discrepancy may arise from differences in center definitions and sample selection. Notably, Oguri et al. (2018) measured the offsets between BCGs and X-ray centers of clusters, with the X-ray temperatures estimated using a mass–temperature relation based on M₅₀₀, which corresponds to more compact regions than the M_200m used in our analysis. Motivated by the relatively large WL–BCG offsets, we further investigated whether this trend reflects genuine physical differences between weak-lensing and baryon-based tracers. As a first step, we compared the optical centers from redMaPPer with those from two other optical catalogs – HSC S19A (Aihara et al. 2022) and WH2022 (Wen & Han 2022). For each HSC or WH2022 cluster, we identified the nearest redMaPPer counterpart using the match_to_catalog_sky method in astropy’s SkyCoord class (Astropy Collaboration 2022), with a maximum separation of 1 arcmin. We find that more than 85% of the matched clusters have projected offsets smaller than 100 kpc, indicating that the optical-optical center differences are negligible.

We then compared the redMaPPer optical centers with those determined from SZ observations (ACT DR5; Naess et al. 2020) and X-ray data (Ghirardini et al. 2024), computing the projected physical offsets using redMaPPer photometric redshifts. In this calculation, we adopted a maximum offset threshold corresponding to a typical cluster size (∼1 h⁻¹ Mpc) to exclude spurious large separations. The resulting offset distributions and their best-fitting two-component models are shown in Fig. 7, along with several results from the literature (Ding et al. 2025; Oguri et al. 2018) for comparison. Although the exact values of f_cen, σ₁, and σ₂ vary among studies, the overall shapes of the offset distributions are broadly consistent.

Fig. 7. Probability distribution functions (PDFs) of the projected offsets between the BCG position and the SZ- or X-ray-determined cluster centers. Black and purple data points represent the offsets between the redMaPPer optical centers and the SZ (ACT DR5) or X-ray cluster centers, respectively. The solid blue and green lines show the best-fitting two-component 2D Gaussian models for the SZ and X-ray cases. The dashed orange and red lines show results from Ding et al. (2025): the orange line corresponds to all 186 optical–SZ cross-matched clusters, while the red line excludes clusters affected by astrophysical effects (e.g., ongoing mergers) or systematic biases in the HSC data and cluster-finding algorithm. The dashed gray line shows the X-ray–optical offset distribution from Oguri et al. (2018), based on HSC Wide S16A clusters. All offsets are shown in physical units of h−1 kpc.

Therefore, the offsets involving SZ or X-ray centers are significantly smaller than those based on weak-lensing-determined centers. This suggests that weak lensing is identifying a different characteristic center–one that traces the projected gravitational potential and is influenced by the full matter distribution along the line of sight weighted by lensing efficiency.

Unlike X-ray or SZ signals, which trace the distribution of hot baryonic gas through electron density and temperature, weak lensing directly probes the total mass, including both baryonic and DM, and is independent of the dynamical state of the cluster. Consequently, in unrelaxed systems or in clusters embedded in complex environments with substantial line-of-sight structures, weak-lensing-derived centers may appear significantly offset from the baryon-based tracers. These intrinsic differences in physical sensitivity and projection response among tracers likely contribute to the broader offset distributions observed in WL-based analyses.

For clusters with the largest centroid offsets, the projected distributions of member galaxies are found to be frequently bimodal, suggesting ongoing mergers or significant dynamical disturbances; however, given the current data quality and the unknown dynamical information of the clusters, we refrain from definitive classification and defer a more detailed analysis to future work.

4.4. Distribution of halo ellipticities

Ellipticity is a key parameter used to characterize the projected shape of DM halos along the line of sight, and its importance in understanding halo structure is well established. Although the spherically symmetric NFW model successfully reproduces the spherical average of the halo density profile, both N-body simulations (Jing & Suto 2002) and observations (Oguri et al. 2010; Liu et al. 2025) consistently show that halos are intrinsically nonspherical.

Moreover, halo shape provides insights into the nature of DM. For instance, halos in self-interacting dark matter (SIDM) models are predicted to be more spherical than those in CDM scenarios, as particle self-interactions tend to isotropize the DM distribution (Gonzalez et al. 2024).

In our analysis, we performed an individual shear map fitting for each cluster to estimate its own ellipticity, and then combined the resulting values to construct the overall probability distribution, as shown in Fig. 8. This approach captures the intrinsic variation in halo shapes across the sample. Errors were estimated using 1000 bootstrap resamplings. The measured ellipticity for the full sample is e = 0.530 ± 0.168, and the mean ellipticity is ⟨e⟩ = 0.505 ± 0.007, which is in good agreement with predictions from N-body simulations (Jing & Suto 2002) and previous observational results (Liu et al. 2025).

Fig. 8. Upper panel: Probability distribution function (PDF) of the ellipticity values derived from individual shear map fittings of each cluster. Blue points represent the observed distribution, with error bars estimated from 1000 bootstrap resamplings. The red curve shows a Gaussian fit to the distribution, centered at e = 0.530 ± 0.168. For comparison, the dashed black line represents the ellipticity distribution reported by Liu et al. (2025), converted to match our ellipticity definition. Lower panel: PDF of the sample-averaged ellipticity, ⟨e⟩, computed from the full cluster ensemble. The red curve shows a Gaussian fit centered at ⟨e⟩ = 0.505 ± 0.007.

In Table 2, we divided the sample into four subsets in the mass–redshift plane based on the mean halo mass and mean redshift. Specifically, clusters are grouped into quadrants: (1) low mass & low redshift, (2) low mass & high redshift, (3) high mass & low redshift, and (4) high mass & high redshift. We computed the average ellipticity in each bin. The results show that halos with higher mass and lower redshift tend to exhibit slightly lower ellipticities, implying that they appear more spherical along the line of sight. This trend contrasts with predictions from DM-only cosmological simulations (Jing & Suto 2002; Lau et al. 2021), which include only gravitational interactions. These simulations typically find that more massive and lower-redshift halos are slightly more elliptical due to persistent anisotropic mass accretion from the surrounding cosmic web. Interestingly, our findings align more closely with the trend for halo masses larger than 10^13.5 h⁻¹ M_⊙ in the hydrodynamical IllustrisTNG simulation (Liu et al. 2025), which incorporates baryonic physics. In this context, more massive halos are expected to assemble earlier and more rapidly, reducing subsequent interactions with their large-scale environments and thus evolving into more isotropic and spherical configurations.

5. Conclusion and discussion

In this paper, we have analyzed 299 individual redMaPPer-selected clusters by constructing their 2D shear maps from the HSC first-year shear catalog with high data quality. We statistically investigated key halo properties, including concentration, ellipticity, and the BCG offset. To mitigate the impact of intrinsic mass uncertainties on the statistical concentration–mass(c–M) relation, we adopted a mass–richness prior from Rykoff et al. (2012) in our fitting procedure. Furthermore, we used the inverse of the posterior variance of each fitted parameter as the statistical weight when computing weighted average values. Our main conclusions are as follows:

• We find that halo concentration decreases with increasing mass in our sample. The derived c–M relation is steeper than those obtained from weak lensing analyses but more consistent with trends seen in strong-lensing-selected cluster studies. This may reflect that our sample selection and fitting method are more sensitive to halo concentration at the high-mass end. In the mass–redshift plane, we observe that halos with lower redshift and higher mass tend to have lower concentrations.

• The distribution of offsets between the observed optical center and the weak-lensing-determined center is well described by a two-component 2D Gaussian model, as defined in Eq. (28). The characteristic offset scales are ${\sigma }_1=115.2_{-12}^{+9.9}{h}^{-1}$ kpc and ${\sigma }_2={409}_{-38}^{+27}{h}^{-1}$ kpc, corresponding to the well-centered and miscentered populations, respectively. The fraction of well-centered clusters is measured as f_cen = 0.43 ± 0.05. We find that halos in the lower mass bin tend to exhibit larger offsets, possibly due to their greater susceptibility to perturbations from surrounding large-scale structures and enhanced lensing projection effects. In comparison, lower-redshift clusters exhibit smaller offsets and higher centering fractions, consistent with the expectation that clusters become more relaxed over cosmic time. The weak-lensing-based offsets are systematically larger than those involving X-ray or SZ-defined centers (Oguri et al. 2018; Ding et al. 2025), likely reflecting differences in tracer sensitivity and center definitions. In particular, weak lensing probes the total projected mass distribution – including DM – while X-ray and SZ observations trace the hot baryonic gas, which is more tightly correlated with the optical BCG position. These intrinsic physical differences among tracers naturally lead to broader offset distributions in WL-based measurements. Our results demonstrate that modeling miscentering with a flexible two-component model significantly improves the reliability of halo property inference. This is crucial for reducing systematics in GGL analyses and for ensuring the robustness of high-precision cosmological probes, such as cluster abundance studies (Ghirardini et al. 2024) and three-by-two-point (3 × 2 pt) analyses (Heymans et al. 2021; Li et al. 2023b; DES Collaboration 2025), where even small biases in halo centering can impact cosmological parameter constraints. A more detailed treatment of these effects will be presented in a forthcoming paper.

• The ellipticity distribution derived from individual cluster fittings has a mean value of e = 0.530 ± 0.168, and the sample-averaged ellipticity is ⟨e⟩ = 0.505 ± 0.007. In the mass–redshift plane, halos with higher mass and lower redshift appear more spherical, exhibiting slightly lower average ellipticities. This trend is opposite to predictions from dark-matter-only N-body simulations but is more closely consistent with the trend for halo masses larger than 10^13.5 h⁻¹ M_⊙ from the IllustrisTNG simulations (Liu et al. 2025) that include baryonic effects. In those models, baryonic cooling contracts matter toward the center, increasing halo roundness (Lau et al. 2012), while isotropic momentum redistribution from merger events also promotes spherical shapes. A similar effect is observed in SIDM models, where particle collisions isotropize the halo structure (Gonzalez et al. 2024). Although our current measurements cannot distinguish between these two physical mechanisms, the results offer meaningful constraints on the nature of DM.

Overall, our results demonstrate that 2D weak-lensing analysis is a powerful tool for investigating the structural properties of individual galaxy clusters, and it holds significant potential for improving the treatment of systematic uncertainties in traditional stacking methods. Future wide-field surveys, such as the Vera C. Rubin Observatory, the Euclid mission, and CSST, will deliver larger samples of clusters and background galaxies, enabling more precise constraints on halo properties and offering new opportunities to probe the fundamental nature of DM.

Acknowledgments

This research is supported by National Key R&D Program of China No. 2022YFF0503403. X.K.L. acknowledges the support from NSFC grant No. 12173033, the grants from the China Manned Space Projects with No. CMS-CSST-2021-B01, the “Yunnan Key Laboratory of Survey Science” with project No. 202449CE340002 and a “Yunnan Provincial Top Team Projects” with project No. 202305AT350002. H.Y.S. acknowledges the support from the Ministry of Science and Technology of China (grant No. 2020SKA0110100), Key Research Program of Frontier Sciences, CAS, Grant No. ZDBS-LY-7013 and Program of Shanghai Academic/Technology Research Leader. Z.H.F. acknowledges the support from NSFC grant No. 11933002, and the grant from the China Manned Space Projects with No. CMS-CSST-2021-A01. We also acknowledge the support from the science research grants from the China Manned Space Project with No. CMS-CSST-2025-A03 and CMS-CSST-2025-A05. Part of the computations in this study were carried out on the servers of the SWIFAR Cosmology Group and the Yunnan University Astronomy Supercomputer. The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from the Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University. This paper makes use of software developed for Vera C. Rubin Observatory. We thank the Rubin Observatory for making their code available as free software at http://pipelines.lsst.io. This paper is based on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by the Subaru Telescope and Astronomy Data Center (ADC) at NAOJ. Data analysis was in part carried out with the cooperation of Center for Computational Astrophysics (CfCA), NAOJ. We are honored and grateful for the opportunity of observing the Universe from Maunakea, which has the cultural, historical and natural significance in Hawaii. The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg, and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

References

All Tables

All Figures

	Fig. 1. Normalized distribution of spectroscopic redshifts for the lens sample (dashed purple) and photometric redshifts for the source sample (solid blue).
In the text

Fig. 2. Upper panels: The left panel shows the observed 2D shear map for an illustrative redMaPPer cluster (M ≃ 1014.73 h−1 M⊙), while the right panel displays the corresponding best-fitting shear field predicted by the elliptical NFW model. Each map covers an area of 18′×18′. In both maps, the orientation and length of the sticks within each 1′×′ pixel indicate the local tangential direction and amplitude of the distortion, respectively. For visualization purposes, both shear fields have been smoothed with a Gaussian kernel with a full width at half maximum (FWHM) of ≈2′ (Oguri et al. 2010). The blue dot marks the weak-lensing-derived center, and the red pentagram denotes the observed center. Lower panel: Posterior constraints on structural parameters of the cluster, including the halo concentration, ellipticity, and projected offsets of the mass centroid. The contours represent the 1σ and 2σ confidence levels.

In the text

Fig. 3. Concentration–mass (c–M) relation. Light gray points represent the best-fitting concentration and mass for individual clusters, while black points show the weighted average values in each mass bin, with error bars indicating the standard deviation. For comparison, the results for strong-lensing selected halos from Meneghetti et al. (2014) are shown by a solid orange line. Weak-lensing-based results from Xu et al. (2021) (upturn model) are plotted as a dashed-dotted blue line, and the power-law model from Duffy et al. (2008) derived from DM-only numerical simulations, is plotted as a dashed red line, respectively.

In the text

Fig. 4. Concentration–mass–redshift (c–M–z) relation. Gray points indicate the positions of individual clusters in the mass–redshift plane. Colored diamonds represent the weighted mean concentration within each bin. Bins containing fewer than ten clusters are excluded from the analysis. Colored contour lines trace surfaces of constant concentration across the plane, and the color bar (labeled “concentration”) gives their numerical values.

In the text

Fig. 5. Upper panel: Probability distribution function (PDF) of the projected offset between the BCG position and the weak-lensing-determined halo center. Red data points show the observed BCG offsets in our sample, with error bars estimated via 1000 bootstrap resamplings. The solid blue and black lines represent the best-fitting single- and two-component 2D Gaussian models, respectively. The dashed purple and orange lines correspond to the well-centered and miscentered components of the two-component model. All offset scale parameters, including σ, σ1, and σ2, are expressed in units of h−1 kpc. Lower panels: Fitting results for the two-component (left) and single-component (right) 2D Gaussian models. Contours indicate 1σ and 2σ confidence levels, with best-fitting values and 68% confidence intervals shown above each 1D marginalized distribution.

In the text

Fig. 6. Probability distribution functions (PDFs) of the projected offsets between the BCG and the weak-lensing-determined halo center for different redshift and mass bins. All error bars are estimated from 1000 bootstrap resamplings. Red and blue points represent the low- and high-redshift subsamples, respectively, with solid lines showing the best-fitting two-component 2D Gaussian models. Corresponding vertical lines indicate the values of σ1 and σ2. Pink and light blue points show the low- and high-mass subsamples, respectively, with dashed lines and vertical markers indicating the best fits and offset scales.

In the text

Fig. 7. Probability distribution functions (PDFs) of the projected offsets between the BCG position and the SZ- or X-ray-determined cluster centers. Black and purple data points represent the offsets between the redMaPPer optical centers and the SZ (ACT DR5) or X-ray cluster centers, respectively. The solid blue and green lines show the best-fitting two-component 2D Gaussian models for the SZ and X-ray cases. The dashed orange and red lines show results from Ding et al. (2025): the orange line corresponds to all 186 optical–SZ cross-matched clusters, while the red line excludes clusters affected by astrophysical effects (e.g., ongoing mergers) or systematic biases in the HSC data and cluster-finding algorithm. The dashed gray line shows the X-ray–optical offset distribution from Oguri et al. (2018), based on HSC Wide S16A clusters. All offsets are shown in physical units of h−1 kpc.

In the text

Fig. 8. Upper panel: Probability distribution function (PDF) of the ellipticity values derived from individual shear map fittings of each cluster. Blue points represent the observed distribution, with error bars estimated from 1000 bootstrap resamplings. The red curve shows a Gaussian fit to the distribution, centered at e = 0.530 ± 0.168. For comparison, the dashed black line represents the ellipticity distribution reported by Liu et al. (2025), converted to match our ellipticity definition. Lower panel: PDF of the sample-averaged ellipticity, ⟨e⟩, computed from the full cluster ensemble. The red curve shows a Gaussian fit centered at ⟨e⟩ = 0.505 ± 0.007.

In the text