# Strong-lensing and kinematic analysis of CASSOWARY 31: Can strong lensing constrain the masses of multi-plane lenses?

H. Wang<sup>1,2,*</sup>①, R. Cañamas<sup>1,2,3</sup>①, S. H. Suyu<sup>1,2,4</sup>, A. Galan<sup>2,1</sup>①, C. Grillo<sup>5,6</sup>①,
G. B. Caminha<sup>2,1</sup>①, and L. Christensen<sup>7</sup>

Received 8 March 2024 / Accepted 2 July 2024

## ABSTRACT

We present a mass measurement for the secondary lens along the line of sight (LoS) from the multi-plane strong lens modeling of the group-scale lens CASSOWARY 31 (CSWA 31). The secondary lens at redshift $z = 1.49$ is a spiral galaxy well aligned along the LoS with the main lens at $z = 0.683$. Using the MUSE integral-field spectroscopy of this spiral galaxy, we measured its rotation velocities and determined the mass from the gas kinematics. We compared the mass estimation of the secondary lens from the lensing models to the mass measurement from kinematics, finding that the predictions from strong lensing tend to be higher. By introducing an additional lens plane at $z = 1.36$ for an overdensity known to be present, we find a mass of $\approx 10^{10} M_{\odot}$ enclosed within 3.3 kpc of the centroid of the spiral galaxy, which approaches the estimate from kinematics. This shows that secondary-lens mass measurements from multiple-plane modeling are affected by systematic uncertainties from the degeneracies between lens planes and the complex LoS structure. Conducting a detailed analysis of the LoS structures is therefore essential to improve the mass measurement of the secondary lens.

**Key words.** gravitational lensing: strong – Galaxy: kinematics and dynamics – galaxies: high-redshift – galaxies: spiral

## 1. Introduction

In the Λ cold dark matter model, galaxies grow hierarchically via the accretion of cold gas from the circumgalactic medium and mergers with nearby galaxies. Measuring robust galaxy masses over broad redshift ranges is key to improving our understanding of this broad picture (e.g., Peng et al. 2010). Yet, at high redshifts ($z \ge 1$), direct measurements are challenging due to observational limitations, and mass estimates not only rely on assumptions derived from our knowledge of the nearby Universe, but are also prone to several systematics. For instance, for gas-rich spiral galaxies, dynamical masses inferred from gas kinematics can be significantly biased in the case of departures from smooth rotation due to secular instabilities (e.g., Yu et al. 2022), merging events, or feedback processes (e.g., Soto et al. 2012).

In this context, strong lensing has been established as a robust, accurate tracer of the total mass and inner structure of foreground galaxies (e.g., Koopmans et al. 2006; Auger et al. 2010; Shajib et al. 2022), galaxy groups, and clusters (e.g., Treu 2010; Meneghetti et al. 2017). Photons traveling from distant sources get deflected by an overdense region associated with the foreground galaxies and form arcs or a full Einstein ring around the deflectors. The locations of the lensed images and the surface

brightness of the extended arcs can be used to constrain the mass distribution of the foreground deflectors.

In most strong lensing scenarios, a single overdense region along the path between the source and observer is considered for modeling. The line-of-sight (LoS) structure is neglected, which is usually a sensible approximation since the main lens contributes dominantly to the deflection angle (e.g., Gavazzi et al. 2008). In some cases, galaxies aligned with the main lens can nonetheless produce non-negligible perturbations, particularly for group- and cluster-scale lenses, which tend to have rich LoS structures (Bayliss et al. 2014). To reconstruct the lensed sources and improve the root mean square (rms) of the predicted image position for such systems, a secondary lens needs to be included when computing multiple-lens-plane models (e.g., Schneider 2019; Wong et al. 2017).

In this work we further explored the constraints on the mass of secondary lenses inferred from multi-plane lens models by conducting a case study based on the group-scale lens CASSOWARY 31 (CSWA 31; Belokurov et al. 2009; Stark et al. 2013). In Wang et al. (2022), through detailed modeling of the multi-scale mass distribution, we found that this system comprises a spiral galaxy at $z = 1.49$ that is well aligned with the main lens at $z = 0.683$, with both structures contributing significantly to the deflection angles of source galaxies at higher redshifts. In addition to strong lensing, gas kinematics inferred

<footnote n="1">1 Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany</footnote>

<footnote n="2">2 Technical University of Munich, TUM School of Natural Sciences, Physics Department, James-Franck Str. 1, 85748 Garching, Germany</footnote>

<footnote n="3">3 Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="4">4 Academia Sinica Institute of Astronomy and Astrophysics (ASIAA), 11F of ASMAB, No. 1, Section 4, Roosevelt Road, Taipei 10617, Taiwan</footnote>

<footnote n="5">5 Dipartimento di Fisica, Università degli Studi di Milano, via Celoria 16, I-20133 Milano, Italy</footnote>

<footnote n="6">6 INAF-IASF Milano, via A. Corti 12, I-20133 Milano, Italy</footnote>

<footnote n="7">7 Cosmic Dawn Center, Niels Bohr Institute, Univ. of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark</footnote>

<footnote>* Corresponding author; wanghan@mpa-garching.mpg.de</footnote>
from integral-field spectroscopy (IFS) with the Multi-Unit Spectroscopic Explorer (MUSE; Bacon et al. 2015) allowed us to infer an independent mass estimate for this secondary lens. While such known systems with two or more deflectors are rare, their number is expected to increase with the next generation of deep, wide-field imaging surveys such as the Legacy Survey of Space and Time, *Euclid*, and *Roman* (e.g., Collett 2015; Mandelbaum et al. 2018). It is therefore becoming essential to investigate multi-plane lens modeling as an alternative mass probe to study the inner structure of galaxies acting as secondary deflectors.

The paper is organized as follows. In Sect. 2 we briefly introduce the properties of CSWA 31 and the datasets. In Sect. 3 we present the independent mass measurements from strong lensing and kinematic modeling. In Sect. 4 we discuss the uncertainties of each method and compare the results. We also present an updated model with an additional lens plane. In Sect. 5 we summarize our findings and discuss the validity of mass predictions for secondary lenses obtained from multi-plane lens modeling. Throughout this work we assume $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_m = 0.3$, and $\Omega_\Lambda = 0.7$. Hence, $1''$ corresponds to $7.08 \text{ kpc}$ at the main lens redshift of $z = 0.683$, and to $8.46 \text{ kpc}$ at the secondary lens redshift of $z = 1.49$.

## 2. Data

CSWA 31 is a group-scale lens with a central ultra-massive brightest group galaxy (BGG) at $z = 0.683$, surrounded by six sets of multiple images from five independent background sources. The high-resolution *Hubble* Space Telescope (HST) imaging in *F160W* and *F438W* filters, the VLT/MUSE spectroscopy of the field (program 0104.A-0830(A), PI: Cañameras), and the joint analysis of these data are described in detail in Wang et al. (2022). The background spiral galaxy at $z = 1.49$ studied in this letter and referred to as S0, is lensed into four multiple images forming an extended arc (see Fig. 1). We identified five additional sets of multiple images – S1, S2, S3, S4, and S5 – at four independent redshifts using MUSE spectroscopic and HST (*F160W*-*F438W*) colors. Image sets S1 and S2 are at redshift $z = 1.49$ and S3, S4 and S5 at redshifts $z = 2.76, 3.43$, and $z = 4.21$, respectively. In addition, 46 group members located within $\pm 3000 \text{ km s}^{-1}$ of the BGG redshift were selected from the MUSE spectroscopic catalog.

The spiral galaxy S0 displays a disk-like morphology in HST/*F160W*, with a bright bulge and three main spiral arms, and is apparently oriented nearly face-on. The best-fit Sérsic profile of image S0(d) shows that this galaxy has an effective radius $R_{\text{eff}}$ of $1.17''$, or about $3.8 \text{ kpc}$ in the source plane (after correcting for magnification), and an axis ratio $q = 0.85$. The properties of S0 are similar to the strongly lensed spiral galaxy at $z \approx 1.5$ hosting SN Refsdal, for which Di Teodoro et al. (2018) inferred 3D kinematic models.

In this work, we focus on the least distorted image S0(d), shown in Fig. 1 (left inset), to conduct a kinematic analysis of S0 in the image plane (see also Di Teodoro et al. 2018). We verify that spatial variations of lensing shear are negligible for S0(d) using the mass model Esr2-MPtest from Wang et al. (2022) (baseline model hereafter). We use this model to map small circles of $0.052''$ radius distributed over the source plane at $z = 1.49$, to the image S0(d). The circular sources get magnified and distorted into ellipses in the image plane, but these ellipses do not vary significantly in size and orientation over S0(d), suggesting overall uniform lensing distortions over this image (see the green ellipses on the inset of Fig. 1 left panel). In addition,

the magnification factor $\mu \approx 7$ of image S0(d) offers a good opportunity to zoom in onto the gas kinematics of this spiral.

The MUSE data cube has a point spread function (PSF) with a full width at half maximum (FWHM) of $0.75'' \times 0.73''$. Its wavelength resolution is FWHM $= 2.35 \AA$, with a channel width of $1.25 \AA$. To focus on the [OII] doublet, the cube was cropped between $3719 \AA$ to $3734 \AA$ in the rest frame of S0. The [OII] line intensity map of S0(d) shows an asymmetrical feature, indicating that ionized gas reservoirs form denser regions along the spiral arms seen in stellar continuum with HST (see Fig. 1). The total [OII] flux over the source is $203.8 \times 10^{-18} \text{ erg s}^{-1}$. The LoS velocity map is derived by fitting double Gaussian functions to each pixel, representing a preliminary step to investigate the presence of velocity gradients within the S0(d) image. Additionally, these maps offer initial values for geometrical parameters in $^3\text{Barolo}$.

## 3. Methodology

To investigate the constraints on secondary lens masses, we inferred two independent mass estimates of the S0 spiral. In Sect. 3.1 we present the kinematic modeling procedure and the mass measurement from the best-fit rotation curve. In Sect. 3.2 we present the multi-plane strong lens modeling and the S0 mass measurement from the lensing approach.

### 3.1. Mass measurement from kinematic modeling

Spiral galaxies at $z \ge 1$ are thought to settle into marginally stable, rotating disks with turbulence (e.g., Bournaud et al. 2009; Nestor Shachar et al. 2023). The apparent velocity gradient in Fig. 1 and the smooth disk-like morphology from HST suggest a rotating velocity field. We adopted $^3\text{Barolo}$ (Di Teodoro & Fraternali 2015) to derive a simple rotating disk model of S0 using the MUSE data of image S0(d), following a similar approach as Di Teodoro et al. (2018). This image is well de-blended from other galaxies in the field and the minor distortions of S0(d) caused by strong lensing vary on scales smaller than the MUSE PSF. The beam-smearing of galaxy kinematics, which is accounted for by $^3\text{Barolo}$, dominates over systematic errors due to lensing shear.

$^3\text{Barolo}$ initially identifies a mask region within the MUSE data cube encompassing the emission lines. Subsequently, it builds a kinematic model convolved with the MUSE PSF to characterize these lines by using tilted rings. In our case, $^3\text{Barolo}$ replicates the [OII] emission doublet at $\lambda = 3726 \AA$ and $3729 \AA$. Each ring is defined with a centroid, position angle, inclination angle *i*, systematic velocity $V_{\text{sys}}$, rotation velocity $V_{\text{rot}}$ and velocity dispersion $\sigma_v$. We built the kinematic model by using four tilted rings, with a $0.4''$ width, and adopting the S0(d) centroid in HST imaging with the S0 position angle determined from the velocity gradient (see Fig. 1) as initial values. We fixed the $V_{\text{sys}}$ in all tilted rings to $0 \text{ km s}^{-1}$. We iteratively refined these values with $^3\text{Barolo}$ and subsequently allowed variations in only the $V_{\text{rot}}$ and $\sigma_v$ at each ring (see Fig. A.1). In the fitting, we adopted the local normalization and a uniform weighting function<fnref n="1" /> to account for the asymmetry in the intensity map of S0(d).

$^3\text{Barolo}$ projects the tilted ring model to the plane of the sky using *i* and compares it to the observational data. The inclination angle is an important parameter, but it is hard to estimate

<footnote n="1">¹ For the “type-two” residual defined in Eqs. (2) and (3b) in Sect. 2.4 of Di Teodoro & Fraternali (2015).</footnote>
<figure><img src="image_1.png" /><figcaption>Fig. 1. HST imaging and kinematics of the lensed galaxy S0. Left panel: HST/WFC3 image of CSWA 31 in the F160W band. The orientation is marked by the black arrows (north upward and east leftward) and kept fixed for other panels. The BGG at z = 0.683 is surrounded by six multiple-image sets at different redshifts marked with circles, using a different color for each set. The image sets S0 and S1 are bright features from the same spiral galaxy at z = 1.49, with a zoomed-in image of S0(d) in the right inset. The green square shows the de-lensed S0 position in the HST imaging. The green ellipses in the zoomed-in panel illustrate the minor variations in distortion over S0(d) caused by strong lensing. Right panel: [OII] intensity map (top) and the LoS velocity map (bottom) for the S0(d) image extracted from MUSE IFS. The FWHM of the MUSE PSF is shown in the upper left in the intensity map. The bold black cross shows the adopted S0(d) centroid determined from 3D Barolo.</figcaption></figure>

in practice due to the unknown intrinsic shape of galaxies. Without losing any generality, we built a range of kinematic models covering i from 5° to 35° (see the discussion in Sect. 4.1). The inclination is not well constrained, and Fig. 2 shows a typical best-fit model per wavelength channel, for i = 20.8°.

Due to the moderate spectral resolution of the MUSE data cube, velocity dispersion values returned by 3D Barolo decrease to 0 km s⁻¹ for the inner radius of 0.2'' and the outer radius of 1.4'', for all kinematic models with inclination angles between 5° and 35°. To address this limitation, we reduced the number of rings to three and obtained the best-fit σv at distances of 0.6'', 0.9'', and 1.2'' from the centroid of image S0(d), finding an average value of σv = 41.3 km s⁻¹. However, the residuals in the channel maps are worse than the model built with four rings (see Fig. A.2). We thus decided to retain the model with four rings and solely adopted the Vrot along the major axis at the second ring at r = 0.6'' and the third ring at r = 1.0'' to determine the dynamical mass of S0 (see Fig. 3).

We determined the total mass of S0 using Vrot via

$$M_{\text{rot}}(r') = \frac{V_{\text{rot}}^2}{G} r', \qquad (1)$$

where r' is the radial distance from the S0 centroid in the source plane at z = 1.49. Strong lensing magnifies a circular source with area πr'^2 by a factor of μ in the image plane (see Table 1). We converted the radii r in the S0(d) image to radii r' in the source plane via

$$r' = \frac{r}{\sqrt{\mu}}. \qquad (2)$$

We obtained a total mass of S0 of ≈ 10¹⁰ M☉ at 3.3 kpc in the source plane, which is about the effective radius of S0 at z = 1.49 (see Fig. 4).

## 3.2. Mass measurement from multi-plane lens modeling

We summarize in Table 1 the different lens models we considered in this work. In this section, we adopt the multi-plane lens models Img-2MP, Img-2MP(L/D) and baseline from Wang et al. (2022). The average magnification factor, μ, of S0(d) has minor variations among these models, and the induced shear is negligible, as shown in Fig. 1 for the baseline model.

These three models consist of two lens planes, and were obtained with the GLEE software (Suyu & Halkola 2010; Suyu et al. 2012). The main lens plane at z = 0.683 incorporates the mass contributions from the BGG, 46 group members confirmed by MUSE spectroscopy, the group halo, and a constant external shear to account for the dense environment surrounding CSWA 31. The mass of the BGG and group members are modeled with truncated dual pseudo-isothermal elliptical (dPIE; Elíasdóttir et al. 2007; Suyu & Halkola 2010) mass profiles. Due to the limited number of constraints from lensing, the Einstein radii and truncation radii of group members are scaled with their luminosities in the F160W band following the approach in Grillo et al. (2015) and Chirivì et al. (2018). The underlying group halo is depicted by a softened power-law elliptical mass profile with an extended core radius (Barkana 1998). The secondary lens plane at z = 1.49 accounts for the mass of the spiral galaxy S0, which serves as a lensed source, being deflected by the main lens. Additionally, as a secondary lens along the LoS,
<figure><img src="image_2.png" /><figcaption>Fig. 2. Comparison between MUSE IFS of image S0(d) (blue contours) and the 3D Barolo model (red contours). The third row in the upper and bottom panels shows the residuals. Contour levels are set at 2σn, where the σ noise = 2 × 10⁻²⁰ erg s⁻¹ cm⁻² Å⁻¹ and n = 1, 2, 4, 8. The green cross shows the centroid position of S0(d). This plot shows the fitting result for the two lines of the [OII] emission doublet at (rest frame) λ = 3726 Å and 3729 Å in the upper and bottom panels, respectively.</figcaption></figure>

S0 deflects the light from the background sources S3, S4, and S5 at higher redshifts. Accounting for this secondary strong deflec- tor (S0) improves the fitting of the image positions of S3, S4, and S5, while also enabling the simultaneous reconstruction of the extended sources S0 and S3 in the baseline model. Excluding the secondary lens plane in the modeling results in the failure to reconstruct S3 and leads to higher errors in the predicted image positions of S3, S4, and S5.

In all three models, S0 is modeled with a spherical isothermal (SIS) profile, which has a 3D density profile described as

$$
\rho(r) = \frac{\sigma_v^2}{2\pi Gr^2}, \qquad (3)
$$

where σ_v is the velocity dispersion, related to the Einstein radius θ_E (for a lensed background source at infinite redshift) via σ_v² = θ_E c² / (4π). Thus, Eq. (3) can be written in terms of θ_E:

$$
\rho(r) = \frac{\theta_E c^2}{8\pi^2 Gr^2}. \qquad (4)
$$

We calculated the mass of S0 enclosed within radius r' in spher-
ical coordinates,

$$
M_{\text{lens}}(r') = \int_{0}^{r'} 4\pi r^2 \frac{\theta_E c^2}{8\pi^2 G r^2} dr, \qquad (5)
$$

such that

$$
M_{\text{lens}}(r') = \frac{\theta_E c^2}{2\pi G} r', \qquad (6)
$$

<figure><img src="image_3.png" /><figcaption>Fig. 3. Histogram of V<sub>rot</sub> weighted by the BIC. The upper and bottom panels show the weighted V<sub>rot</sub> distribution at 2.0 kpc and at 3.3 kpc from the galaxy centroid in the source plane, respectively. The 1σ uncertainties on V<sub>rot</sub> are determined by drawing 5000 samples of the inclination angles in the range [5°, 35°].</figcaption></figure>

where r' is the radial distance in the source plane. Equation (6)
is valid in the inner region of the galaxy within the halo
truncation.

The lens models Img-2MP, Img-2MP(L/D), and baseline
adopted in this study rely on different sets of observational con-
straints to ascertain the best-fit values of the Einstein radius θ<sub>E,S0</sub>
(see Table 1). The model Img-2MP adopts image positions of
six multiple-image sets as constraints. In model Img-2MP(L/D),
Wang et al. (2022) combines the stellar kinematics of the BGG
with the multiple image positions, to perform a joint strong
lensing and dynamical modeling. The baseline model uses
extended arcs of S0 and S3 in addition to the image positions
of three multiple-image sets S2, S4, and S5 as constraints. In
Fig. 4 we plot the S0 mass distribution as a function of the
source plane radius for these three models. The S0 mass pre-
dicted from the baseline model is three times smaller than
that for the other two, and we further analyze this discrepancy
in Sect. 4.2.

4. Mass comparison and discussion

Strong lens models presented in Sect. 3.2 tend to assign signif-
icantly more mass to S0 than the modeling of gas kinematics.
We discuss this discrepancy by highlighting the main uncertain-
<table><caption>**Table 1.** Strong lens models used in this study, their different constraints, the best-fit values of $\theta_{E,S0}$, and the average magnification factors, $\mu$, weighted by the surface brightness over the S0(d) image.</caption><thead><tr><th>Model name</th><th>Constraints</th><th>$\theta_{E,S0}$</th><th>$\mu$</th></tr></thead><tbody><tr><td>Img-2MP</td><td>6 lensed-image sets</td><td>4.4"</td><td>5.5</td></tr><tr><td>Img-2MP (L/D)</td><td>6 sets + kinematics of BGG</td><td>4.6"</td><td>6.9</td></tr><tr><td>baseline</td><td>3 sets + extended arcs of S0, S3</td><td>1.65"</td><td>6.7</td></tr><tr><td>Esr2-3MP</td><td>3 sets + extended arcs of S0, S3</td><td>0.85"</td><td>7.3</td></tr></tbody></table>

**Notes.** The first three models share the same configuration and mass parametrization. In the last model Esr2-3MP, we accounted for the additional perturbation caused by foreground galaxies at $z = 1.36$ and performed lens modeling with three lens planes.

ties from each approach in Sects. 4.1 and 4.2. In Sect. 4.2.1 we discuss the results from the two-plane lens models that are presented in Wang et al. (2022). Then, in Sect. 4.2.2, we present an additional model with a third lens plane to further explore the impact of LoS structures.

## 4.1. Uncertainties from the kinematic modeling

The unknown inclination angle $i$ of the spiral has the largest contribution to the error budget on the mass inferred from kinematic modeling. Since S0 is observed nearly face-on, a small variation of $i$ leads to significant changes in the rotation velocity along the major axis, as $V_{rot} = V_{LoS}/\sin i$. To obtain a sensible estimation of the uncertainty caused by $i$, we randomly drew 5000 samples from a broad range of angles $i$ over [5°, 35°] and fit the MUSE data cube for each sample. We performed a model comparison by making use of the Bayesian information criterion (BIC). Since we used the same data cube and only allowed the $V_{rot}$ and $\sigma_v$ to vary in each sample, the BIC in our case is computed by the log-likelihood term. We followed the approach of Birrer et al. (2019) defining the BIC weight $w_i$ with respect to the minimal BIC value $BIC_{min}$ in our sample,

$$w_i = \frac{\exp\{-(BIC_i - BIC_{\text{min}})/2\}}{\sum_{i=1}^{5000} \exp\{-(BIC_i - BIC_{\text{min}})/2\}} \quad (7)$$

We then determined the median $V_{rot}$ and the $1\sigma$ uncertainty from the 16th, 50th, and 84th of the weighted $V_{rot}$ at radii of 0.6" and 1.0" (on the image plane) from the centroid of image S0(d), which correspond to 2.0 kpc and 3.3 kpc in the source plane (see Fig. 3).

We obtained the median values of $V_{rot} = 123 \text{ km s}^{-1}$ and $\sigma_v = 40 \text{ km s}^{-1}$ at 3.3 kpc from the galaxy center, and a ratio $V_{rot}/\sigma_v$ of 3.1 showing that this spiral galaxy is indeed rotationally-supported. Median values are broadly consistent with the evolution of ionized gas kinematics of star-forming galaxies at $z \approx 1-2$ in the literature (e.g., Wisnioski et al. 2015; Übler et al. 2019; Genzel et al. 2020; Price et al. 2021; Puglisi et al. 2023). With a similar modeling of the [OII] doublet, Di Teodoro et al. (2018) characterized the kinematics of the spiral hosting SN Refsdal up to 8 kpc from the galaxy center, probing into the flat portion of the rotation curve. While their similar $V_{rot}$ and lower $\sigma_v \approx 27 \text{ km s}^{-1}$ indicate stronger rotation support, the properties of both spirals at $z \approx 1.5$ are broadly consistent.

We used the value of $V_{rot}$ at 3.3 kpc to place S0 on the stellar-mass Tully-Fisher relation (TFR). We obtained an estimate of the galaxy stellar mass by fitting the g-, r-, i-, and z-band photometry of image S0(d) from the DESI Legacy Surveys DR10, together

<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Comparison of mass distributions for the spiral galaxy S0 at source plane from strong lensing and kinematic modeling. The solid lines show the S0 mass estimated from lens models using Eq. (6), given different datasets as constraints (see Table 1). The statistical $1\sigma$ uncertainty of the S0 mass from lens models is estimated from the 16th, 50th, and 84th percentiles of the posterior probability distribution functions (PDFs). The $1\sigma$ uncertainty along the x-axis of the S0 mass is determined from the magnification modification using Eq. (2). The gray squares show the S0 mass estimated from kinematic modeling using Eq. (1). The radial distance of the gray squares is divided by the $\mu$ from the baseline model. The error bar shows $1\sigma$ uncertainties associated with the range of inclination angles.</figcaption></figure>

with the F160W band from HST with CIGALE (Boquien et al. 2019) and following the approach described in Wang et al. (2022). This yields $M_* = (2.1 \pm 0.4) \times 10^{10} M_{\odot}$, corrected by $\mu = 6.7$ from the baseline model, showing that S0 follows the stellar-mass TFR at $z \approx 0$ from McGaugh & Schombert (2015) (see Fig. A.3). The stellar mass is also comparable to the total mass at $R_{eff}$ from the kinematic analysis, showing that the baryonic component dominates the inner mass budget.

Moreover, we tested the impact of the mask region detected by 3D Barolo on the fitting results. We used different detection modes and noise thresholds to infer the mask region and to reproduce the fitting. The resulting variations in output velocities are minor and within the $1\sigma$ range above.

As stated in Sect. 3.1, we neglected the effect of pressure support to the modeling, assuming that the gas kinematics of S0 follow perfect circular motions. The low intrinsic velocity dispersion suggests that the ionized gas in S0 is less turbulent than in other disks at $z \gtrsim 1$ (e.g., Swinbank et al. 2011), and asymmetric drift corrections due to pressure support should only have minor impact on the measured $V_{rot}$. Other feedback processes could also produce radial motions for the stellar and gas components in S0. In particular, large-scale outflowing winds are a common property for star-forming galaxies at $z \gtrsim 1$ (e.g., Newman et al. 2012; Bradshaw et al. 2013), and can influence their internal mass distributions. However, converting the [OII] flux of image S0(d) to a star formation rate, using the Kennicutt (1998) relation and a Salpeter (1955) initial mass function, would place S0 near, or even slightly below the star formation main sequence at $z \approx 1.5$ (e.g., Popesso et al. 2023). In consequence, while error bars shown in Fig. 4 might be broadened by such feedback effects, they are unlikely to play a significant role on the dynamical mass measurement of S0.
## 4.2. Uncertainties from strong lens modeling

The uncertainties in the S0 mass measurement through the lensing approach arise from degeneracies between the primary and secondary lens planes, as well as from LoS effects.

### 4.2.1. Two-lens-plane modeling

The mass distribution of S0 from the two-plane lens models described in Sect. 3.2 can be classified into two groups. Firstly, the model Img-2MP based on image positions and the model Img-2MP (L/D) combining with the BGG stellar kinematics both predict an elevated mass of ≈10¹¹ M☉ at 3.3 kpc. This falls on the high-mass end of the RC41 sample at z ≈ 0.65-2.45 from Genzel et al. (2017). Secondly, the baseline model using image positions and the extended arcs as constraints predicts a lower θ<sub>E,S0</sub> of 1.65". Because the mass degeneracies between the main lens plane at z = 0.683, and S0 at z = 1.49, are partially broken by the additional constraints from extended arcs, θ<sub>E,S0</sub> decreases from ≈ 4" to 1.65".

To further probe this difference, we checked the impact of changing the prior on θ<sub>E,S0</sub>, from the broad flat prior between 0" and 5" assumed in the baseline model, to an informative prior motivated by the kinematic analysis. Assuming that the S0 mass from dynamical and lens modeling are identical, we obtain

$$ \theta_E = \frac{2\pi}{c^2} V_{\text{rot}}^2. \qquad (8) $$

In this case, the best-fit θ<sub>E,S0</sub> in the baseline model corresponds to a constant V<sub>rot</sub> of 338.3 km s<sup>-1</sup> as a function of radius, a value above the 1σ upper limit from ³DBarolo (see Fig. 3). On the other hand, the median rotation velocities correspond to θ<sub>E</sub> ≈ 0.2". We assigned a Gaussian prior G (0.2", 0.07") to θ<sub>E,S0</sub> and recomputed the baseline model. The best-fit θ<sub>E,S0</sub> only decreases mildly to 1.62".

Consequently, parameter degeneracies cannot entirely explain the tension between the two modeling approaches. Since the S0 mass is also affected by the choice of the mass parameterization, we also modeled S0 with the elliptical power law and dPIE profiles, finding that the resulting S0 mass remains similar to that of the baseline. We thus consider that some additional sources along the LoS discarded from the model might be projected onto the mass plane at z = 1.49 and artificially boosting the S0 mass estimate.

### 4.2.2. Three-lens-plane modeling

To further investigate the impact of LoS effects on the S0 mass measurement from strong lensing, we built a model with a third lens plane at z = 1.36. Using the MUSE redshift catalog from Wang et al. (2022), we find 21 galaxies at z = 0.7-0.85 and seven at z = 1.36. Two galaxies at z = 1.36 are within 15" of the BGG centroid and have F160W < 22 mag, while none of the galaxies at z = 0.7-0.85 reach F160W < 22 mag within this aperture. Together with the distance to the extended arcs, this suggests that galaxies at z = 1.36 are stronger LoS perturbers. We assigned an SIS profile at z = 1.36 and performed the lens modeling with three planes at z = 0.683 (main lens plane), z = 1.36 (perturbers), and z = 1.49 (S0), adopting the same constraints as in the baseline model. The best-fit model results in θ<sub>E,S0</sub> = 0.85" and θ<sub>E,z</sub> = 1.36 = 0.59" (Einstein radius of perturbers at z = 1.36).

The enclosed mass of S0 within 3.3 kpc from this new model (Esr2-3MP; see Table 1) is close to the 68% percentile from the

kinematic measurement (see the brown line in Fig. 4). We projected the observed image positions of S0 at z = 0.683 onto the newly added lens plane at z = 1.36, finding that the projected images of S0 are close to the best-fit perturber position (see Fig. 5). The closer the perturber at z = 1.36 is to the projected S0 image positions, the greater its best-fit mass. When the perturber moves away from the projected S0 image positions, the value of θ<sub>E,S0</sub> falls back to that of the baseline model. In the Esr2-3MP model, we assigned broad flat priors in the range [0", 5"] for θ<sub>E,S0</sub> and θ<sub>E,z=1.36</sub>, and obtain strong degeneracies between the mass profiles in the three lens planes (see Fig. B.1). The θ<sub>E,S0</sub> is negatively correlated with the Einstein radius of the group-halo θ<sub>E,GH</sub> and with θ<sub>E,z=1.36</sub>, and is positively correlated with the Einstein radius of the BGG θ<sub>E,BGG</sub>. Our model Esr2-3MP provides comparable fitting results as baseline. It improves the modeling of extended arcs for S0 and S3 by O(100) in the log likelihood<fnref n="2" />, which indicates the goodness-of-fit (Suyu et al. 2006) but in the meantime, the χ<fnref n="2" /> of image position modeling worsens for S2, S4, and S5 by comparable amounts.

To further explore the mass distribution at z = 1.36, we performed the three-lens-plane modeling with either a constant mass sheet (Esr2-3MP (MS)) or a constant mass sheet with an external shear (Esr2-3MP (MS + Shear)) replacing the SIS. The former shows similar fitting results as Esr2-3MP, while the latter improves both the image position and extended image modeling by O(100). Notably, Esr2-3MP (MS + Shear) significantly improves the image position fitting of the multiple-image set S4 at z = 4.21. However, this model has nearly negligible mass on the mass sheet at z = 1.36, resulting in a substantial S0 mass prediction with θ<sub>E,S0</sub> = 2.01" (see Table B.1); this translates to V<sub>rot</sub> = 382 km s<sup>-1</sup>, which is higher than the V<sub>rot</sub> inferred from the kinematic model. The critical curves for the lens plane at z = 1.36 in these three models almost completely overlap while the caustics exhibit similar shapes but are scaled differently (see Figs. B.2 and B.3). The rescaling of caustics between Esr2-3MP and Esr2-3MP (MS) can be explained by a generalized mass-sheet transformation in a multi-plane lens scenario (Schneider 2019) that leaves all observables but the time delay invariant, corresponding to a uniform isotropic scaling in each lens and source plane. It is also not surprising that the critical and caustic curves of Esr2-3MP (MS) and Esr2-3MP (MS + Shear) for lens plane at z = 1.36 are nearly identical because the Esr2-3MP (MS) yields a best-fit mass sheet with 0.03 and Esr2-3MP (MS + Shear) has essentially zero mass sheet and a minor external shear with the strength of 0.04. In the Esr2-3MP (MS + Shear) model, the best-fit position angles of the external shear at z = 1.36 and the main lens plane at z = 0.683 are 49.0° and -37.9° counterclockwise from the x-axis, respectively. These directions of mass concentrations from the shear are consistent with the overdensity of galaxies in the environment of CSWA 31 (see Fig. B.2).

These models can be considered as additional tests illustrating the differences in the predicted S0 mass when adding a third lens plane. Overall, the model Esr2-3MP provides comparable fitting results as baseline and is also consistent with the kinematics modeling. In this sense, a concentrated mass profile performs better than adding a global change in the third plane. Nevertheless, the value of θ<sub>E,z=1.36</sub> = 0.59" is lower than for typical dark-matter halos and the best-fit perturber centroid does not fall in the vicinity of HST-detected galaxies. This shows that the

<footnote n="2">² After marginalizing over the source intensity pixel parameters, i.e., the log likelihood of the lens mass parameters is the log evidence of the source intensity reconstruction given the lens mass parameters.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Positions of the seven perturbers at *z* = 1.36 (filled brown circles) and the caustics at *z* = 1.36 based on model Esr2–3MP (dashed brown lines). Offsets are relative to the coordinate origin at the bottom-left corner of the HST/F160W image. The filled brown square shows the centroid of the SIS mass profile in model Esr2–3MP. The solid brown and green circles show the centroid of the SIS mass profile and the lensed source S0 projected onto the lens plane at *z* = 1.36, respectively. The position of the SIS at *z* = 1.36 is close that of the projected S0 images, showing that this perturber affects the light propagation of S0 along the LoS.</figcaption></figure>

actual LoS structure in CSWA 31 is complex and cannot be fully captured by simple mass profiles, i.e., SIS, constant mass sheet or shear. Possible asymmetries in the overall group-scale mass distribution could also play a role.

## 5. Summary

In this work we have explored the constraints on the total masses of secondary lenses obtained from multi-plane strong lensing modeling. We performed a case study of a lensed spiral galaxy at *z* = 1.49 in CSWA 31, which is along the LoS to the BGG and its group halo at *z* = 0.683. This system forms a rare, well-suited configuration with multiple image sets at different redshifts, providing tight constraints on the lens model and adequate data to infer the maximum circular velocity of the spiral.

We measured the total mass of S0 from independent strong lensing and kinematic modeling in order to assess the robustness of secondary lens mass measurements (see Sect. 4). Relying on the detailed strong lensing analysis of Wang et al. (2022), we measured the S0 mass with a range of models based on different constraints. We also modeled the velocity gradient using the least distorted image S0(d) with $\mu \approx 7$, using the state-of-the-art software 3DBarolo. We obtain a median $V_{\text{rot}} = 123^{+117}_{-38.3}$ km s$^{-1}$ at 3.3 kpc, consistent with other kinematic models of high-redshift disks of similar stellar masses (e.g., Genzel et al. 2017, 2020). The total masses of S0 derived via the two approaches are only compatible when adding a third lens plane with a concentrated mass distribution to account for LoS perturbers.

We conclude that the secondary-lens measurement from the multiple-plane modeling includes high systematic uncertainties

from degeneracies between lens planes and the complex LoS structure. Leveraging extended arcs as additional constraints can reduce the degeneracies. The secondary lens contributes less to the deflection angle of light from the source galaxies compared to the main lens, and its best-fit properties can easily be affected by other perturbers along the LoS since they project their mass onto the secondary lens plane. Accounting for additional structures along the LoS is therefore key to improving the mass measurement of the secondary lens.

*Acknowledgements.* We thank E. Di Teodoro for helpful discussions and support on 3DBarolo, and the anonymous referee for the constructive feedback that improved our paper. HW, RC and SHS thank the Max Planck Society for support through the Max Planck Research Group and Max Planck Fellowship for SHS. AG acknowledges funding and support by the Swiss National Science Foundation (SNSF). C.G. acknowledges support through grant MIUR2020 SKSTHZ.

## References

Auger, M. W., Treu, T., Bolton, A. S., et al. 2010, ApJ, 724, 511

Bacon, R., Brinchmann, J., Richard, J., et al. 2015, A&A, 575, A75

Barkana, R. 1998, ApJ, 502, 531

Bayliss, M. B., Johnson, T., Gladders, M. D., Sharon, K., & Oguri, M. 2014, ApJ, 783, 41

Belokurov, V., Evans, N. W., Hewett, P. C., et al. 2009, MNRAS, 392, 104

Birrer, S., Treu, T., Rusu, C. E., et al. 2019, MNRAS, 484, 4726

Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103

Bournaud, F., Elmegreen, B. G., & Martig, M. 2009, ApJ, 707, L1

Bradshaw, E. J., Almaini, O., Hartley, W. G., et al. 2013, MNRAS, 433, 194

Chirivi, G., Suyu, S. H., Grillo, C., et al. 2018, A&A, 614, A8

Collett, T. E. 2015, ApJ, 811, 20

Di Teodoro, E. M., & Fraternali, F. 2015, MNRAS, 451, 3021

Di Teodoro, E. M., Grillo, C., Fraternali, F., et al. 2018, MNRAS, 476, 804

Elásdóttir, Á., Limousin, M., Richard, J., et al. 2007, ArXiv e-prints [arXiv:0710.5636]

Gavazzi, R., Treu, T., Koopmans, L. V. E., et al. 2008, ApJ, 677, 1046

Genzel, R., Förster Schreiber, N. M., Übler, H., et al. 2017, Nature, 543, 397

Genzel, R., Price, S. H., Übler, H., et al. 2020, ApJ, 902, 98

Grillo, C., Suyu, S. H., Rosati, P., et al. 2015, ApJ, 800, 38

Kennicutt, R. C., Jr 1998, ARA&A, 36, 189

Koopmans, L. V. E., Treu, T., Bolton, A. S., Burles, S., & Moustakas, L. A. 2006, ApJ, 649, 599

Mandelbaum, R., Eifler, T., Hložek, R., et al. 2018, The LSST Dark Energy Science Collaboration (DESC) Science Requirements Document

McGaugh, S. S., & Schombert, J. M. 2015, ApJ, 802, 18

Meneghetti, M., Natarajan, P., Coe, D., et al. 2017, MNRAS, 472, 3177

Nestor Shachar, A., Price, S. H., Förster Schreiber, N. M., et al. 2023, ApJ, 944, 78

Newman, S. F., Genzel, R., Förster-Schreiber, N. M., et al. 2012, ApJ, 761, 43

Peng, Y.-J., Lilly, S. J., Kovač, K., et al. 2010, ApJ, 721, 193

Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526

Price, S. H., Shimizu, T. T., Genzel, R., et al. 2021, ApJ, 922, 143

Puglisi, A., Dudzevičiūtė, U., Swinbank, M., et al. 2023, MNRAS, 524, 2814

Salpeter, E. E. 1955, ApJ, 121, 161

Schneider, P. 2019, A&A, 624, A54

Shajib, A. J., Vernardos, G., Collett, T. E., et al. 2022, ArXiv e-prints [arXiv:2210.10790]

Soto, K. T., Martin, C. L., Prescott, M. K. M., & Armus, L. 2012, ApJ, 757, 86

Stark, D. P., Auger, M., Belokurov, V., et al. 2013, MNRAS, 436, 1040

Suyu, S. H., & Halkola, A. 2010, A&A, 524, A94

Suyu, S. H., Marshall, P. J., Hobson, M. P., & Blandford, R. D. 2006, MNRAS, 371, 983

Suyu, S. H., Hensel, S. W., McKean, J. P., et al. 2012, ApJ, 750, 10

Swinbank, A. M., Papadopoulos, P. P., Cox, P., et al. 2011, ApJ, 742, 11

Treu, T. 2010, ARA&A, 48, 87

Übler, H., Genzel, R., Wisnioski, E., et al. 2019, ApJ, 880, 48

Wang, H., Cañameras, R., Caminha, G. B., et al. 2022, A&A, 668, A162

Wisnioski, E., Förster Schreiber, N. M., Wuyts, S., et al. 2015, ApJ, 799, 209

Wong, K. C., Suyu, S. H., Auger, M. W., et al. 2017, MNRAS, 465, 4895

Yu, S.-Y., Xu, D., Ho, L. C., Wang, J., & Kao, W.-B. 2022, A&A, 661, A98
## Appendix A: Kinematics of S0

We present the kinematic modeling results from 3D Barolo and the stellar-mass TFR in Figs. A.1, A.2, and A.3.

<figure><img src="image_6.png" /><figcaption>**Fig. A.1.** Position-velocity diagram along the major axis of S0(d) for the [OII] doublet. The contour levels and color coding are the same as in Fig. 2. The yellow points show the position on the S0(d) image where we have secure velocity measurements. We set the $V_{\text{LOS}} = 0$ km s⁻¹ at the center of [OII]\lambda3726 Å emission line.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. A.2.** Channel maps for the kinematic models with three rings. The contour colors and levels are the same as in Fig. 2. The residual of the [OII]\lambda3729 Å emission line (bottom row) is worse than the models with four rings.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. A.3.** Stellar-mass TFR. The green points show the sample of galaxies at $z \approx 0$ from McGaugh & Schombert (2015, their Table 7), and the dashed line shows the corresponding best-fitting relation. The red point shows the S0 spiral at $z = 1.49$ from this work with $1\sigma$ uncertainty, which is consistent with the local TFR.</figcaption></figure>

## Appendix B: Three-lens-plane models

We present the degeneracies of the mass parameters in model Esr-3MP in Fig. B.1. We list the best-fit parameters in two complementary models Esr2-3MP (MS) and Esr-3MP (MS + Shear) in Table B.1 and we show their critical curves and causatics for the lens plane at $z = 1.36$ in Figs. B.2 and B.3.

<figure><img src="image_9.png" /><figcaption>**Fig. B.1.** Joint posterior PDFs for the model Esr2-3MP based on two extended sources and image positions. We show the parameters representing the lens masses at three multiple planes. The $\theta_{E,GH}$ and $\theta_{E,BGG}$ are the Einstein radius of the group halo and the BGG at the main lens plane at $z = 0.683$. The $x_{z=1.36}$ and $y_{z=1.36}$ are the centroid of the perturber, and the $\theta_{E,z=1.36}$ is its Einstein radius. The $\theta_{E,S0}$ is the Einstein radius of S0 at $z = 1.49$. The three shaded areas on the joint PDFs show the 68.3%, 95.4%, and 99.7% credible regions. The 1D histograms show the marginalized PDFs for the selected mass parameters, and the vertical lines mark the $1\sigma$ confidence intervals. The $\theta_{E,S0}$ shows strong degeneracies with $\theta_{E,BGG}$, $\theta_{E,GH}$, and $\theta_{E,z=1.36}$.</figcaption></figure>
<table><caption>**Table B.1.** Complementary three-lens-plane models with the constant mass sheet at $z = 1.36$, their constraints, the best-fit values of $\theta_{E,S0}$, and the average magnification factor, $\mu$, weighted by the surface brightness over the S0(d) image.</caption><thead><tr><th>Model name</th><th>Constraints</th><th>θ<sub>E,S0</sub></th><th>μ</th></tr></thead><tbody><tr><td>Esr2-3MP (MS)</td><td>3 sets + extended arcs<br>of S0, S3</td><td>1.51"</td><td>6.6</td></tr><tr><td>Esr2-3MP (MS+Shear)</td><td>3 sets + extended arcs<br>of S0, S3</td><td>2.10"</td><td>6.5</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Fig. B.2.** Critical curves for the lens plane at $z = 1.36$ in three-lens-plane models. The Esr2-3MP model uses an isothermal mass profile to model the mass distribution at $z = 1.36$. Meanwhile, the Esr2-3MP (MS) model uses a constant mass sheet, and the Esr2-3MP (MS + Shear) model uses a constant mass sheet with an external shear. The purple and gray arrows indicate the mass concentrations predicted by the external shear at $z = 1.36$ and $z = 0.683$, respectively, in the model Esr2-3MP (MS + Shear). The critical curves for these models nearly overlap, with the orange and blue lines largely hidden beneath the purple line.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. B.3.** Caustics for lens plane at $z = 1.36$ in three-lens-plane models.</figcaption></figure>