© The Authors 2024

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 ≥ 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 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₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7. Hence, 1″ corresponds to 7.08 kpc at the main lens redshift of z = 0.683, and to 8.46 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 ±3000 km s⁻¹ of the BGG redshift were selected from the MUSE spectroscopic catalog.

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 3DBarolo.

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_eff of 1.17″, or about 3.8 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 ≈ 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−MP_test 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 μ ≃ 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″ × 0.73″. Its wavelength resolution is FWHM = 2.35 Å, with a channel width of 1.25 Å. To focus on the [OII] doublet, the cube was cropped between 3719 Å to 3734 Å 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 × 10⁻¹⁸ erg s⁻¹. 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 ^3DBarolo.

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 ≥ 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 ^3DBarolo (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 ^3DBarolo, dominates over systematic errors due to lensing shear.

^3DBarolo 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, ^3DBarolo replicates the [OII] emission doublet at λ1 = 3726 Å and 3729 Å. Each ring is defined with a centroid, position angle, inclination angle i, systematic velocity V_sys, rotation velocity V_rot and velocity dispersion σ_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_sys in all tilted rings to 0 km s⁻¹. We iteratively refined these values with ^3DBarolo and subsequently allowed variations in only the V_rot and σ_v at each ring (see Fig. A.1). In the fitting, we adopted the local normalization and a uniform weighting function¹ to account for the asymmetry in the intensity map of S0(d).

^3DBarolo 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 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°.

Fig. 2. Comparison between MUSE IFS of image S0(d) (blue contours) and the 3DBarolo 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−20 erg s−1 cm−2 Å−1 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.

Due to the moderate spectral resolution of the MUSE data cube, velocity dispersion values returned by ^3DBarolo 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 V_rot 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).

Fig. 3. Histogram of Vrot weighted by the BIC. The upper and bottom panels show the weighted Vrot distribution at 2.0 kpc and at 3.3 kpc from the galaxy centroid in the source plane, respectively. The 1σ uncertainties on Vrot are determined by drawing 5000 samples of the inclination angles in the range [5°, 35°].

We determined the total mass of S0 using V_rot via

$$\begin{array}{c}\hfill M_{\text{rot}}(r\prime )=\frac{V_{\text{rot}}^2}{G}r\prime ,\end{array}$$(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′² 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

$$\begin{array}{c}\hfill r\prime =\frac{r}{\sqrt{\mu }}.\end{array}$$(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).

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σ 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σ 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 μ from the baseline model. The error bar shows 1σ uncertainties associated with the range of inclination angles.

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, S0 deflects the light from the background sources S3, S4, and S5 at higher redshifts. Accounting for this secondary strong deflector (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

$$\begin{array}{c}\hfill \rho (r)=\frac{{{\sigma }_{\text{v}}}^2}{2\pi G{r}^2},\end{array}$$(3)

where σ_v is the velocity dispersion, related to the Einstein radius θ_E (for a lensed background source at infinite redshift) via ${{\sigma }_{\text{v}}}^2=\frac{{\theta }_{\text{E}}{c}^2}{4\pi }$. Thus, Eq. (3) can be written in terms of θ_E:

$$\begin{array}{c}\hfill \rho (r)=\frac{{\theta }_{\text{E}}{c}^2}{8{\pi }^{2}G{r}^2}.\end{array}$$(4)

We calculated the mass of S0 enclosed within radius r′ in spherical coordinates,

$$\begin{array}{c}\hfill M_{\text{lens}}(r\prime )={\int }_0^{r\prime }4\pi r^2\frac{{\theta }_{\text{E}}{c}^2}{8{\pi }^{2}G{r}^2}\mathrm{d}r,\end{array}$$(5)

such that

$$\begin{array}{c}\hfill M_{\text{lens}}(r\prime )=\frac{{\theta }_{\text{E}}{c}^2}{2\pi G}r\prime ,\end{array}$$(6)

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 constraints to ascertain the best-fit values of the Einstein radius θ_E, S0 (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 predicted 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 significantly more mass to S0 than the modeling of gas kinematics. We discuss this discrepancy by highlighting the main uncertainties 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 σ_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,

$$\begin{array}{c}\hfill w_i=\frac{exp\{-({\mathrm{BIC}}_i-{\mathrm{BIC}}_{\min })/2\}}{{\sum }_{i=1}^{5000}exp\{-({\mathrm{BIC}}_i-{\mathrm{BIC}}_{\min })/2\}}.\end{array}$$(7)

We then determined the median V_rot and the 1σ 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 km s⁻¹ and σ_v = 40 km s⁻¹ at 3.3 kpc from the galaxy center, and a ratio V_rot/σ_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 ≈ 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 σ_v ≃ 27 km s⁻¹ indicate stronger rotation support, the properties of both spirals at z ≈ 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 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 ± 0.4)×10¹⁰ M_⊙, corrected by μ = 6.7 from the baseline model, showing that S0 follows the stellar-mass TFR at z ≃ 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 ^3DBarolo 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σ 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 ≳ 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 ≳ 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 ≈ 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 θ_E, S0 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, θ_E, S0 decreases from ≃4″ to 1.65″.

To further probe this difference, we checked the impact of changing the prior on θ_E, S0, 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

$$\begin{array}{c}\hfill {\theta }_{\text{E}}=\frac{2\pi }{c^2}{V}_{\text{rot}}^2.\end{array}$$(8)

In this case, the best-fit θ_E, S0 in the baseline model corresponds to a constant V_rot of 338.3 km s⁻¹ as a function of radius, a value above the 1σ upper limit from ^3DBarolo (see Fig. 3). On the other hand, the median rotation velocities correspond to θ_E ≃ 0.2″. We assigned a Gaussian prior 𝒢 (0.2″,  0.07″) to θ_E, S0 and recomputed the baseline model. The best-fit θ_E, S0 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 θ_E, S0 = 0.85″ and θ_E, z = 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 θ_E, S0 falls back to that of the baseline model. In the Esr2−3MP model, we assigned broad flat priors in the range [0″, 5″] for θ_E, S0 and θ_{E, z = 1.36}, and obtain strong degeneracies between the mass profiles in the three lens planes (see Fig. B.1). The θ_E, S0 is negatively correlated with the Einstein radius of the group-halo θ_E, GH and with θ_{E, z = 1.36}, and is positively correlated with the Einstein radius of the BGG θ_E, BGG. Our model Esr2−3MP provides comparable fitting results as baseline. It improves the modeling of extended arcs for S0 and S3 by 𝒪(100) in the loglikelihood², which indicates the goodness-of-fit (Suyu et al. 2006) but in the meantime, the χ² of image position modeling worsens for S2, S4, and S5 by comparable amounts.

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.

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 𝒪(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 θ_E, S0 = 2.01″ (see Table B.1); this translates to V_rot = 382 km s⁻¹, which is higher than the V_rot 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 θ_{E, z = 1.36} = 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 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 μ ≃ 7, using the state-of-the-art software ^3DBarolo. We obtain a median $V_{\mathrm{rot}}={123}_{-38.3}^{+117}$ng 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.

Acknowledgments

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

Appendix A: Kinematics of S0

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

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 VLoS = 0 km s−1 at the center of [OII]λ3726 Å emission line.

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]λ3729 Å emission line (bottom row) is worse than the models with four rings.

Fig. A.3. Stellar-mass TFR. The green points show the sample of galaxies at z ≃ 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σ uncertainty, which is consistent with the local TFR.

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 caustics for the lens plane at z = 1.36 in Figs. B.2 and B.3.

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 θE, GH and θE, BGG are the Einstein radius of the group halo and the BGG at the main lens plane at z = 0.683. The xz = 1.36 and yz = 1.36 are the centroid of the perturber, and the θE, z = 1.36 is its Einstein radius. The θ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σ confidence intervals. The θE, S0 shows strong degeneracies with θE, BGG, θE, GH, and θE, z = 1.36.

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.

Fig. B.3. Caustics for lens plane at z = 1.36 in three-lens-plane models.

All Tables

All Figures

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 3DBarolo.

In the text

Fig. 2. Comparison between MUSE IFS of image S0(d) (blue contours) and the 3DBarolo 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−20 erg s−1 cm−2 Å−1 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.

In the text

	Fig. 3. Histogram of Vrot weighted by the BIC. The upper and bottom panels show the weighted Vrot distribution at 2.0 kpc and at 3.3 kpc from the galaxy centroid in the source plane, respectively. The 1σ uncertainties on Vrot are determined by drawing 5000 samples of the inclination angles in the range [5°, 35°].
In the text

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σ 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σ 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 μ from the baseline model. The error bar shows 1σ uncertainties associated with the range of inclination angles.

In the text

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.

In the text

	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 VLoS = 0 km s−1 at the center of [OII]λ3726 Å emission line.
In the text

	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]λ3729 Å emission line (bottom row) is worse than the models with four rings.
In the text

	Fig. A.3. Stellar-mass TFR. The green points show the sample of galaxies at z ≃ 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σ uncertainty, which is consistent with the local TFR.
In the text

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 θE, GH and θE, BGG are the Einstein radius of the group halo and the BGG at the main lens plane at z = 0.683. The xz = 1.36 and yz = 1.36 are the centroid of the perturber, and the θE, z = 1.36 is its Einstein radius. The θ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σ confidence intervals. The θE, S0 shows strong degeneracies with θE, BGG, θE, GH, and θE, z = 1.36.

In the text

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.

In the text

	Fig. B.3. Caustics for lens plane at z = 1.36 in three-lens-plane models.
In the text