# HerS-3: An Exceptional Einstein Cross Reveals a Massive Dark Matter Halo

P. Cox, K. Butler, C. Keeton, L. Eid, E. Borsato, T. Bakx, R. Neri, B. Jones, P. Prajapati, A. Baker, et al.

► To cite this version:

P. Cox, K. Butler, C. Keeton, L. Eid, E. Borsato, et al. HerS-3: An Exceptional Einstein Cross Reveals a Massive Dark Matter Halo. The Astrophysical Journal, 2025, 991 (1), pp.53. 10.3847/1538-4357/adf204. insu-05281537

HAL Id: insu-05281537

https://insu.hal.science/insu-05281537v1

Submitted on 24 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# HerS-3: An Exceptional Einstein Cross Reveals a Massive Dark Matter Halo

P. Cox<fnref n="1" />, K. M. Butler<fnref n="2" />, C. R. Keeton<fnref n="3" />, L. Eid<fnref n="3" />, E. Borsato<fnref n="4" />, T. J. L. C. Bakx<fnref n="2" />, R. Neri<fnref n="5" />, B. M. Jones<fnref n="6" />, P. Prajapati<fnref n="7" />, A. J. Baker<fnref n="3" />,<fnref n="8" />, S. Berta<fnref n="5" />, A. Cooray<fnref n="9" />, E. M. Corsini<fnref n="4" />,<fnref n="10" />, L. Marchetti<fnref n="11" />,<fnref n="12" />, A. Omont<fnref n="1" />, A. Beelen<fnref n="13" />, R. Gavazzi<fnref n="13" />, D. Ismail<fnref n="14" />, R. J. Ivison<fnref n="15" />,<fnref n="16" />, M. Krips<fnref n="5" />, M. D. Lehnert<fnref n="1" />,<fnref n="17" />, H. Messias<fnref n="18" />,<fnref n="19" />, D. Riechers<fnref n="6" />, C. Vlahakis<fnref n="20" />, A. Weiss<fnref n="7" />, P. van der Werf<fnref n="21" />, and C. Yang<fnref n="2" />

Received 2025 March 5; revised 2025 July 16; accepted 2025 July 18; published 2025 September 16

## Abstract

We present a study of HerS-3, a dusty star-forming galaxy at $z_{spec} = 3.0607$, which is gravitationally amplified into an Einstein cross with a fifth image of the background galaxy seen at the center of the cross. Detailed 1 mm spectroscopy and imaging with NOEMA and the Atacama Large Millimeter/submillimeter Array resolve the individual images and show that each of the five images display a series of molecular lines that have similar central velocities, unambiguously confirming that they have identical redshifts. The Hubble Space Telescope F110W image reveals a foreground lensing group of four galaxies with a photometric redshift $z_{phot} \sim 1.0$. Lens models that only include the four visible galaxies are unable to reproduce the properties of HerS-3. By adding a fifth massive component, lying southeast of the brightest galaxy of the group, the source reconstruction is able to match the peak emission, shape, and orientation for each of the five images. The fact that no galaxy is detected near that position indicates the presence of a massive dark matter halo in the lensing galaxy group. In the source plane, HerS-3 appears as an infrared luminous starburst galaxy seen nearly edge on. The serendipitous discovery of this exceptional Einstein cross offers a potential laboratory for exploring at small spatial scales a nuclear starburst at the peak of cosmic evolution and studying the properties of a massive dark matter halo associated with the lensing galaxy group.

*Unified Astronomy Thesaurus concepts: Strong gravitational lensing (1643)*

## 1. Introduction

Strong gravitational lensing occurs when light from a distant galaxy passes by a massive galaxy (or a group of galaxies) that distorts spacetime and causes the path of light of the background source to bend, making it appear brighter and producing multiple images, arcs, or a ring. The resulting amplification enables studies of galaxies in the early Universe at spatial scales below 100 pc with current facilities. Such systems are powerful astrophysical laboratories for constraining the properties of dark matter in galaxies, groups, and clusters (e.g., Y. D. Hezaveh et al. 2016; D. Gilman et al. 2019; A. C. Sengül & C. Dvorkin 2022;

K. T. Inoue et al. 2023; P. Natarajan et al. 2024; A. M. Nierenberg et al. 2024; W. Sheu et al. 2024). Multiply imaged sources with time variability related to quasars or supernova explosions can also be used to monitor time delays and measure the Hubble constant (e.g., S. Refsdal 1964; C. S. Kochanek & P. L. Schechter 2004; G. C. F. Chen et al. 2019; T. Schmidt et al. 2023).

For a nearly perfect alignment and an elliptical lens mass distribution, the background source will appear as a quadruply imaged system. The archetypel of such systems is the Einstein cross, where four distinct images of the background source form a cross-like pattern with a high degree of symmetry. It remains a challenge finding and confirming such rare sources, and, as of today, the majority of cases have been discovered using optical or near-infrared imaging.

The first calculations of gravitational lensing for stars (as both lens and source) were published by O. Chwolson (1924)

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Sorbonne Université, UPMC Paris 6 and CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98b bd. Arago, 75014 Paris, France; pierre.cox@iap.fr</footnote>

<footnote n="2">² Department of Space, Earth & Environement, Chalmers University of Technology, Chalmersplatsen 4, Gothenburg SE-412 96, Sweden</footnote>

<footnote n="3">³ Department of Physics and Astronomy, Rutgers University—New Brunswick, 136 Frelinghuysen Rd., Piscataway, NJ 08854-8019, USA</footnote>

<footnote n="4">⁴ Dipartimento di Fisica e Astronomia "G. Galilei," Università di Padova, vicolo dell'Osservatorio 3, I-35122 Padova, Italy</footnote>

<footnote n="5">⁵ Institut de Radioastronomie Millimétrique (IRAM), 300 rue de la Piscine, 38406 Saint-Martin-d'Hères, France</footnote>

<footnote n="6">⁶ Institut für Astrophysik, Universität zu Köln, Zülpicher Strasse 77, D-50937 Köln, Germany</footnote>

<footnote n="7">⁷ Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany</footnote>

<footnote n="8">⁸ Department of Physics and Astronomy, University of the Western Cape, Robert Sobukwe Rd., Bellville 7535, Cape Town, South Africa</footnote>

<footnote n="9">⁹ University of California Irvine, Department of Physics & Astronomy, FRH 2174, Irvine, CA 92697, USA</footnote>

<footnote n="10">¹⁰ INAF—Osservatorio Astronomico di Padova, vicolo dell'Osservatorio 5, I-35122 Padova, Italy</footnote>

<footnote n="11">¹¹ IDIA & Departement of Astronomy, University of Cape Town, 7701 Rondebosch, Cape Town, South Africa</footnote>

<footnote n="12">¹² Istituto Nazionale di Astrofisica, Istituto di Radioastronomia—Italian ARC, Via Piero Gobetti 101, 40129 Bologna, Italy</footnote>

<footnote n="13">¹³ Aix-Marseille Université, CNRS and CNES, LAM, 38 rue Frédéric Joliot-Curie, 13388 Marseille, France</footnote>

<footnote n="14">¹⁴ Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, 67000 Strasbourg, France</footnote>

<footnote n="15">¹⁵ Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="16">¹⁶ School of Cosmic Physics, Dublin Institute for Advanced Studies, 31 Fitzwilliam Pl., Dublin D02 XF86, Ireland</footnote>

<footnote n="17">¹⁷ Université Lyon 1, ENS de Lyon, Centre de Recherche Astrophysique de Lyon (UMR5574), 69230 Saint-Genis-Laval, France</footnote>

<footnote n="18">¹⁸ European Southern Observatory, Alonso de Córdova 3107, Casilla 19001, Vitacura, Santiago, Chile</footnote>

<footnote n="19">¹⁹ Joint ALMA Observatory, Alonso de Córdova, 3107, Vitacura, Santiago 763-0355, Chile</footnote>

<footnote n="20">²⁰ National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22901, USA</footnote>

<footnote n="21">²¹ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>
and, later, by A. Einstein (1936), who noted that “there is no great chance of observing this phenomenon.” This led F. Zwicky (1937) to point out that extragalactic nebulae (i.e., galaxies) “offer a much better chance than stars for the observation of gravitational lens effects.” Forty years later, the first known lens was discovered in the optical by D. Walsh et al. (1979), who showed that the quasar QSO-0956+561, at $z = 1.405$, was lensed in two images separated by $6''$. A second discovery of a gravitationally amplified source was published soon thereafter, the quasar PG 1115+080, a quadruple lensed system at $z = 1.722$ (R. J. Weymann et al. 1980; E. K. Hege et al. 1981). The third detection was made a few years later by J. Huchra et al. (1985) who reported a quasar lensed in a cross pattern, QSO-2237+0305 at $z = 1.695$, known as the Einstein Cross (see also H. K. C. Yee 1988; G. Adam et al. 1989). Many other examples of quasars lensed by foreground galaxies were found since then, including further cases of Einstein crosses. J. R. Lucey et al. (2018) reported two diamond shaped systems of quadruply imaged quasars; L. Wisotzki et al. (2002) discovered a lensed $z = 1.69$ quasar (HE 0435−1223) with four images in a cross-shaped arrangement around a massive elliptical galaxy at $z = 0.45$ (see also K. C. Wong et al. 2017); one Einstein cross was identified in a sample of six quadruple-image lensed quasars by K. C. Wong et al. (2020); and D. Stern et al. (2021) presented a dozen quadruply imaged systems discovered in Gaia Data Release 2, including three Einstein crosses.

Einstein crosses produced by galaxy–galaxy strong lensing are equally rare with only a dozen confirmed cases among a few hundred known lensed galaxies. K. U. Ratnatunga et al. (1995) reported the serendipitous discovery of two cases of such lenses using the Hubble Space Telescope (HST). A. S. Bolton et al. (2006) identified a cruciform configuration of an $Ly\alpha$ emitter galaxy at $z = 2.701$ lensed by a luminous red galaxy at $z \sim 0.33$. M. Pettini et al. (2010) presented the case of CSWA 20, a blue star-forming galaxy at $z = 1.433$ lensed by a luminous red galaxy at $z = 0.741$ into a four-imaged configuration reminiscent of an Einstein cross, with a large separation ($\sim 6''$) between the images. D. Bettoni et al. (2019) demonstrated the lens nature of a cross-like system by deep spectroscopic observations, showing that it consisted of a Lyman-break galaxy at $z = 3.03$ lensed by a galaxy at $z = 0.556$. N. R. Napolitano et al. (2020) discovered two Einstein crosses in the footprint of the Kilo-Degree Survey, and A. Cikota et al. (2023) reported on a strong lensing system, DESI 253.2534+26.8843, a starburst galaxy at $z = 2.597$ lensed into an Einstein cross by an elliptical galaxy at $z = 0.636$.

Strongly lensed systems have also been discovered using submillimeter facilities, mostly through follow-up observations of infrared luminous dust-enshrouded starburst galaxies in the early Universe selected from the Herschel, Planck, and South Pole Telescope surveys (J. D. Vieira et al. 2010; S. Eales et al. 2010; S. J. Oliver et al. 2012; Planck Collaboration et al. 2015). The majority of these gravitational lensed systems, now in excess of 300 cases, appear as complete or partial Einstein rings or extended arcs, with no case known to date of an Einstein cross. One famous example of a lensed submillimeter galaxy is SDP.81 at $z = 3.042$, magnified by a galaxy at $z = 0.3$ (M. Negrello et al. 2010; A. Omont et al. 2013), which has been studied in exquisite sensitivity and angular resolution (23 mas) with the Atacama Large Millimeter/submillimeter Array

(ALMA) to reveal the full details of a nearly complete Einstein ring (ALMA Partnership et al. 2015). These high-quality ALMA data have also provided evidence of dark matter substructure in the foreground lensing galaxy (Y. D. Hezaveh et al. 2016). Another example is the Cloverleaf, a quadruply imaged starburst at $z = 2.55$, harboring a broad absorption line quasar, lensed by a $z \sim 1.7$ galaxy (R. Barvainis et al. 1994; J. P. Kneib et al. 1998; G. Chartas et al. 2004), in which HST and ALMA images revealed a partial Einstein ring linking the four image components (V. Chantry & P. Magain 2007; L. Zhang et al. 2023). Other quadruple lensed sources, which were found in the submillimeter, all show lensing arcs connecting the images, e.g., HLSW-01 at $z = 2.95$ (A. Conley et al. 2011), G09v.1.97 at $z = 3.63$ (C. Yang et al. 2019), and J0314−4452 at $z = 2.93$ and J0150−5924 at $z = 2.78$ (D. A. Riechers 2025).

Detailed follow-up observations of high-$z$ galaxies detected in far-infrared and submillimeter continuum surveys are only possible when robust spectroscopic redshifts are available. Using the IRAM Northern Extended Millimeter Array (NOEMA), the $z$-GAL spectroscopic redshift survey of bright Herschel-selected dusty star-forming galaxies (DSFGs) has recently delivered precise redshifts for 165 individual galaxies, representing the largest flux-limited sample of DSFGs with unambiguous redshifts obtained to date (R. Neri et al. 2020; P. Cox et al. 2023). The $z$-GAL survey provides a foundation for exploring in detail the properties of these high-$z$ galaxies.

In a recent follow-up study to trace fueling and feedback in $z$-GAL selected galaxies (K. M. Butler et al. 2025, in preparation), one gravitationally amplified galaxy stood out for its unusual lensed morphology. The source, HerS-3 at $z = 3.0607$, was revealed to be an exceptional Einstein cross that, uniquely compared to all previous examples, shows a central fifth image.

In this paper, we present a multiwavelength study of HerS-3 and its surroundings, including submillimeter, radio, and optical/near-infrared observations, focusing on the origin of the remarkable lensed image of this galaxy. In Section 2, we describe the series of observations that are used in this study. In Section 3, we report the main results from the multi-wavelength observations for the molecular emission lines as well as the optical/near-infrared, submillimeter, and radio continuum emission. Section 4 presents the gravitational lensing model and the results of the analysis, including the implied mass distribution of the foreground lensing system. Section 5 outlines the implications of the results of this study, addressing the main properties of HerS-3 in the source plane informed by the lensing model (Section 5.1), discussing the nature of the additional mass component needed to reproduce the characteristics of the five images of the Einstein cross, and presenting them in a broader context by highlighting the potential use of this cosmic laboratory (Sections 5.2 and 5.3). Finally, Section 6 summarizes the main findings of this paper and outlines future observational prospects.

Throughout this paper, we adopt a spatially flat Lambda cold dark matter (ΛCDM) cosmology with $H_0 = 67.4 \text{ km s}^{-1} \text{ Mpc}^{-1}$ and $\Omega_M = 0.315$ (Planck Collaboration et al. 2020). With these parameters, $1''0$ corresponds to $\sim 7.8 \text{ kpc}$ at $z = 3.06$ in the source plane and the luminosity distance to the source is $D_L = 2.7 \times 10^4 \text{ Mpc}$; in the lens plane, at $z \sim 1.0$ (see Section 3.1.2), $1''0$ corresponds to $\sim 8.2 \text{ kpc}$.
## 2. Observations

### 2.1. NOEMA

We used NOEMA to target high-frequency molecular lines in HerS-3, redshifted into the 1 mm band. The observations were carried out under project W21DF (PI: P. Cox). The project was completed on 2022 January 15 with 12 antennas using the intermediate C-configuration. The source HerS-3 was observed for a total observing time of 1.13 hr, which resulted in a coherent set of high-sensitivity observations in Band 3 covering the two strongest ground-state $OH^+$ lines, corresponding to the redshifted transitions of $OH^+(1_1-0_1)$ and $OH^+(1_2-0_1)$ ($\nu_{\text{rest}} = 1033.118$ GHz and 971.803 GHz, respectively), together with the $^{12}$CO(9–8) and $H_2O(2_{02}-1_{11})$ emission lines ($\nu_{\text{rest}} = 1036.912$ GHz and 987.927 GHz, respectively), and the underlying dust continuum emission.

Observing conditions were excellent with an atmospheric phase stability rms of 20° on the longest baseline (400 m) and 1.5 mm of precipitable water vapor. The correlator was operated in the low-resolution mode to provide spectral channels with a nominal resolution of 2 MHz. The NOEMA antennas were equipped with 2SB receivers that cover a spectral window of 7.744 GHz in each sideband and polarization. We covered the frequency range from 237.4 to 245.2 GHz in the lower sideband (LSB) and 252.8–260.5 GHz in the upper sideband (USB).

The flux calibrators used were MWC349 and LkHα101, and the phase calibrator was 0106+013. The data were calibrated, averaged in polarization, mapped, and analyzed in the GILDAS software package. The absolute flux calibration is accurate to within 10%. The sensitivities reached in the 40 MHz ($\sim$50 km s⁻¹) channels are 0.90 mJy beam⁻¹ for the USB and 0.88 mJy beam⁻¹ for the LSB.

A continuum uv table of HerS-3 was produced by averaging the visibilities from the line-free sections of the USB and LSB. The USB continuum visibilities and their weights were scaled based on an assumed spectral index of 2.28 and aligned in the uv plane to the reference frequency of 239.301 GHz. The continuum flux density of HerS-3 was sufficiently strong to allow phase self-calibration down to the thermal noise level of the visibilities. A continuum dirty map was constructed using calibrated visibilities, applying a Briggs robust weighting parameter of 0.1. The map was cleaned down to the thermal noise level with the CLARK algorithm, using a mask to exclude regions without continuum emission. This process achieved high image fidelity and an imaging dynamic range, defined as the ratio of peak flux density to rms noise in the map, which exceeded 30 in the dust continuum emission. The resulting map has an angular resolution of 0.92 × 0.52 (a position angle (PA) of 23°) and a 1σ sensitivity of 69 µJy beam⁻¹.

For our spectral analysis of the data, we created separate uv tables of the USB and LSB. We imaged both sidebands using the Hogbom algorithm and a natural weighting scheme to maximize the signal-to-noise ratio (SNR). A mask was applied to exclude regions outside the continuum and line emission regions before cleaning down to the thermal noise of each sideband. The resulting data cubes are binned in frequency by a factor of 30, resulting in channel widths of 60 MHz or 70.7 km s⁻¹ at the redshifted rest frequency of the $OH^+(1_1-0_1)$ line.

### 2.2. ALMA

The source HerS-3 was part of the ALMA Cycle 9 project ADS/JAO.ALMA 2022.1.00145.S (PI: T. Bakx) to obtain short Band 7 continuum observations of Herschel-selected high-z galaxies with robust spectroscopic redshifts. The source was observed at 291.8 GHz for a total on-source time of $\sim$3.0 minutes with 4 × 128 dual polarization channels covering a bandwidth of 7.875 GHz with a 4 GHz sideband separation. The observations were done on 2023 June 1. The bandpass and phase calibrators used were J0006–0623 and J0125–0005, respectively. Observing conditions were excellent with a median precipitable water vapor at zenith of 0.67 ± 0.03 mm, a phase stability of 9.°1 rms, and a system temperature $T_{\text{sys}} = 91$ K in configuration C43-6 with a baseline length between 27 m and 3637 m. The observations were carried out with 45 antennas and a maximum baseline of 1 km, yielding an angular resolution of 0.′12 × 0.′10 using a robust = 2.0 weighting and a continuum sensitivity of 63 µJy beam⁻¹ (1σ). The ALMA data were reduced and imaged using the Common Astronomy Software Application<fnref n="22" /> (CASA; v6.5.4.9; J. P. McMullin et al. 2007; CASA Team et al. 2022).

The continuum flux densities of each of the five images were extracted from the primary-beam-corrected fields, using apertures that fit all of the continuum emission and extend two beams beyond the 2σ contours, in an effort to recover all the flux. The errors on these extracted flux densities were then calculated from the per-beam error, multiplied by the square root of the number of beams within the aperture. Finally, an additional 10% error was added in quadrature to account for calibration errors, as recommended by the ALMA Technical Handbook (P. Cortés et al. 2024). Although the images are extended, the peak flux positions are used as the central positions listed in Table 1.

### 2.3. VLA

We observed the radio continuum emission of HerS-3 using the NSF's Karl G. Jansky Very Large Array (VLA) as part of a 6 GHz continuum survey of all the z-GAL sources (program IDs VLA/22A-211 and VLA/23A-030; PI: T. Bakx). A complete description of the results of this survey will be presented in T. Bakx et al. (2025, in preparation). The observations were done in the A-configuration in order to achieve an angular resolution of $\sim$0.′4. We used a 4 GHz bandwidth and 3 bit samplers to maximize bandwidth and continuum sensitivity, using the standard 1 MHz ($\sim$50 km s⁻¹) frequency resolution (dual polarization) in broadband mode. The data were acquired in 2023 November under stable atmospheric conditions and HerS-3 was observed for a total time of 40 minutes (27 minutes on-source integration). 3C48 was used for bandpass and flux calibration, and J0115-0127 as the gain calibrator.

The data were reduced, calibrated, and imaged using CASA v6.5.4.9 (J. P. McMullin et al. 2007; CASA Team et al. 2022). The accuracy of the flux density calibration is within 10%–15%, when comparing the measured fluxes for the calibrators from our observations to the calibrator models. The continuum map was created and cleaned using the TCLEAN function of CASA with natural weighting.

<footnote n="22">²² http://casa.nrao.edu</footnote>
<table><caption>**Table 1**

Coordinates and Continuum Flux Densities of the Individual Images of the HerS-3 Einstein Cross</caption><thead><tr><th rowspan="3">Image</th><th colspan="2">Coordinates (J2000)</th><th colspan="2">NOEMA</th><th rowspan="2">ALMA<br>292 GHz</th><th rowspan="2">VLA<br>6 GHz</th></tr><tr><th>R.A.</th><th>Decl.</th><th>241 GHz</th><th>258 GHz</th></tr><tr><th colspan="2">(HH:MM:SS.SS)</th><th></th><th>(mJy)</th><th></th><th></th></tr></thead><tbody><tr><td>SW</td><td>01:27:54.01</td><td>00:49:37.1</td><td>4.05 ± 0.16</td><td>5.81 ± 0.18</td><td>5.70 ± 0.39</td><td>0.075 ± 0.015</td></tr><tr><td>NE</td><td>01:27:54.37</td><td>00:49:42.1</td><td>4.27 ± 0.17</td><td>6.23 ± 0.19</td><td>8.48 ± 0.59</td><td>0.046 ± 0.013</td></tr><tr><td>E</td><td>01:27:54.26</td><td>00:49:38.5</td><td>2.40 ± 0.08</td><td>2.37 ± 0.09</td><td>2.17 ± 0.31</td><td>0.018 ± 0.007</td></tr><tr><td>C</td><td>01:27:54.15</td><td>00:49:40.0</td><td>1.82 ± 0.08</td><td>2.11 ± 0.09</td><td>3.07 ± 0.33</td><td>0.048 ± 0.014</td></tr><tr><td>W</td><td>01:27:54.06</td><td>00:49:41.1</td><td>1.40 ± 0.08</td><td>0.80 ± 0.06</td><td>1.82 ± 0.33</td><td>&lt;0.0012</td></tr></tbody></table>

Note. The coordinates correspond to the center of each image of the HerS-3 Einstein cross, as derived from the ALMA 292 GHz dust continuum. The flux densities listed are the measured values with the 1σ uncertainty. In case of a nondetection, the flux density is given as a 3σ upper limit.

## 2.4. HST Imaging

The HST imaging of HerS-3 is part of a snapshot observations project (ID 15242; PI: L. Marchetti), which was carried out from 2017 October to 2018 June. The observations were performed using the wide YJ filter F110W of the Wide Field Camera 3 (WFC3) infrared channel. HerS-3 was observed with a total exposure time of 711.74 s. The HST data were obtained from the Mikulski Archive for Space Telescopes (MAST) and the specific observations here analyzed can be accessed via DOI: 10.17909/vg9e-2714. However, these HST data were rerReduced, as explained in E. Borsato et al. (2024), adopting, in particular, a different drizzling strategy than the standard one used for the MAST data products by setting the output pixel scale and linear drop size equal to 0.0"064 pixel⁻¹ and 0.8, respectively. Details of the observations, data reduction, and computation of the background, noise map, source mask, and point-spread function (PSF) are provided in E. Borsato et al. (2024). The final HST pixel scale is 0.0"064. Importantly, astrometric corrections were performed by matching the Gaia early Data Release 3 sources (Gaia Collaboration et al. 2021) with their nearest HST counterparts. We found three Gaia sources for which we computed the mean absolute offsets $⟨d_N = 0.0"36⟩$ and $⟨d_E = 0.0"17⟩$ along the north and east directions. Computing the rms of the residual offsets after correcting for the mean shifts, we obtain a value of rms = 0.0"022. Since errors of both the Gaia Data Release 3 (DR3) astrometry and the relative HST sources are much smaller than the residual rms, we adopted the latter as the error on the absolute astrometry of the HST data.

## 2.5. Subaru and SDSS Imaging

In order to recover the optical rest-frame emission of the lensing galaxies, we used stacked Sloan Digital Sky Survey (SDSS) and Subaru/Hyper Supreme-Cam (HSC) observations. We downloaded Subaru/HSC data from the Subaru HSC public DR3 page. We retrieved ~10 × 10 arcmin² cutouts around the centroid of the SPIRE detection of HerS-3 in the g, r, i, z, and Y bands. We also retrieved similar cutouts of the noise maps and, through the HSC PSF picker (J. Bosch et al. 2018), cutouts of PSF models near the position of the HerS-3 system. The Subaru/HSC pixel scale is 0.0"168. To add some constraints in the blue, we used stacked SDSS data. We downloaded the deep stacks computed by L. Jiang et al. (2014)<fnref n="23" /> in the u band. Each stacked image is a ~14 × 18

arcmin² band covering the SDSS Stripe 82. The details on the stacking procedure are available in L. Jiang et al. (2014). To construct a noise map and the PSF for the SDSS data, we followed the same procedure applied in E. Borsato et al. (2024). The SDSS pixel scale is 0.0"396.

## 3. Results

The combination of high-quality multiwavelength observations in the submillimeter, radio, and near-infrared allow us to study the configurations and properties of both the lensed background galaxy, HerS-3, and the foreground lensing galaxy group. In this section, we will describe the field around HerS-3 and the properties of the foreground galaxies (Section 3.1), then present the detailed characteristics of the Einstein cross (Section 3.2), analyze the spectral energy distribution (SED) of HerS-3 (Section 3.3), and outline the properties of the molecular lines in Section 3.4.

### 3.1. Overall View of the HerS-3 Field

#### 3.1.1. The Einstein Cross

In the z-GAL survey, HerS-3 was observed at 3 and 2 mm at moderate angular resolution (~5"). The source was found to be extended over ~10" in the northeast-southwest direction, with two distinct peaks detected in the $^{12}$CO(3–2) and (5–4) emission lines, yielding a redshift of $z_{spec} = 3.0607$ for both components (P. Cox et al. 2023). The new NOEMA observations, done at 1 mm with a beam of 0.0"92 × 0.0"52, reveal that HerS-3 is an Einstein cross that includes a bright central fifth image (Figure 1; the coordinates of each of the continuum images are listed in Table 1).

The NE and SW images, which are the brightest, are separated by 7"5. This is an unusually wide separation as compared to other Einstein crosses or quadruple imaged quasars, where the distance between the images is typically of a few arcseconds—notable exceptions are the lensed quadrupole quasars J1004+412 and J165105–041725 that have maximum separations between the components of 14"6 and 10", respectively (N. Inada et al. 2003; D. Stern et al. 2021). The three images along the southeast-northwest direction have a smaller separation and extend over ~4". It is also noteworthy that the alignment of the NE and SW images with the image in the middle is slightly bent.

All five images, which have various molecular emission lines detected, are confirmed to be at the same redshift ($z_{spec} = 3.0607$) and relative velocities, demonstrating that they originate from the same physical source (Figure 2 and

<footnote n="23">²³ Available in http://das.sdss.org/ge/sample/stripe82/.</footnote>
<figure><img src="image_4.png" /></figure>

**Figure 1.** The 1.2 mm dust continuum emission of HerS-3 observed with NOEMA and obtained with a robust weighting of 0.1 (shown in yellow contours) reveals that the background starburst galaxy at $z = 3.0607$ is lensed into an Einstein cross with a fifth image in the middle of the cross. The HST F110W image shows the foreground group of four lensing galaxies (labeled G1, G2, G3, and G4) at $z \sim 1.0$—see Section 3.1.2. The $0.92 \times 0.52$ (PA = 23°) NOEMA beam is shown in the lower-left corner. The continuum contours are plotted starting at $5\sigma$ in steps of $3\sigma$ (the $1\sigma$ noise level is 69 µJy beam$^{-1}$). The NOEMA continuum map has been corrected for primary-beam attenuation.

Section 3.4). This is the first time that a gravitationally amplified high-$z$ galaxy is seen in an Einstein cross system with an image of the background galaxy at the center of the cross and that an Einstein cross is detected at submillimeter wavelengths.

### 3.1.2. *The Foreground Group of Galaxies*

The HST F110W image (at $1.15 \ \mu\text{m}$) shows that the $25'' \times 25''$ field around HerS-3 contains dozens of galaxies with various levels of brightness (Appendix A and Figure 15). The near-infrared emission is dominated by the three close by and brightest galaxies in the field, which are aligned along the

minor axis of the Einstein cross, with one galaxy to the southeast (G1) and the two others to the northwest (G2 and G3) of the central image (Figure 1). The galaxies G2 and G3, which are heavily blended with their diffuse emission overlapping, are each separated from G1 by $2.7''$. All three galaxies show faint, extended diffuse emission in the HST image. Based on the photometry available from the HST, SDSS, and Subaru/HSC data, the photometric redshifts of the galaxies G1, G2, and G3 are estimated to be in the range $0.97 < z_{\text{phot}} < 1.03$, indicating that they are part of a galaxy group. The galaxy G4, which is located $13.5''$ southeast of G1, also belongs to this group with a best-fitting photometric
displays a low median SNR computed for the unmasked pixels of 0.41 and was flagged as having a high fraction of points more than 5σ away from the best-fitting model. For these reasons, this spectroscopic redshift is only indicative, although we note that the value remains consistent with the results obtained from the SED analysis (Figure 17 in Appendix A).

Based on these findings, HerS-3 appears to be lensed by a group of at least four galaxies, G1, G2, G3, and G4, at $z_{phot} \sim 1.0$. By fitting the SEDs and modeling the surface brightness of each of these four lensing galaxies (Appendix A), we derived their photometric properties that are listed in Table 2. The galaxies G1, G2, G3, and G4 have low SFRs, $0.02 \lesssim \text{SFR}/(M_{\odot} \text{yr}^{-1}) \lesssim 1.5$, and stellar masses of $\sim 10^{11} M_{\odot}$, showing that they are massive quenched systems.

Investigating the field within a $\sim 60''$ ($\sim 500$ kpc) radius from G1 to search for potential additional members belonging to the lensing group, we found a total of 10 galaxies with photometric redshifts in the range $0.9 < z < 1.1$, which are consistent with those of the main lensing galaxies G1, G2, G3, and G4 (Appendix A and Figure 18). These other candidate members of the galaxy group have either low stellar masses or are located at far distances from the main mass distribution of the group. Therefore, these additional galaxies belonging to the group are unlikely to contribute significantly to the lensing of the background source, HerS-3.

### 3.2. Detailed View of the Einstein Cross

The $0''92 \times 0''52$ beam of the NOEMA data gives first indications that each of the images of the Einstein cross are extended in the east-west direction, in particular the NE and SW components (Figure 1). There is also evidence of changes in the orientation from one image to the other, again most clearly seen in the brightest NE and SW components.

These characteristics are revealed in greater detail in the higher angular resolution ($0''12 \times 0''10$) ALMA continuum Band 7 (1 mm) image, which provides crucial information on the spatial extent and orientation of the individual images of the Einstein cross (Figure 3). The centroids of all five images are well aligned with those of the 1 mm NOEMA image. The NE and SW images are elongated with extensions of $1''4$ and $1''2$, respectively, and widths of $0''3$, both showing a central emission peak, and orientations of PA $\sim 108^\circ$ and $\sim 126^\circ$, respectively. Of the three east-west oriented images, the central one is the strongest with a clear peak and emission extending over $\sim 0''6$ with a PA of $\sim 138^\circ$; the E component is weaker, extending over $\sim 0''8$ with a PA of $\sim 146^\circ$; and the W image is the weakest of all five components with an extension comparable to the central image.

The VLA 6 GHz continuum shows four of the five images (the weaker W component has only an upper limit), with characteristics that are comparable to those seen in the ALMA continuum image (Figure 4).

Table 1 lists the centroids (derived from the ALMA image) and the continuum flux densities of each of the images of the HerS-3 Einstein cross as measured by NOEMA, ALMA, and the VLA.

### 3.3. Continuum Spectral Energy Distribution

Together with the Herschel and z-GAL continuum flux densities (D. Ismail et al. 2023), the new NOEMA and ALMA data in the range between 240 and 290 GHz enable for a

continuous sampling of the far-infrared and submillimeter SED of HerS-3. In addition, the VLA continuum data at 6 GHz and 28 GHz (P. Prajapati et al. 2025, in preparation) provides crucial information on the origin of the radio emission in HerS-3. Figure 5 shows all the continuum data available for HerS-3 from infrared to radio wavelengths, and Table 3 lists the continuum flux densities.

We fit the continuum data using two different SED dust models, namely a single-temperature MBB in the optically thin approximation and an MBB model in its general form (S. Berta et al. 2021; D. Ismail et al. 2023, and references therein). The effect of the CMB, which produces a steepening of the millimeter slope, has been taken into account following the prescriptions of E. da Cunha et al. (2013). We adopt the $\kappa_\nu$ renormalization by B. T. Draine et al. (2014). The best-fit solutions are shown in Figure 5. The radio continuum data point is compared, a posteriori, to a power-law spectrum $S_\nu \propto \nu^{-0.8}$, which represents the synchrotron emission normalized on the basis of the radio/far-infrared correlation (B. Magnelli et al. 2015; J. Delhaize et al. 2017). The continuum at 6 GHz lies slightly below the expected star formation synchrotron power law based on the J. Delhaize et al. (2017) relation. There is therefore no indication of a radio-loud active galactic nucleus (AGN) in HerS-3. The 28 GHz continuum flux density, on the other hand, is in excess with respect to the J. Delhaize et al. (2017) synchrotron emission. To account for this, we add a free-free component in the form $S_\nu \propto \nu^{-0.1}$.

When adopting an optically thin (or general form, in parenthesis) MBB to model the dust continuum of HerS-3, we derive values for the following properties, where $\mu$ is the gravitational lens magnification: the spectral index is $\beta = 2.31 \pm 0.08$ ($2.67 \pm 0.16$); the dust temperature $T_{dust} = 34.1 \pm 1.1$ K ($60.0 \pm 0.08$ K); the dust mass $\mu M_{dust} = 7.43 \pm 0.28 \times 10^9 M_{\odot}$ ($5.08 \pm 0.96 \times 10^9 M_{\odot}$); and the far-infrared luminosity (from 50 to 1000 $\mu$m) $\mu L_{FIR} = 5.28 \pm 0.09 \times 10^{13} L_{\odot}$ ($5.31 \pm 0.11 \times 10^{13} L_{\odot}$), yielding an estimate of the total infrared luminosity (from 8 to 1000 $\mu$m) of $\mu L_{IR} = 7.54 \pm 0.12 \times 10^{13} L_{\odot}$ ($7.58 \pm 0.16 \times 10^{13} L_{\odot}$) by dividing $\mu L_{FIR}$ by 0.70 (S. Berta et al. 2023). These dust properties are consistent with the values that were derived by D. Ismail et al. (2023).

### 3.4. Molecular Emission and Absorption Lines

All five images of the Einstein cross of HerS-3 reveal a similarly rich spectrum, displaying the emission lines of $^{12}$CO(9–8) and H$_2$O($2_{20}-1_{11}$), and the OH$^+$($1_1-0_1$) and OH$^+$($1_2-0_1$) lines that are both seen in absorption. In addition, there is a tentative detection of NH$_2$ ($2_{20}-2_{11}$) in the global spectrum of HerS-3 (Figure 2).

We performed spectral fitting of the spectra in each individual image and the stacked global HerS-3 spectrum. Using the *curve_fit* procedure, we simultaneously fitted each of the molecular emission and absorption lines with a single Gaussian component and the underlying continuum with a linear curve. To reduce the number of free parameters and improve our sensitivity to the line parameters, we anchored the central velocity and line width of the two OH$^+$ transitions to one another. The same was done for the molecular emission lines, as they are expected to trace the same molecular gas reservoir. The best-fit spectra are plotted in Figure 2, and the best-fit parameters for the global spectrum are presented in Table 4.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The 291.8 GHz dust continuum emission of HerS-3 observed with ALMA in Band 7 provides a detailed view of the extent and orientation of each of the five images of the Einstein cross. The continuum contours are plotted starting at $3\sigma$ in steps of $2\sigma$ levels (with $1\sigma = 63 \mu$Jy beam$^{-1}$). The $0.12' \times 0.10' \ (PA = 74^\circ)$ beam is shown in the lower-left corner of the panel.</figcaption></figure>

The $^{12}$CO(9–8) spectra observed in each of the five images do not show any significant variation in their profile or central velocity, unambiguously demonstrating that the five images have the same redshift ($z = 3.0607$). The global spectra of the $^{12}$CO(9–8) and H$_2$O(2$_{02}-1_{11}$) emission lines have an FWHM of $320 \pm 28 \text{ km s}^{-1}$, and are, within the errors, at the expected redshifted frequencies (Table 4).

The peaks of both OH$^+$ absorption lines are blueshifted by ~350 km s$^{-1}$ with respect to the systemic velocity of the main molecular gas reservoir as traced by the CO emission lines (Table 4). Detected in absorption against the dust continuum, the gas seen in OH$^+$ is therefore interpreted as tracing material that is outflowing from the starburst along the line of sight.

The NOEMA data provide a first view of the molecular gas morphology and kinematics of HerS-3 as traced by the $^{12}$CO(9–8) emission line, which is the strongest line detected in the 1 mm spectra. The morphology of the molecular gas is shown in Figure 2, which displays the velocity-integrated moment-0 map of the $^{12}$CO(9–8) emission line. To investigate the velocity–position information of the molecular gas, we first binned the continuum-subtracted uv table presented in

Section 2.1 into channels of $\Delta v = 118 \text{ km s}^{-1}$. We then used the 2D uv plane source fitting procedure UV FIT to simultaneously fit 2D circular Gaussians to all five images in each velocity bin. This provides position and flux information of the $^{12}$CO(9–8) emission for each of the images in every velocity bin. The $^{12}$CO(9–8) emission spectra for all five images, extracted from the 2D Gaussian fitting, are shown in the left panel of Figure 6, which highlights their similarity in spectral shape. In the right panel of Figure 6, we show the shift in the central position of the $^{12}$CO(9–8) emission between the four highest-SNR velocity bins for each image. The arrows indicate the direction of the spatial shifts from one bin to the next, starting from 140 to –330 km s$^{-1}$, i.e., from red- to blueshifted emission.

In the NE and SW images, where the CO(9–8) emission is marginally resolved, the central position displays a general shift toward the southeast from red- to blueshifted velocities for both images. The high-resolution ALMA dust continuum map is shown underneath the arrows indicating the spatial drift in Figure 6. With this overlay, we see that the southeasterly shift in CO(9–8) emission from red to blue velocities follows
<figure><img src="image_6.png" /><figcaption>**Figure 4.** The 6 GHz radio continuum emission of HerS-3 observed with the VLA. The continuum contours are plotted starting at $2\sigma$ in steps of $1\sigma = 6~\mu\text{Jy}\text{ beam}^{-1}$. The peak 6 GHz continuum flux is $45~\mu\text{Jy}\text{ beam}^{-1}$. The $0.4'' \times 0.4''$ (PA = 130°) beam is shown in the lower-left corner.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** SED of HerS-3 from infrared to radio wavelengths, showing the contributions of the starburst, free-free, and star formation synchrotron emissions. The data include continuum flux densities from Herschel SPIRE (orange dots), ALMA (green dot), NOEMA (red dots), and VLA (pale blue dots)—the error bars correspond to $1\sigma$ uncertainties (Table 3). The fits are based on the modified blackbody (MBB) thin and general models (green and red dashed lines, respectively), including corrections for the effects of the cosmic microwave background (CMB; see text for details). The blue lines represent the synchrotron emission due to star formation, with a fixed spectral index $\alpha_{\text{synchr}} = -0.8$, normalized on the basis of the two radio/far-infrared correlations derived by J. Delhaize et al. (2017, lower line), and B. Magnelli et al. (2015, upper line). The dotted blue line shows the free-free emission, with a spectral index $\alpha_{\text{ff}} = -0.1$, normalized to fit the 28 GHz data after subtracting the J. Delhaize et al. (2017) synchrotron emission. The solid black line is the sum of the dust, the synchrotron (J. Delhaize et al. 2017), and the free-free components in the radio regime.</figcaption></figure>

the major axis of the galaxy's dust continuum emission. This may indicate rotation in the galaxy's disk with the northwest side of the galaxy rotating away from the observer and the southeast side moving toward the observer. High-SNR spectral line data at higher spatial resolution is needed to confirm this result and further investigate the kinematics of HerS-3.

<table><caption>**Table 3**

Far-infrared, Submillimeter, and Radio Continuum Flux Densities of HerS-3</caption><thead><tr><th>λ (mm)</th><th>ν (GHz)</th><th>S<sub>ν</sub> (mJy)</th><th>Reference</th></tr></thead><tbody><tr><td></td><td></td><td>Herschel</td><td></td></tr><tr><td>0.25</td><td></td><td>253 ± 6</td><td>(1)</td></tr><tr><td>0.35</td><td></td><td>250 ± 6</td><td>(1)</td></tr><tr><td>0.50</td><td></td><td>191 ± 7</td><td>(1)</td></tr><tr><td></td><td></td><td>ALMA</td><td></td></tr><tr><td>1.02</td><td>292</td><td>21.2 ± 1.4</td><td>(2)</td></tr><tr><td></td><td></td><td>NOEMA</td><td></td></tr><tr><td>1.16</td><td>258</td><td>18.6 ± 0.1</td><td>(2)</td></tr><tr><td>1.24</td><td>241</td><td>14.2 ± 0.1</td><td>(2)</td></tr><tr><td>1.94</td><td>154</td><td>3.1 ± 0.5</td><td>(3)</td></tr><tr><td>2.15</td><td>139</td><td>2.2 ± 0.4</td><td>(3)</td></tr><tr><td>2.91</td><td>103</td><td>0.30 ± 0.07</td><td>(3)</td></tr><tr><td>3.15</td><td>95</td><td>&lt;0.39</td><td>(3)</td></tr><tr><td>3.44</td><td>87</td><td>0.43 ± 0.08</td><td>(3)</td></tr><tr><td>3.79</td><td>79</td><td>&lt;0.21</td><td>(3)</td></tr><tr><td></td><td></td><td>VLA</td><td></td></tr><tr><td>11</td><td>28</td><td>0.10 ± 0.03</td><td>(4)</td></tr><tr><td>50</td><td>6</td><td>0.19 ± 0.02</td><td>(2)</td></tr></tbody></table>

**Note.** References: (1) H. Nayyeri et al. (2016); (2) This paper (3); D. Ismail et al. (2023); and (4) P. Prajapati et al. (2025, in preparation). The wavelengths and frequencies are the observed ones. The flux densities listed are the measured values with $1\sigma$ uncertainties. In the case of a nondetection, the flux density is given as a $3\sigma$ upper limit.

<table><caption>**Table 4**

Molecular Absorption and Emission Line Properties of HerS-3, Obtained by Line Fitting of the Global Spectrum</caption><thead><tr><th></th><th>S (Jy km s<sup>-1</sup>)</th><th>V (km s<sup>-1</sup>)</th><th>FWHM (km s<sup>-1</sup>)</th><th>Notes</th></tr></thead><tbody><tr><td colspan="5">Absorption Lines</td></tr><tr><td>OH<sup>+</sup>(1<sub>1</sub>-0<sub>1</sub>)</td><td>-3.78 ± 0.81</td><td>-350 ± 32</td><td>321 ± 78</td><td></td></tr><tr><td>OH<sup>+</sup>(1<sub>2</sub>-0<sub>1</sub>)</td><td>-2.87 ± 0.72</td><td>"</td><td>"</td><td>(1)</td></tr><tr><td colspan="5">Emission Lines</td></tr><tr><td><sup>12</sup>CO(9-8)</td><td>11.3 ± 0.95</td><td>-9.3 ± 12</td><td>320 ± 28</td><td></td></tr><tr><td>H<sub>2</sub>O(2<sub>02</sub>-1<sub>11</sub>)</td><td>5.7 ± 0.87</td><td>"</td><td>"</td><td>(2)</td></tr><tr><td>NH<sub>2</sub> (2<sub>20</sub>-2<sub>11</sub>)</td><td>2.3 ± 0.84</td><td>"</td><td>"</td><td>(2)</td></tr></tbody></table>

**Note.** The values listed are the primary-beam corrected line flux (S), the velocity shift (V) with respect to the redshift of HerS-3 ($z = 3.0607$) and the FWHM as derived from the line fitting of the global spectrum (Figure 2). (1) The absorption depth, central velocity, and line width of the two $OH^+$ transitions are anchored to one another. (2) The same procedure is applied to the central velocity and line width of the emission lines but their fluxes are left as free parameters. The uncertainties are taken from those given by the fitting procedure in *curve_fit* and the $NH_2$ emission is considered to be a tentative detection.

## 4. Gravitational Lensing Model

In order to recover the source-plane morphology and the intrinsic properties of HerS-3, and to explore the lensing configuration, we derived a lens model using the high angular resolution ALMA continuum data, with initial lens parameters informed by the HST, Subaru/HSC, and SDSS data. We performed the lens modeling and source reconstruction using the *lensmodel* package (C. R. Keeton 2001) with the pixel-based source reconstruction code *pixsrc* (A. S. Tagore &
<figure><img src="image_8.png" /><figcaption>**Figure 6.** 2D uv plane fitting results of the ¹²CO(9–8) emission line in bins of 118 km s⁻¹ for the five HerS-3 images. Left: ¹²CO(9–8) emission line flux as a function of velocity offset from the systemic velocity of HerS-3. Each of the five images are highlighted by a different shade of gray, revealing their ¹²CO(9–8) spectra. Right: movement between velocity bins of the ¹²CO(9–8) emission flux centers of each image overlaid on the ALMA 291.8 GHz dust continuum emission map where the contours are plotted starting at 3σ in steps of 2σ (with 1σ = 63 µJy beam⁻¹). The colored arrows indicate the direction of movement from bin to bin within 140 and −330 km s⁻¹, corresponding to the velocity ranges shaded in the same colors in the left panel. The NOEMA and ALMA beams are shown together in the lower-left corner and the axes centered on the NOEMA phase center as in Figure 2.</figcaption></figure>

C. R. Keeton 2014; A. S. Tagore & N. Jackson 2016; see Appendix B.2).

For an efficient and thorough analysis of the lensing configuration of the foreground group of galaxies, we conducted the lens modeling in a two-stage approach. First, we treated the five images of the Einstein cross as point sources, consistent with their compact sizes and the fact that they are well separated, to assess the constraints imposed by the positions of the lens images on possible solutions and identify the most promising lens models. We took into account the properties of the group of four lensing galaxies, such as their redshifts and shape priors (see Section 3.1.2 and Appendix B), assuming that the galaxies G1, G2, G3, and G4 represent the four mass components of the lensing configuration. Second, after identifying the most promising model, we reconstructed the full extended source structure and compared the predicted size, orientation, and brightness of each individual image to what is observed in the dust continuum emission with ALMA. Details of the lens modeling setup, including lists of parameters, priors, and constraints, are provided in Appendix B.1.

As described in Appendix B.2, the lens models that only include galaxies G1, G2, G3, and G4 are unable to reproduce the lensing properties observed for HerS-3. With mild priors on the shapes of the model galaxies, the model can match the lensed image positions only by making the model galaxy G3 much more elongated than observed (Figure 19(a)). Even then, the model cannot match the orientation of the extended structure in the SW image (Figure 20). If we impose stronger priors on the galaxy shapes, the model cannot reproduce the observed image positions, with offsets up to 0.′′2 (Figures 19(b) and 21). While we show results assuming that all four galaxies lie at the same redshift ($z_{phot} \approx 1$), the

conclusions hold for models where the foreground galaxies have different redshifts.

The failure of those lens models, and the fact that there is no other close by galaxy at the same redshift in the field, led us to consider the possibility of an additional fifth mass component that could represent a group dark matter halo. We try two models for this additional mass component. The first is a cored isothermal sphere with density profile $\rho(r) = \sigma^2/(2\pi G(s^2 + r^2))$ where s is the core radius, which corresponds to a scaled surface mass density for lensing of $\kappa(R) = \Sigma(R)/\Sigma_{crit} = b/(2\sqrt{s^2 + R^2})$ (see Appendix B.1 for more details about the parameters). The second model is a spherical Navarro–Frenk–White (NFW; J. F. Navarro et al. 1996) halo with density profile $\rho(r) = \rho_s r_s^3/(r(r_s + r)^2)$ where $r_s$ is the scale radius and $\rho_s$ is the scale density; the lensing properties of this model are given by M. Bartelmann (1996). The additional mass component is assumed to lie at the same redshift as the four visible galaxies. Details of this new lens model are given in Appendix B.3.

With the additional mass component, the lens model and observed configuration are in excellent agreement, with an rms offset between the observed and model brightness peaks of better than 0.′′02. Figure 7 depicts the model with a cored halo by showing the surface mass density map, the lensing time-delay surface, and the lensing critical curves and caustics. (The model uncertainties are listed in Table 11 in Appendix B.3.) The location of the additional component is constrained to lie 3′′–5′′ to the southeast of galaxy G1. The NE and SW images both lie outside the lensing critical curve and at minima of the time-delay surface, implying that they have positive parity. The three components E, C, and W all lie inside the critical curve and at saddle points of the time-delay surface, so they are expected to have negative parity. These results have
<figure><img src="image_9.png" /><figcaption>**Figure 7.** The best single-plane lens model with a cored halo as the additional mass component. All mass components are assumed to be at the same redshift ($z_{phot} \approx 1$). (a) The gray scale shows the surface density map; the galaxies G1, G2, and G3 are labeled (G4 is out of the frame to left); and the center of the additional mass component is marked with a yellow ★. The green + symbols show the observed positions of the five images of the Einstein cross, while the red × symbols show the positions of the images predicted by the lens model; the observed and predicted positions are in excellent agreement. The lensing critical curves are shown in cyan. (b) The color map shows the lensing time-delay surface, assuming a lens redshift $z_l = 1$ along with the measured source redshift $z = 3.0607$. The lensed images appear at stationary points of this surface; NE and SW lie at local minima, while E, C, and W lie at saddle points. The lensing critical curves are shown again for completeness. (c) The lensing caustics in the source plane; note that the spatial scale differs from the other panels. The inferred source position is indicated with the red × symbol.</figcaption></figure>

implications on the kinematics of the images as discussed in Section 5.1.

Strictly speaking, the lens model predicts two “central” images corresponding to local maxima of the time-delay surface, in addition to the five observed images. However, these “central” images are extremely demagnified: they could easily be $10^4$ time fainter than the NE image, with specific

values depending on the assumed central densities of the foreground lensing galaxies.

The full pixel-based source reconstruction shown in Figure 8 provides further support for this lens model. The model image configuration matches the observations in terms of the peak emission, brightness distribution, shape, and orientation for each of the five lensed images.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Results of the pixel-based source reconstruction of HerS-3 using the lens model with a cored halo from Figure 7. Top left: the ALMA 292 GHz continuum image, masked around the five lensed images. To guide the eye, green circles mark the locations of the peak brightnesses in the observed map. Top right: model image, lensed and convolved with the beam. Green circles are the same as in the data image. Bottom left: reconstructed source-plane image. The scale bar marks 1 kpc in the source plane. The green ellipse has been visually estimated to represent the extent of the brightness distribution; it has a semimajor axis of 0.15 or 1.2 kpc in the source plane. Bottom right: data minus model residuals. Black polygons indicate the masks used for the *pixsrc* analysis. Yellow lines in the color bar represent the noise level of 63 µJy beam⁻¹ measured outside of the masked regions.</figcaption></figure>

Figures 9 and 10 show the corresponding results for a lens model with an NFW halo for the additional mass components. Here the rms offset between the observed and model brightness peaks is again better than 0.02. Many of the features are qualitatively similar to the cored halo model: the halo is centered 3"–7" to the southeast of galaxy G1; the lensing critical curves extend from the southeast to the northwest; and the NE and SW images are again local minima of the time-delay surface while E, C, and W are all saddle points. The pixel-based source reconstruction again matches all the properties observed in the ALMA 292 GHz data. Broadly speaking, then, we conclude that the lensing configuration indicates the presence of a mass component to the southeast, but it cannot precisely constrain the internal structure of this structure.

We briefly explored two ways to generalize the cored halo models. First, we allowed the additional dark matter halo to be elliptical rather than spherical. There is a mild preference (~1.5σ) for this halo to be slightly elongated in the north–south direction; the other parameters, including those we use to draw conclusions, did not change substantially. Second, we treated galaxies G1, G2, and G3 with power-law mass profiles that were allowed to deviate from an isothermal profile. There

is a moderate preference (~2σ) for G2 to have a density profile that is steeper than isothermal; again the main lens modeling conclusions were largely unchanged. We defer a more thorough examination of complex lens models to future work.

We can use the lens models to predict various properties of the system. First, we can compute the projected lensing masses of the model galaxies, including both the baryon and dark matter components, within an aperture $R$ to be $M_{\text{lens}} = \pi b R \Sigma_{\text{crit}}$, where $b$ is the Einstein radius and $\Sigma_{\text{crit}}$ is the critical surface density for lensing. We use $R = 3''$ to match the aperture used to compute stellar masses. Table 5 shows that both the cored and NFW halo lens models predict lensing masses around $10^{12} M_{\odot}$ for galaxies G1 and G4, and slightly lower values for G2 and G3, all of which are reasonable in comparison with the stellar masses.

We can similarly compute the projected lensing mass of the large dark matter halo in each model. Figure 11 shows the enclosed mass curve, $M(r)$ as a function of $r$, for both the cored halo and NFW halo models. Over the range of radii spanned by the lensed images, the enclosed dark matter mass ranges from $\log M_{\text{DM}}/M_{\odot} = 12.2$ up to 13.0. In both models the enclosed mass curve continues to grow; while these larger radii are not well constrained by the lensing data, it is
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Similar to Figure 7, but for the lens model with an NFW halo for the additional mass component.</figcaption></figure>

interesting to see that the cored and NFW halo lens model yield similar results. Because the cored halo model has an isothermal profile, it is straightforward to convert the Einstein radius parameter into an effective velocity dispersion of $\sigma = 633_{-129}^{+185}$ km s⁻¹ (see Appendix B.1 for the relation). However, that value may not be especially enlightening as it depends on radii larger than the range spanned by the lensed images. There is no analogously simple quantity for an NFW halo.

We can also compute the lensing magnification factors, which are given in Table 6. The total magnification of all images is around 19 for the cored halo model and 17 for the NFW halo model, with some nonnegligible statistical

uncertainties. There is a degeneracy such that lens models can make the halo more massive and the galaxies less massive, or vice versa, over some reasonable range. When the halo is more massive, it contributes a higher surface mass density at the positions of the lensed images, which leads to higher magnification factors. Figure 12 illustrates this degeneracy by showing the relation between the total lensing magnification and the lensing mass of galaxy G1.

Finally, we can predict the lensing time delays, which are listed, in days, for both the cored halo and the NFW halo models in Table 7. In most models, the leading image is the NE one, with the SW image arriving a few days later, in between $7.9_{-11.3}^{+12.2}$ and $2.6_{-15.4}^{+14.6}$ days for the cored halo and NFW
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Similar to Figure 8, but for the lens model with the NFW halo from Figure 9. Here the ellipse in the source panel has a semimajor axis of 0.19 or 1.5 kpc.</figcaption></figure>

<table><caption>**Table 5**
Predicted Lensing Masses</caption><thead><tr><th rowspan="3">Galaxy</th><th colspan="2">Cored Halo Model</th><th colspan="2">NFW Halo Model</th></tr><tr><th>Best</th><th>MCMC</th><th>Best</th><th>MCMC</th></tr><tr><th colspan="2">log(M<sub>lens</sub>)</th><th colspan="2">log(M<sub>☉</sub>)</th></tr></thead><tbody><tr><td>G1</td><td>11.85</td><td>11.99<sup>+0.08</sup><sub>-0.15</sub></td><td>12.01</td><td>12.02<sup>+0.06</sup><sub>-0.06</sub></td></tr><tr><td>G2</td><td>11.53</td><td>11.65<sup>+0.11</sup><sub>-0.13</sub></td><td>11.67</td><td>11.68<sup>+0.08</sup><sub>-0.08</sub></td></tr><tr><td>G3</td><td>11.57</td><td>11.65<sup>+0.09</sup><sub>-0.12</sub></td><td>11.68</td><td>11.67<sup>+0.08</sup><sub>-0.08</sub></td></tr><tr><td>G4</td><td>11.86</td><td>11.98<sup>+0.13</sup><sub>-0.15</sub></td><td>12.02</td><td>12.05<sup>+0.12</sup><sub>-0.11</sub></td></tr></tbody></table>

**Note.** For each type of halo model, we report the projected lensing mass (including both baryon and dark matter components) within 3" (assuming a lens redshift $z_l = 1$). Uncertainties indicate the 68% confidence interval from Markov Chain Monte Carlo (MCMC) sampling of the posterior. For comparison, the stellar masses derived from the optical and near-infrared data are given in Table 2.

halo models, respectively. It should be noted that in some cases, the SW image arrives ahead of the NE image, resulting in negative time delays. For the three images E, C, and W, which are located at saddle points of the time-delay surface (Figures 7 and 9), the time delays relative to the NE image are predicted to be much longer, amounting to about 1 yr with the C, W, and E images in order of arrival time.

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** The curves show the projected enclosed dark matter halo mass, $M(r)$, as a function of $r$, for the two lens models using a cored isothermal sphere or a spherical NFW halo (blue and orange lines, respectively). Points mark the locations of the lensed images (relative to the halo center), to indicate the range of radii probed by the lensing data. Between the two models, the points have different $r$ values because the halos have different central positions.</figcaption></figure>

## 5. Discussion

In this section, we first present the general properties of the DSFG HerS-3 (Section 5.1) derived from all the currently available submillimeter to radio data and based on the results
<table><caption>**Table 6**
Predicted Lensing Magnifications</caption><thead><tr><th rowspan="2">Image</th><th colspan="3">Cored Halo Model</th><th colspan="3">NFW Halo Model</th></tr><tr><th colspan="3">Eigenvalues</th><th colspan="3">Eigenvalues</th></tr><tr><th></th><th></th><th></th><th>Magnification</th><th></th><th></th><th>Magnification</th></tr></thead><tbody><tr><td>SW</td><td>2.74<sup>+0.91</sup><sub>-0.38</sub></td><td>2.16<sup>+0.67</sup><sub>-0.31</sub></td><td>5.91<sup>+4.47</sup><sub>-1.51</sub></td><td>2.63<sup>+0.39</sup><sub>-0.32</sub></td><td>2.08<sup>+0.28</sup><sub>-0.23</sub></td><td>5.46<sup>+1.66</sup><sub>-1.20</sub></td></tr><tr><td>NE</td><td>3.81<sup>+0.23</sup><sub>-0.37</sub></td><td>1.99<sup>+0.68</sup><sub>-0.31</sub></td><td>7.64<sup>+5.67</sup><sub>-5.21</sub></td><td>3.68<sup>+0.57</sup><sub>-0.42</sub></td><td>1.90<sup>+0.27</sup><sub>-0.25</sub></td><td>7.00<sup>+2.16</sup><sub>-1.64</sub></td></tr><tr><td>E</td><td>-0.88<sup>+0.15</sup><sub>-0.32</sub></td><td>1.77<sup>+0.75</sup><sub>-0.42</sub></td><td>-1.55<sup>+0.52</sup><sub>-1.39</sub></td><td>-0.86<sup>+0.15</sup><sub>-0.19</sub></td><td>1.43<sup>+0.24</sup><sub>-0.20</sub></td><td>-1.24<sup>+0.30</sup><sub>-0.44</sub></td></tr><tr><td>C</td><td>-1.36<sup>+0.26</sup><sub>-0.47</sub></td><td>1.46<sup>+0.53</sup><sub>-0.24</sub></td><td>-1.97<sup>+0.61</sup><sub>-1.68</sub></td><td>-1.29<sup>+0.21</sup><sub>-0.25</sub></td><td>1.34<sup>+0.19</sup><sub>-0.15</sub></td><td>-1.74<sup>+0.44</sup><sub>-0.53</sub></td></tr><tr><td>W</td><td>-1.17<sup>+0.31</sup><sub>-0.44</sub></td><td>1.57<sup>+0.43</sup><sub>-0.18</sub></td><td>-1.85<sup>+0.63</sup><sub>-1.32</sub></td><td>-1.07<sup>+0.21</sup><sub>-0.24</sub></td><td>1.55<sup>+0.22</sup><sub>-0.17</sub></td><td>-1.66<sup>+0.43</sup><sub>-0.53</sub></td></tr><tr><td>Total</td><td></td><td></td><td>19.00<sup>+14.66</sup><sub>-5.44</sub></td><td></td><td></td><td>17.06<sup>+5.21</sup><sub>-4.01</sub></td></tr></tbody></table>

**Note.** For each type of halo model, we report the two eigenvalues of the lensing magnification tensor (which indicate the linear magnification factors in the principal directions; a negative value indicates parity reversal) and the overall magnification factor (the product of the eigenvalues). Uncertainties indicate 68% confidence intervals from MCMC sampling of the lens models.

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Relationship between the projected lensing mass of galaxy G1 (horizontal axis) and the total lensing magnification (vertical axis) for a selection of lens models drawn from the posterior parameter distributions for the cored halo models and NFW halo models.</figcaption></figure>

of the lensing model described in Section 4. In Section 5.2, we review the characteristics of the lensing galaxy group placing the findings of this study in a broader context, in particular, by outlining what future observations are needed to further constrain the properties of the massive dark matter halo, and in Section 5.3, we describe the potential use of HerS-3 to measure the Hubble constant.

### 5.1. The Starburst HerS-3

In the source plane, HerS-3 appears as a nearly edge-on galaxy (Figures 8 and 10) with a bright central region surrounded by two weaker emission peaks. The size of the galaxy (major and minor axis), estimated from the 1 mm dust continuum emission, is $0''30 \times 0''12$ (~2.3 $\times$ 0.9 kpc$^2$) or $0''38 \times 0''15$ (~2.9 $\times$ 1.2 kpc$^2$) for the models with the cored halo and the NFW profile, respectively, with PA ~ 120°. Based on the results of the source-plane reconstruction, we can derive the inclination angle from the minor to major axis ratio $b/a$ of the disk as $i = \cos^{-1}(b/a)$, yielding $i \sim 66^\circ$ for both the cored halo and NFW profile models. The rotation in the galaxy's disk, traced in the $^{12}$CO(9-8) emission line (Figure 6), indicates a gradient from red to blueshifted gas

along the west-east axis, which is best seen in the NE and SW images that have the highest-SNR data. This is consistent with the predictions of the lensing model as illustrated in the left panel of Figure 13. Higher angular resolution observations are required to further explore the kinematics of this starburst and reveal the details of the velocity field, in particular in the E, C, and W images. In addition, the blueshifted absorption lines of $OH^+$, which are detected in each of the five images, indicate a strong outflow activity with gas expelled from the starburst with velocities of at least ~350 km s$^{-1}$.

HerS-3 appears therefore to be a nuclear starburst with a morphology that is comparable to the nearby galaxy NGC 253, known to possess a starburst-driven wind (A. D. Bolatto et al. 2013; J. P. Pérez-Beaupuits et al. 2018). As in NGC 253, the outflow axis in HerS-3 is most likely perpendicular to the active nucleus, specifically along the northeast-southwest axis where the parity reversal is occurring (middle panel of Figure 13). Observing this outflow at higher angular resolution, using molecular absorption lines (such as $OH^+$ or $OH$) or atomic emission lines (e.g., [C II] 158 µm), will allow us to explore in greater detail the morphology of the wind activity in HerS-3, akin to what has been achieved in other high-$z$ active star-forming galaxies (e.g., K. M. Butler et al. 2021; J. S. SPIEKER et al. 2020a, 2020b; R. Herrera-Camus et al. 2021), and to trace the parity reversal expected in each of the three saddle-point images, providing a basis to further refine the lensing model.

<table><caption>**Table 7**
Predicted Lensing Time Delays</caption><thead><tr><th rowspan="3">Image</th><th>Cored Halo Model</th><th>NFW Halo Model</th></tr><tr><th colspan="2">Δt</th></tr><tr><th colspan="2">(days)</th></tr></thead><tbody><tr><td>SW</td><td>7.9<sup>+12.2</sup><sub>-11.3</sub></td><td>2.6<sup>+14.6</sup><sub>-15.4</sub></td></tr><tr><td>E</td><td>391.1<sup>+76.2</sup><sub>-99.9</sub></td><td>403.8<sup>+65.7</sup><sub>-54.0</sub></td></tr><tr><td>C</td><td>273.6<sup>+53.3</sup><sub>-71.2</sub></td><td>275.9<sup>+50.3</sup><sub>-38.2</sub></td></tr><tr><td>W</td><td>344.5<sup>+68.4</sup><sub>-85.6</sub></td><td>349.7<sup>+57.2</sup><sub>-44.3</sub></td></tr></tbody></table>

**Note.** Predicted time delays relative to the NE image, assuming a lens redshift $z_l = 1$. In most models, NE is the leading image, but, in some cases, the SW image arrives slightly ahead of NE and those correspond to negative time delays.
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Left: model lensed velocity maps assuming a simple rotating disk galaxy with a flat rotation curve and a PA of 120° based on our source reconstruction. The insets, displaying each image in detail, are 0.5 across. The observed changes in the ¹²CO(9–8) velocity distributions from image to image (as presented in Figure 6) are compatible with the model predictions. Middle: velocity map for a model of outflowing gas assumed to be perpendicular to the disk, showing that the axis along which the parity is flipped runs from upper left to lower right. The lensing critical curves are shown in black. Right: in this model, the outflowing gas is predicted to cross two lensing caustics, creating an additional redshifted filament between G2 and G3 in the velocity map. We assume here that the redshifted outflow lies to the north and the blueshifted outflow lies to the south, but the situation could be reversed depending on the galaxy’s inclination. This figure shows results for the cored halo lens model, but the results for the NFW halo model are very similar.</figcaption></figure>

### 5.1.1. Molecular Gas Excitation and Properties

To investigate the CO line excitation and estimate the physical conditions of the molecular gas in HerS-3, we modeled the CO line fluxes using a large velocity gradient (LVG) treatment of statistical equilibrium and radiative transfer that is described in Appendix C. In total, we have measurements of four ¹²CO transitions for HerS-3, namely: $J=9-8$ (11.3 ± 0.95 Jy km s⁻¹; Table 4), $J=5-4$ and $3-2$ (8.27 ± 0.68 and 5.25 ± 0.36 Jy km s⁻¹, respectively; P. Cox et al. 2023), and $J=1-0$ (0.71 ± 0.08 Jy km s⁻¹; P. Prajapati et al. 2025, in preparation). The velocity-integrated line fluxes here used are from the global spectrum of HerS-3 and have not been corrected for gravitational magnification. As the CO spectral line energy distribution (SLED) of HerS-3 displays a strong ¹²CO(9–8) emission line, we used a two-excitation component model, assigning different parameters to each of the components. The two-component best-fitting solution of the CO SLED of HerS-3 is shown in Figure 14. In Appendix C, the posterior probability distributions of the molecular gas density $n_{H_2}$, the gas temperature $T_{kin}$, and the CO column density per velocity $N_{CO}/dv$ of the source for the cold and warm components are presented in Figure 23. The fit reproduces the overall SLED very well, explaining the observed ¹²CO(3–2), (5–4), and (9–8) relative line fluxes, which a single component model could not have achieved. The physical properties derived from the model are listed in Table 8.

The molecular gas in HerS-3 is significantly excited with conditions similar to those of the Cloverleaf (A. Weiß et al. 2003; C. M. Bradford et al. 2009; D. A. Riechers et al. 2011), HFLS-3 (D. A. Riechers et al. 2013), and NGC 6240 (R. Meijerink et al. 2013). Excitation through the hard radiation (including X-rays) associated with a luminous AGN is likely not at play as there is no indication for an AGN in HerS-3 (see Section 3.3). It is therefore more plausible that the excitation of the dense molecular gas in HerS-3 is dominated by infrared photons and collisions (shocks) originating in the massive, intense starburst and the outflowing gas traced in the OH⁺ absorption lines. This is in line with the CO luminosity to infrared luminosity ratio of $L_{CO}/L_{IR} \sim 1.6 \times 10^{-4}$ in HerS-3, where $L_{CO}$ has been

extrapolated, up to $J=11-10$, from the CO SLED (Figure 14). This ratio is higher than for Mrk 231 ($7 \times 10^{-5}$; P. P. van der Werf et al. 2010), where the excitation is dominated by X-rays from the AGN, and the maximum ratio of $L_{CO}/L_{IR} \sim 10^{-4}$ predicted for models of photon and X-ray dominated regions (R. Meijerink et al. 2013, and references therein). As in NGC 6240, where the ratio is $7 \times 10^{-4}$ (R. Meijerink et al. 2013), the excitation of the interstellar medium in HerS-3 is therefore most likely dominated by shocks, which heat the gas and not the dust. Measuring additional, intermediate $J$ level CO transitions would be useful to further constrain the nature of the excitation of the molecular gas in HerS-3.

The ¹²CO(1–0) emission line flux enables us to estimate the line luminosity, $L'_{CO}(1-0)$, and the molecular gas mass, $M_{H_2}$, associated with HerS-3, without the uncertainties related to the conversion when using higher-$J$ CO lines. We derive the values of $L'_{CO}(1-0)$ and $M_{H_2}$ using the following equations (e.g., C. L. Carilli & F. Walter 2013):

$$L'_{CO} = 3.25 \times 10^7 \frac{S_{CO} \Delta v_{CO} D_L^2}{\nu_{CO,rest}^2 (1+z)} \text{ K km s}^{-1} \text{ pc}^{-2}, \quad (1)$$

$$M_{H2} = \alpha_{CO} \times L'_{CO}(1-0) M_{\odot}, \quad (2)$$

where $S_{CO} \Delta v_{CO}$ is the measured flux density of the line in units of Jy km s⁻¹, $D_L$ the luminosity distance in Mpc units, $\nu_{CO,rest}^2$ the rest frequency of the line in units of GHz, and $\alpha_{CO}$ is the ¹²CO(1–0)-to-$H_2$ conversion factor. We here adopt a value of $\alpha_{CO} = 3 M_{\odot}$ (K km s⁻¹ pc²)⁻¹ as derived by L. Dunne et al. (2021), and we multiply it by the usual factor of 1.36 to account for the helium contribution. The ¹²CO(1–0) line luminosity is $\mu L'_{CO} = 3.07 \pm 0.35 \times 10^{11}$ K km s⁻¹ pc⁻², and the molecular gas mass $\mu M_{H_2} = 1.25 \pm 0.14 \times 10^{12}$ M$_{\odot}$. The resulting molecular gas to dust mass ratio is $\delta_{GDR} \sim 168$ for the optically thin MBB dust mass (and 246, in the case of the general MBB), which remains consistent with the mean ratio of $128^{+54}_{-35}$ derived by L. Dunne et al. (2021) for a sample of Herschel-selected starburst galaxies (see also L. Dunne et al. 2022; S. Berta et al. 2023; F. Stanley et al. 2023).

From these properties, we can estimate the star formation efficiency (SFE) of HerS-3, which is a measure of the SFR per
<figure><img src="image_17.png" /><figcaption>**Figure 14.** LVG two-excitation component model fit to the ¹²CO SLED of Hers-3. The velocity-integrated line fluxes (not corrected for gravitational magnification) and their errors are shown as black dots and the solid orange curve shows the best fit from the two-component (MCMC) model corresponding to the maximum posterior possibility. The orange dashed and red dashed-dotted lines show the warmer and cooler component fits, respectively. The derived physical properties of the molecular gas of HerS-3 are listed in Table 8.</figcaption></figure>

<table><caption>**Table 8**

Physical Properties of the Molecular Gas of HerS-3 Derived from the ¹²CO Spectral Line Energy Distribution Two-component Large Velocity Gradient Model</caption><thead><tr><th>Component</th><th>log(n<sub>H2</sub>)<br>(cm<sup>-3</sup>)</th><th>log(T<sub>kin</sub>)<br>(K)</th><th>log<sub>10</sub>(N<sub>CO/dv</sub>)<br>(cm<sup>-2</sup> km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Low-excitation</td><td>3.08<sup>+0.92</sup><sub>-1.08</sub></td><td>1.21<sup>+0.16</sup><sub>-0.09</sub></td><td>17.22<sup>+0.95</sup><sub>-0.73</sub></td></tr><tr><td>High-excitation</td><td>3.22<sup>+0.89</sup><sub>-0.74</sub></td><td>2.33<sup>+0.25</sup><sub>-0.26</sub></td><td>18.24<sup>+0.78</sup><sub>-0.69</sub></td></tr></tbody></table>

**Note.** The values listed are the gas density ($n_{H2}$), the kinetic temperature ($T_{kin}$), and the CO column density ($N_{CO/dv}$), where the errors are $3\sigma$, as derived from the two-component LVG model.

unit gas, and the gas depletion time, $t_{depl}$, which is the inverse of the SFE, indicating the time it would take to deplete the gas reservoir with a constant SFR and no other ongoing processes, such as gas accretion or outflows. The SFE can be expressed as $SFE = SFR/M_{H2} \text{ yr}^{-1}$, where SFR ($M_{\odot} \text{ yr}^{-1}$) = $1.09 \times 10^{-10} L_{IR} (L_{\odot})$ and $L_{IR}$ is the rest-frame (8–1000 $\mu$m) infrared luminosity (R. C. J. Kennicutt 1989; corrected for a G. Chabrier 2003 initial mass function). For HerS-3, using the infrared luminosity derived in Section 3.3, we estimate $\mu$SFR = $8.2 \times 10^3 M_{\odot} \text{ yr}^{-1}$ and SFE = $6.5 \times 10^{-9} \text{ yr}^{-1}$, yielding a depletion time of $t_{depl}$ = 152 Myr. This depletion time is significantly shorter than typical main-sequence galaxies and is instead consistent with infrared luminous starburst galaxies at high $z$ (see, e.g., S. Berta et al. 2023; F. Stanley et al. 2023, and references therein).

### 5.1.2. Intrinsic Properties of HerS-3

Depending on the type of dark matter halo profile adopted in the lensing model, the overall gravitational magnification factor changes from $\mu = 19_{-5}^{+14}$ for the cored halo model to $\mu = 17_{-4}^{+5}$ for the NFW halo model (Table 6). As explained in Section 4, this difference is due to the degeneracy between the halo mass and the galaxy masses and deeper observations at higher angular resolution of both the lensing group and HerS-3

will be needed to further constrain the lensing properties and, hence, reduce the significant errors on the magnification factor. Adopting a value of $\mu \sim 18$, based on the current estimates, we derive the following intrinsic properties of HerS-3: $L_{IR} \sim 4.2 \times 10^{12} L_{\odot}$, for the total infrared luminosity; SFR $\sim 455 M_{\odot} \text{ yr}^{-1}$, for the SFR; $M_{H2} \sim 6.9 \times 10^{10} M_{\odot}$, for the molecular gas mass; and $M_{dust} \sim 4.1 \times 10^8 M_{\odot}$ (or $2.8 \times 10^8 M_{\odot}$, in the case of the general MBB), for the dust mass.

### 5.2. Lensing Group and Dark Matter Halo

The observed properties of the five images of the HerS-3 Einstein cross can be precisely matched by a lensing model that includes, in addition to the group of four visible galaxies at $z_{phot} \sim 1$, a fifth massive component at the same redshift, located $\sim 3''$ to $7''$ (i.e., $\sim 25$ to 60 kpc) southeast of the brightest galaxy G1 (Section 4). The fact that no galaxy is seen at that position in the HST F110W image points to the presence of a massive dark matter halo associated with the lensing group. Heuristically, the location to the southeast of this massive component is indicated by the curvature of the arc defined by the NE, C, and SW images. As noted above, the need for an additional halo seems robust, but conclusions about its properties depend on assumptions, such as the radial profile (which we have shown explicitly by comparing cored and NFW halo models) and our simplifying approach that the halo is spherical. Even so, it is worthwhile to examine what we can infer about the system from the current data and lens models.

Since we assumed in the lensing model that the dark matter halo is at the same redshift as the visible galaxies, we also tested the effect of placing the halo at different redshifts. Taking the best models of each class (cored and NFW halo), we explored the redshift range $z = 0.6$ to 1.2, and reoptimized all of the model parameters. For both models, only minor changes resulted at the different redshifts in the core density and the mass of the dark matter halo.

As reported in Appendix A, the four visible galaxies G1, G2, G3, and G4 lensing HerS-3 are part of a larger group containing 10 more detected members at $z_{phot} \sim 1$ within a $1'$ (500 kpc) radius from the brightest galaxy G1. The stellar mass of the four galaxies located at the center of the group is estimated to be $M_* \simeq 3.6 \times 10^{11} M_{\odot}$ (Table 2), accounting for 60% of total stellar mass of the group (Appendix A). Compared to their total (baryonic and dark matter) predicted by the lensing model of $M_{tot} = 3 \times 10^{12} M_{\odot}$ (Table 5), the ratio of the stellar to total mass for the four galaxies G1–G4 is $M_*/M_{tot} \sim 12\%$.

Galaxy groups are small aggregates of galaxies, comprising several to $\sim 50$ galaxies within $\sim 1$ Mpc, that are embedded within their dark matter halos. They are less massive than galaxy clusters, which have masses up to a few $\times 10^{13} M_{\odot}$, and have thus smaller lensing cross sections than clusters. Being more numerous than clusters, groups of galaxies are known to be potential powerful systems to probe the properties and structure of dark matter halos (e.g., D. C. Fox & U.-L. Pen 2001). The total mass of the group, containing the galaxies lensing HerS-3, is in line with the median mass of known galaxy groups, $M_{group} \sim 3 \times 10^{12} M_{\odot}$ (e.g., M. Einasto et al. 2024). The lensing group is also a relatively compact group with a radius of $\approx 500$ kpc, where the four most massive galaxies (G1, G2, G3, and G4) are quenched systems with low SFRs ($0.02 \le SFR/(M_{\odot} \text{ yr}^{-1}) \le 1.5$; see Section 3.1.2).
It is interesting to note that the position of the massive dark matter halo component, needed to explain the observed properties of the HerS-3 Einstein cross, is located in between the galaxies G1–G3 and G4, at a position approximately corresponding to the center of mass of the entire group. From the estimated lensing mass of the dark matter halo, $log(M_{DM}/M_\odot) \sim 13.0$, using the mass enclosed within the largest sampled radius (Figure 11), the $z \sim 1.0$ galaxy group lensing HerS-3 follows the stellar-to-halo mass relation derived in the $log(M_{DM}/M_*) - log(M_{DM})$ plane for $z \sim 1$ galaxies with $log(M_*/M_{DM}) \le -1.2$ and $log(M_{DM}/M_\odot) \sim 13$ (see Figure 6 in G. Girelli et al. 2020).

Based on the current available observational data, the lensing model, which successfully reproduces the properties of the five images of the HerS-3 Einstein cross, reveals a massive dark matter halo at the center of a galaxy group, with no counterpart seen in the HST near-infrared image. The four brightest lensing galaxies G1–G4 are distributed around the dark matter halo at distances ranging from ~25 to 80 kpc. All other members of the group are located at much larger distances (see Appendix A). This configuration is reminiscent of a nonrelaxed galaxy group, suggesting that the group at $z_{phot} = 1$ lensing HerS-3 is not virialized. However, additional observations will be needed to further constrain the lensing model and, henceforth, the properties of both the massive dark matter halo and the galaxy group, so as to better understand its evolutionary state.

Finally, as an example of further refinements in the analysis of the lensing configuration, our lens models indicate that the center of the background source, HerS-3, lies 0.6–0.8 kpc outside a lensing caustic (see panel (c) in each of Figures 7 and 9). If the source exhibits an outflow to the northeast that is ≳1 kpc in extent, as expected based on the outflow activity traced by the $OH^+$ absorption (Section 5.1), that emission would cross two lensing caustics, creating an additional lensed arc-like image in between the galaxies G2 and G3 (see Figure 13). This is a clear prediction of current lens models that could directly be tested with future high angular resolution observations of molecular or atomic tracers of the outflow activity in HerS-3, either with the James Webb Space Telescope (JWST) via atomic emission lines to trace the extended ionized gas of the outflow, via $H_2$ rovibrational lines to probe the hot turbulent gas, or using ALMA to detect the high-velocity emission wings of the [C II] emission line. Observing this arc-like image and deriving its characteristics will provide the basis to refine the lens model and constrain the properties of the dark matter halo by distinguishing between the core halo or NFW profile.

### 5.3. HerS-3 as a Cosmological Probe

HerS-3 provides a great cosmological probe for measuring gravitational time delays and derive the value of the Hubble constant (see, e.g., P. L. Kelly et al. 2023; T. Schmidt et al. 2023). The active star-forming activity of this galaxy makes the event of a supernova explosion possible or could result in intrinsic starburst variability that could be followed to measure the time delay of the brightness variations between the five observed images in the plane of the lens. Both lensing models predict lensing time delays between the five images of the Einstein cross in HerS-3 that are comparable, in particular for the E, C, and W images, and within a factor of 3 for the SW image relative to the NE one. Whereas the delays are in the

range of a few days for the NE and SW images (with arrival times that can swap between them), the E, C, and W images have delays that amount to about a year, with the central image predicted to appear first, with a time delay of about 275 days relative to the NE image. Long-term observations and monitoring will therefore be needed to follow the changes and the corresponding time delays in the E and W images, once an event or indications of variability have been detected in the NW or SE image.

## 6. Conclusions

We have presented results based on NOEMA, ALMA, and VLA observations of the DSFG HerS-3 at $z=3.0607$, which was found to be gravitationally amplified in an Einstein cross with a fifth bright image at the center. The exceptional lensing configuration of HerS-3 has never been seen before and represents the first detection of an Einstein cross at submillimeter wavelengths. HerS-3 is lensed by four galaxies located at the center of a larger galaxy group at $z_{phot} \sim 1.0$, seen in the HST F110W image, and a massive dark matter halo, which is inferred from the lensing model in order to reproduce all the properties of the five images of the Einstein cross. Our main findings are the following.

1. All five images of the Einstein cross of HerS-3 have been detected with NOEMA in various molecular lines, including the $^{12}$CO(9–8) and $H_2O(2_0-1_{11})$ emission lines, which have similar line profiles and central velocities, and the $OH^+(1_1-0_1)$ and $OH^+(1_2-0_1)$ lines seen in absorption, unequivocally indicating that they share similar redshifts and are therefore the lensed images of HerS-3. The images have wide separations, extending to 7"/5 for the NE and SW images, which are unusually large compared to other Einstein crosses. The high angular resolution (0"/1) 1 mm dust continuum measured with ALMA reveals that each image is extended in the east-west direction with PAs and brightnesses changing from one image to the other. The NOEMA and ALMA data are complemented with VLA imaging data tracing the radio continuum at 6 GHz.

2. At a redshift of $z_{phot} \sim 1.0$, the lensing group is composed of three, close by, massive ($M_* \sim 10^{11} M_\odot$), quenched systems (G1, G2, and G3), with a fourth galaxy (G4) 13"/5 to the southeast of G1 with a stellar mass comparable to G1. The Einstein cross is centered in between the three close by galaxies of the group. These lensing galaxies belong to a larger group with 10 more members within an ~500 kpc radius around G1.

3. Taking only into account the four visible massive galaxies G1–G4 that are close to HerS-3 and located at the center of the galaxy group, the lensing models were unable to reproduce the properties of the five images of the Einstein cross observed in the high angular resolution ALMA continuum image. The fact that there is no other visible galaxy at the same redshift in the nearby field around the lensing group (the other galaxies of the group, which are far away and much less massive, are unlikely to contribute significantly to the lensing of HerS-3), led us to consider the possibility of a massive additional component that could be the group dark matter halo. Only by adding this massive component, which is constrained to lie ~3"–7" (or 25–60 kpc) to the southeast
of G1, does the source reconstruction match the properties of the five images, namely their position, brightness, spatial extent, and orientation.

4. Lensing models adopting a cored isothermal sphere and a spherical NFW halo for the dark matter halo yield equally good results for the source reconstruction. However, the internal structure of the dark matter halo cannot be constrained with the current data. In particular, the distribution of the surface density map is different for each of the dark matter halo profiles. As a result, the halo mass may have any value from $log(M_{DM}/M_{\odot}) = 12.2$ to 13.0, and the total lensing magnification is in between $\mu \sim 17$ and 19 (with significant uncertainties) for the NFW and cored halo profiles, respectively.

5. In the source plane, HerS-3 appears as a starburst galaxy with a highly inclined rotating disk traced in $^{12}$CO(9–8). Blueshifted absorption lines of $OH^+$ are detected in each image, revealing a strong outflow activity with gas expelled at velocities of at least 350 km s$^{-1}$ from the nuclear region. The molecular gas reservoir of HerS-3, with an estimated mass of $M_{H_2} \sim 6.9 \times 10^{10} M_{\odot}$, has an average density of $n_{H_2} \sim 3000 \text{ cm}^{-3}$ and high-excitation conditions with a kinetic temperature $T_{kin} \sim 270$ K, likely dominated by shocks. With a total infrared luminosity of $L_{IR} \sim 4.2 \times 10^{12} L_{\odot}$ and an SFE $= 6.5 \times 10^{-9} \text{ yr}^{-1}$, the depletion timescale is $t_{depl} = 152$ Myr, placing HerS-3 in the regime of infrared luminous starbursts.

The remarkable lensing configuration of HerS-3 and the results that have been obtained to date provide a foundation to study in greater depth the properties of this system. Future higher angular resolution and deeper observations, with facilities such as ALMA in the submillimeter and the JWST in the near and mid-infrared, will enable us to further explore the morphology and the kinematics of HerS-3 down to spatial scales of $\le 100$ pc in the source plane, and to obtain up to 3 orders of magnitude better near-infrared images than the HST of the lensing galaxy group and the environment of the dark matter halo, which could reveal weaker components of the galaxy group and additional lensing effects. The prospect of such high-quality data opens the way for further improvements in the lensing model, including testing observational predictions related to the parity reversal in the E, C, and W images of the Einstein cross and the effect expected if the powerful wind of HerS-3 does cross two lensing caustics. Henceforth, it will be possible to better constrain the mass distribution of the galaxy group lensing HerS-3 and understand its evolutionary state, derive with greater precision the properties of the dark matter halo, such as its profile and associated mass, place stringent limits on any associated stellar emission, and, eventually, show evidence of substructures.

The HerS-3 system with its exceptional Einstein cross with a fifth central image has been revealed by the present study to be a unique astrophysical laboratory to explore, at small spatial scales, a nearly edge-on dusty starburst galaxy at the peak of cosmic evolution and, importantly, to study the characteristics of the $z_{phot} \sim 1.0$ galaxy group lensing HerS-3 and the properties of the associated massive dark matter halo.

## Acknowledgments

This work is based on observations carried out under project number W21DF with the IRAM NOEMA Interferometer.

IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). This paper makes use of the following ALMA data: ADS/JAO.ALMA 2022.1.00145.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. This work is also based on observations carried out under program VLA/22A-211 using the National Radio Astronomy Observatory’s (NRAO) Karl G. Jansky Very Large Array (VLA). The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The anonymous referee is thanked for providing useful comments that helped to improve the contents of this paper. R. Eliot and D. Herbey are kindly thanked for their initial help with the IRAM data reduction. We would also like to thank L.-D. Roger for useful discussions about some aspects of this study. This work benefited from the support of the project Z-GAL ANR-AAPG2019 of the French National Research Agency (ANR). E.B. and E.M.C. are supported by Padua University grants Dotazione Ordinaria Ricerca (DOR) 2022–2024 and by the Istituto Nazionale di Astrofisica (INAF) grant Progetto di Ricerca di Interesse Nazionale (PRIN) 2022 2022383WFT “SUNRISE” (CUP C53D23000850006). T.J.L. C.B. gratefully acknowledges the financial support from the Knut and Alice Wallenberg foundation through grant no. KAW 2020.0081. A.J.B. acknowledges support from the Radcliffe Institute for Advanced Study at Harvard University. D.R. gratefully acknowledges support from the Collaborative Research Center 1601 (SFB 1601 subprojects C1, C2, C3, and C6) funded by the Deutsche Forschungsgemeinschaft (DFG)-500700252. L.M. acknowledges financial support from the Inter-University Institute for Data Intensive Astronomy (IDIA), a partnership of the University of Cape Town, the University of Pretoria, and the University of the Western Cape, and from the South African Department of Science and Innovation’s National Research Foundation under the ISARP RADIOMAP Joint Research Scheme (DSI-NRF grant No. 150551) and the CPRR Projects (DSI-NRF grant No. SRUG2204254729).

## Appendix A

### Photometric Redshift of the Foreground Group of Lensing Galaxies

In this section, we derive the characteristics of the foreground group of lensing galaxies from the near-infrared HST F110W image (Figure 15) and the available SDSS and Subaru/HSC photometry. Extracting the lens’ photometry is necessary to characterize the properties of each member of the galaxy group, such as their fluxes, surface brightness distributions, SEDs, and photometric redshifts. However, this extraction is complicated when sources are heavily blended, making simple aperture photometry techniques not straightforward (e.g., galaxies G2 and G3 in Figure 15). To address this issue, we applied surface brightness modeling techniques to fully disentangle the 2D emission of each galaxy.

Some key ingredients are necessary to model the surface brightness distribution of a galaxy: these are the characteristic PSF of the observations, a noise map, and a source mask covering all sources but those that will be modeled.
<figure><img src="image_18.png" /><figcaption>**Figure 15.** The HST F110W image showing the group of the four massive lensing galaxies (G1, G2, G3, and G4) at $z_{phot} \sim 1$.</figcaption></figure>

Additionally, given the large number of sources present in an HST-like field of view (FOV), it is useful to have a segmentation image where each source detected above a given threshold is labeled with a nonzero integer. This segmentation image will not be used to compute photometry but to estimate the blending between sources and the sources' centroids to be able to automatically perform the fitting. We detected sources up to ~60" (~500 kpc) from the main lensing group of G1, G2, and G3. For the HST image, the PSF model and noise map are available from E. Borsato et al. (2024), and for the Subaru data these were both downloaded from the archive. Whereas for the SDSS *u* band, we only downloaded the available noise map. In order to compute the missing input ingredients, we first produced the segmentation image. We stacked the Subaru/HSC data computing a mean Subaru image and its associated noise map. We convolved the stacked images and noise map with a small Gaussian kernel to reduce the pixel-to-pixel variations ensuring well defined detection segments. We used the *detect.sourcesphotutils* task to identify all sources having at least 10 connected pixels above SNR = 2. We decided to base the detection on the Subaru data and not on the HST data in order (i) to ensure a reasonable detection threshold across most of the available bands and (ii) to avoid

including in the segmentation map all those HST sources that are too faint to have a Subaru/SDSS counterpart and are thus not usable for the SED fitting. This kind of detection approach severely suffers from source blending, i.e., when two sources are close enough to be connected above the threshold level and are identified as one source. To partially address this issue, we applied the *deblend.sourcesphotutils* task that uses multithresholding and watershed segmentation to separate blended source segments. The resulting separation is sharp and does not include contamination between the deblended sources. In order to make full use of the HST image quality, we first projected the Subaru-based segmentation image to the HST scale and then performed the deblending on it. To produce the source mask, we first enlarged each deblended segment and then computed their intersections. When two or more sources intersected, we defined a group mask as the merged enlarged intersecting segments. With the groups defined, we computed the source centers in the HST image through *photutilssource catalogue* given the HST-deblended segmentation image. Lastly, we made cutouts of each group around their group mask, masking all sources that were partially included in the cutout. We then projected each cutout and mask from the HST image to the other bands. With
the source mask computed, we measured and subtracted the background using the photutils class background2D. This class estimates the background map using (sigma-clipped) statistics in grid boxes covering the input image. The background is estimated as the median value computed in the masked and sigma-clipped box region. We adopted 30 × 30 pixel boxes. By combining the boxes, it is possible to have a 2D low-resolution map of the image background. Then, the final background map is calculated by interpolating the low-resolution background map.

We performed the surface brightness modeling using PyAutoGalaxy (J. Nightingale et al. 2023). For each group, the procedure adopted for the fitting is as follows: We first modeled the HST data; when the group was made up of multiple galaxies, we performed a first preparatory fit, where we masked all but one of the blended sources fitting it with a Sérsic profile iteratively. We used the results of these fits to initialize the simultaneous fit of all group galaxies by taking Gaussian priors centered on the maximum likelihood values of the common parameters and with a width corresponding to that of the posterior distribution scaled to allow an adequate sampling of the parameter space (see Section 6.1 in J. W. Nightingale et al. 2018, for details).

With the HST fit performed, we moved to fitting all the other bands. To simplify the fitting procedure, we now assumed that the galaxies' relative positions, flattening, and PAs are the same in all bands and equal to the maximum likelihood HST results. However, the effective surface brightness, effective radius, and Sérsic index were left free but given narrow priors as before. We introduced two new parameters to model possible residual astrometric offsets between the HST image and each other given band. This approach seemed to work well for most cases, while significantly reducing the computing time and complexity of the model. As a concluding step, we inspected all the surface brightness fits, editing the model and refitting the data where significant residuals were found, similarly to what was done in E. Borsato et al. (2024). We found that it was enough to add unresolved components in the HST band for the more complex sources. The substructures that required the addition of complexity in the HST image were in most cases washed away by the lower resolution and sensitivity of both the Subaru and SDSS data.

Notably, some galaxies were undetected in the SDSS u and Subaru/HSC g bands. We randomly placed 3" apertures in the SDSS FOV to compute upper limits in these bands. The images were masked to remove the source's emission. Then, we computed the flux distribution through these apertures and estimated the standard deviation, weighting each flux by the number of unmasked pixels in the aperture. To measure the fluxes of the lenses, we integrated the surface brightness profiles up to a circularized radius of 3". This was done for two reasons: first, to provide a straightforward comparison to the upper limits in the SDSS u and SDSS g bands and, second, to limit the flux from the tails of the surface brightness model of the southeast lens in the HST data. The 3" aperture was large enough to include all significant model contributions given the constraints possible with the data available for a galaxy at $z \sim 1$. If the flux estimate is compatible with zero, i.e., $F_d - 3*\sigma_{F_d} \le 0$, we adopt as the upper limit the 95th percentile of the flux distribution. In Figure 16, we show the results of the surface brightness modeling for the four lensing galaxies G1, G2, G3, and G4, done in the SDSS u filter, the

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Results of the surface brightness modeling for the lensing galaxies G1, G2, G3 (upper panel), and G4 (lower panel) in the SDSS u, Subaru/HSC g, r, i, z, and Y, and HST/WFC3 F110W filters. In each row, we show the same 5"/6 × 5"/6 cutouts of the observed data, best-fitting models, and residuals. The data, models, and residuals are all normalized to the corresponding noise map. The color map is defined to be in linear gray scale where the SNR is between -3 and +3, and in logarithmic scale from an SNR of three to the maximum observed in the cutout otherwise. The red dashed regions mark the portions of the data that have been masked during the fit and that were therefore not used in the minimization process.</figcaption></figure>

Subaru/HSC g, r, i, z, and Y filters, and the HST/WFC3 F110W filter. In Table 9, we report the flux densities of the G1, G2, G3, and G4 lensing galaxies derived from the SDSS, Subaru/HSC, and HST observations.

We fit the lenses' SEDs using the bagpipes code<fnref n="25" />(A. C. Carnall et al. 2018). We parameterized the star formation history of each lens with an exponential burst, and we included a variable dust absorption following the Calzetti absorption law (D. Calzetti 2001), as well as a variable stellar metallicity and redshift. The wavelength coverage of our observations and their broadband nature do not allow us to constrain the main group lenses' (G1, G2, G3, and G4) stellar metallicity and dust absorption well. This effect produced broad posterior distributions reaching very high and unlikely metallicities ($\gtrsim 2 Z_\odot$) and dust absorptions ($\gtrsim 1$ mag). These fits generally produced photometric redshift posteriors $\sim 1$, low SFRs of $0.01 \lesssim \text{SFR}/(M_\odot \text{yr}^{-1}) \lesssim 1$, and stellar masses of $\sim 10^{11} M_\odot$, confirming that these galaxies are massive (at least in their stellar component) quenched systems. Due to degenerations between the SED parameters in the fit, a better constraint of the stellar metallicities and dust absorptions is necessary to avoid affecting the photometric redshift solution. Hence, as priors for the metallicities and dust absorptions, we adopted the values measured by A. C. Carnall et al. (2022) for a sample of 114 $z \sim 1$ passive and massive galaxies with $\log_{10}(M_*/M_\odot) > 10.8$ that matched reasonably well with our lensing galaxies. These priors force the metallicity (Z) and dust extinction ($A_V$) of our lenses to span reasonable values. Under these assumptions, the median photometric redshifts we found are $z_{G1} = 1.03_{-0.03}^{+0.03}$, $z_{G2} = 0.99_{-0.05}^{+0.04}$, $z_{G3} = 0.97_{-0.03}^{+0.04}$, and $z_{G4} = 1.06_{-0.01}^{+0.01}$. The uncertainty intervals are computed from the 16th and 84th percentiles of the full posterior distribution.

<footnote n="25">²⁵ Available at https://github.com/ACCarnall/bagpipes.git</footnote>
<table><caption>**Table 9**
Flux Densities of the G1, G2, G3, and G4 Lensing Galaxies Derived from the SDSS, Subaru/HSC, and HST Observations</caption><thead><tr><th>Lens</th><th>SDSS <i>u</i></th><th>Subaru <i>g</i></th><th>Subaru <i>r</i></th><th>Subaru <i>i</i></th><th>Subaru <i>z</i></th><th>Subaru <i>Y</i></th><th>HST F110W</th></tr></thead><tbody><tr><td>G1</td><td>&lt;1.21</td><td>0.62 ± 0.07</td><td>1.81 ± 0.11</td><td>5.53 ± 0.15</td><td>11.96 ± 0.27</td><td>16.65 ± 0.45</td><td>21.21 ± 0.06</td></tr><tr><td>G2</td><td>&lt;1.21</td><td>0.18 ± 0.04</td><td>0.64 ± 0.07</td><td>2.06 ± 0.11</td><td>4.28 ± 0.24</td><td>7.08 ± 0.35</td><td>8.17 ± 0.03</td></tr><tr><td>G3</td><td>&lt;1.21</td><td>&lt;0.30</td><td>0.77 ± 0.08</td><td>2.51 ± 0.14</td><td>5.56 ± 0.28</td><td>8.87 ± 0.43</td><td>9.85 ± 0.03</td></tr><tr><td>G4</td><td>&lt;1.21</td><td>0.42 ± 0.05</td><td>2.69 ± 0.09</td><td>7.99 ± 0.13</td><td>17.23 ± 0.25</td><td>25.86 ± 0.46</td><td>30.37 ± 0.10</td></tr></tbody></table>

Note. All the flux densities and upper limits are in units of µJy. The upper limits are computed at 3σ.

All galaxies have redshift distributions that are contained within $z=0.9$ and $z=1.1$. Given the differences between the group members and possible residual systematics related to the 2D deblending and SED modeling, we assume the group to be at a redshift of $z_{phot} \sim 1$. The best-fitting SEDs are shown in Figure 17, and the properties of the galaxies are listed in Table 2.

For the remaining galaxies within a radius of ~60" (~500 kpc) from the main lensing group, given that many of the neighboring sources were fainter, we performed three additional SED fits for each galaxy, where we reduced the model complexity by fixing the metallicity to subsolar, solar, and supersolar values, i.e., $Z/Z_{⊙} = 0.5$, $Z/Z_{⊙} = 1$, and $Z/Z_{⊙} = 1.5$, respectively. Again, if all the best-fitting SEDs predicted a photometric redshift $z_{phot} \sim 1$ and low specific SFRs, we used the A. C. Carnall et al. (2022) templates as above. To decide which best-fitting SED model to use, we prioritized the A. C. Carnall et al. (2022) SED when available, otherwise we used the all-parameters-free models, unless

either the metallicity or the dust extinction reached extreme values. Then, we selected as candidate group members all the galaxies that had the chosen best-fitting median photometric redshift between 0.9 and 1.1. This resulted in a total of 14 candidate members, including G1, G2, G3, and G4, with a total stellar mass of $M_{*} \sim 6 \times 10^{11} M_{⊙}$. A summary plot of the galaxy group with all the candidate members is shown in Figure 18. It is important to note that most of the stellar mass of this galaxy group at $z_{phot} \sim 1$ is accounted for by the four galaxies G1, G2, G3, and G4, which contain ~60% of the total mass within 500 kpc and ~90% within 250 kpc (where only two additional candidate group members were identified). These findings suggest that the lensing galaxies G1–G4 belong to a more extended and complex galaxy group. This result is in line with the lensing modeling analysis, which indicates the presence of a massive dark matter halo, which is located in between G1 and G4 and is close to the mass center of the group, in order to reproduce the observed properties of the HerS-3 Einstein cross (see Section 4).
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Best-fitting SEDs for the lensing galaxies G1, G2, G3, and G4 in blue, orange, green, and red, respectively (from top to bottom). In each panel, the SED is shown in black, while the 1σ and 3σ uncertainty intervals are highlighted with dashed and more transparent bands, respectively. The observed flux densities are displayed with filled black symbols when there is a detection, and hollow symbols when an upper limit is used. Each data set is represented by a different symbol (shown in the upper left box) with the triangles tracing the SDSS data, squares tracing the Subaru/HSC data, and pentagons showing the HST/WFC3 data. The posterior distributions of the photometric redshifts of the lensing galaxies are shown in separate insets in each of the corresponding panels.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Galaxies that are candidate group members in the field around G1, and have photometric redshifts in the range 0.9 < z<sub>phot</sub> < 1.1, are shown in orange and overplotted on the HST image. The size of each candidate member galaxy is normalized to its best-fitting stellar mass. The continuous red circle corresponds to a radius of 60'' (~500 kpc), while the dashed circle corresponds to a radius of 20'' (~150 kpc), with both circles centered on the position of the galaxy G1. The image is shown in gray scale for pixels having −3 < SNR < 5 and in viridis in the range 5 ≤ SNR < 100 in order to highlight the sources in the HST image.</figcaption></figure>

## Appendix B
### Details of the Lens Model

We conducted the lens modeling in two stages. First, because the ALMA images are well separated and fairly compact, we treated them as point images for initial parameter space exploration. This approach is computationally efficient and facilitates a thorough exploration. Once we identified promising models, we reconstructed the full extended source structure using the `pixsrc` software by A. S. Tagore & C. R. Keeton (2014) and A. S. Tagore & N. Jackson (2016).

### B.1. Lens Model Setup: Parameters and Constraints

For the point image analysis, we used the positions in Table 1 and computed offsets relative to galaxy G1 (see

coordinates in Table 2). We assumed position uncertainties of 0.''03. We used the ALMA 292 GHz flux densities (and their uncertainties) as constraints on the relative brightnesses.

We treat the galaxies G1, G2, G3, and G4 as singular isothermal ellipsoids. These models implicitly assume that the galaxies have dark matter halos that extend beyond the baryonic components, and the baryons and dark matter combine to have a total density profile of the form $\rho(r) = \sigma^2/(2\pi Gr^2)$, which corresponds to a scaled surface mass density for lensing of

$$ \kappa(R) = \frac{\Sigma(R)}{\Sigma_{\text{crit}}} = \frac{b}{2R} \quad \text{where} \quad b = 4\pi \left(\frac{\sigma}{c}\right)^2 \frac{D_{ls}}{D_s}. \quad (B1) $$
<figure><img src="image_22.png" /><figcaption>**Figure 19.** The gray scale shows surface mass density maps for lens models that only include galaxies G1, G2, G3, and G4, all assumed to lie at the same redshift ($z_{phot} \sim 1$). The green + symbols show the observed image positions, while the red × symbols show the positions of the images predicted by the lens model. Panel (a) shows results with mild shape priors; the model G3 is highly elongated, making this model implausible. Panel (b) shows results with strong shape priors; here the galaxy shapes are reasonable, but this model cannot reproduce the image positions well. In both panels, galaxy G4 is outside of the frame to the left.</figcaption></figure>

Here $\sigma$ is the velocity dispersion, $D_s$ and $D_{ls}$ are the angular diameter distances from the observer or lens to the source, respectively, and $\Sigma_{crit}$ is the critical surface mass density for lensing. We can generalize to an elliptical model by replacing $R$ with an elliptical radius.

Each isothermal mass component is described by five parameters:<fnref n="26" /> two position coordinates, the Einstein radius parameter $b$, an ellipticity, and a PA. We constrain the galaxy positions using the offsets given in Table 2 with uncertainties of 0.07/03. For a fixed aperture $R$, the lensing mass of an isothermal ellipsoid is $M_{lens} = \pi b R \Sigma_{crit}$, where $\Sigma_{crit}$ is the critical density for lensing. Therefore we place lognormal priors on the relative Einstein radii using the relative stellar masses (from Table 2), with uncertainties of 0.1 dex. We fit the ellipticity and PA using the quasi-Cartesian components $e_c = e \cos 2\theta_e$ and $e_s = e \sin 2\theta_e$. In initial modeling, we found that the model galaxies sometimes became highly elongated. In an attempt to limit such behavior, we imposed mild Gaussian priors on $e_c$ and $e_s$ using the observed ellipticities and PAs, with uncertainties of 0.20 in each component (see Table 10). Last but not least, we included external shear characterized by the quasi-Cartesian components $\gamma_c$ and $\gamma_s$ to represent tidal effects from other mass in the vicinity of the lens galaxies or along the line of sight. We imposed Gaussian priors on the shear of $\gamma_{c,s} = 0.00 \pm 0.05$ for mild priors and $0.00 \pm 0.02$ for strong priors.

Since the precise lens redshifts ($z_l$) are unknown, we initially assume that all four galaxies are at the same redshift. In that case the modeling can be done in angular units without assuming a specific value for the lens redshift, but we do need to assume $z_l$ in order to convert from angular units to physical units. We explored the possibility that the galaxies have

<table><caption>**Table 10**
Galaxy Shape Priors</caption><thead><tr><th>Galaxy</th><th>e<sub>c</sub></th><th>e<sub>s</sub></th><th>Standard priors</th><th>Strong priors</th></tr></thead><tbody><tr><td>G1</td><td>-0.04</td><td>-0.17</td><td>±0.20</td><td>±0.05</td></tr><tr><td>G2</td><td>0.09</td><td>-0.38</td><td>±0.20</td><td>±0.05</td></tr><tr><td>G3</td><td>0.05</td><td>0.57</td><td>±0.20</td><td>±0.05</td></tr></tbody></table>

**Note.** The quasi-Cartesian shape parameters are $e_c = e \cos 2\theta_e$ and $e_s = e \sin 2\theta_e$, where $e$ is the ellipticity and $\theta_e$ is the PA measured east of north. Note that G4 is assumed to be spherical. Columns 4 and 5 give the standard deviation of Gaussian priors on $e_c$ and $e_s$ in the standard and strong cases, respectively. We also include Gaussian priors on external shear of the form $\gamma_{c,s} = 0.00 \pm 0.05$ for standard priors, and $0.00 \pm 0.02$ for strong priors.

different redshifts; while some details changed, the overall conclusions remained the same, so we defer further consideration of multiplane lens models until more information is known about the redshifts.

To summarize, the minimal models have a total of 20 mass model parameters (18 for the four galaxies and two for external shear) and three source parameters (position and brightness). There are a total of 34 constraints or priors on the models: 15 constraints from the image positions and brightnesses, eight additional constraints from the galaxy positions, three priors on the relative Einstein radii, and eight shape priors.

## B.2. Models with Just Four Galaxies

We first consider lens models that include four mass components representing the four observed galaxies G1, G2, G3, and G4 (plus external shear). Figure 19 depicts the results by showing the surface mass density map as well as the observed and predicted image positions. The only way for this type of model to match the images reasonably well is to make

<footnote n="26">²⁶ We treated G4 as spherical for simplicity, so it has three parameters rather than five.</footnote>
<figure><img src="image_23.png" /><figcaption>**Figure 20.** Results of pixel-based source reconstruction of HerS-3 using the lens model with just the four main galaxies (G1, G2, G3, and G4) and standard priors (as shown in Figure 19(a)). Top left: the ALMA 292 GHz continuum image, masked around the five lensed images. To guide the eye, circles mark the image coordinates given in Table 1. Top right: model image, lensed and convolved with the beam. Green circles represent the same image peak brightness positions as in the data image. Bottom left: reconstructed source-plane image. The scale bar marks 1 kpc in the source plane. Bottom right: data minus model residuals. Black polygons indicate the masks used for the *pixsrc* analysis. Yellow lines in the color bar represent the noise level of 63 µJy beam⁻¹ measured outside of the masked regions.</figcaption></figure>

one of the model galaxies highly elongated (panel (a)). This model has a total $\chi^2 = 104$, with contributions of 12 from the image positions (corresponding to an rms offset in the image positions of 0.05), 40 from the image brightnesses, and 52 from the priors (compared to 11 degrees of freedom). Because the model G3 is unreasonably elongated, we consider strengthening the shape priors by reducing the uncertainties on the $e_c$ and $e_s$ shape components to 0.05. The resulting model has more plausible galaxy shapes but it cannot reproduce the image positions very well (panel (b)); this model has $\chi^2 = 300$, with 90 from the image positions (an rms offset of 0.13), 52 from the image brightnesses, and 158 from the priors.

Even though these models have clear shortcomings, it is still worthwhile to see how they handle the full ALMA 292 GHz data. We use the pixel-based source reconstruction code *pixsrc* from A. S. Tagore & C. R. Keeton (2014). This code reconstructs the source surface brightness on a grid in the source plane in a way that minimizes a penalty function with two terms. The first term is a goodness of fit comparing the model image (lensed and convolved with the beam) to the ALMA data. The second “regularization” term penalizes

deviations from a smooth surface brightness distribution to prevent the model from overfitting noise. Note that we apply *pixsrc* to the real space data (not the *uv* plane) and adopt the noise scaling technique of D. A. Riechers et al. (2008) to account for spatially correlated noise (see also C. E. Sharon et al. 2019).

We construct a data mask around each of the five observed images of the background galaxy, including a small margin of background in each region as required by *pixsrc*. Outside of the masked regions, we measure a noise level of 63 µJy beam⁻¹. With the noise scaling approach, we find that increasing the noise level in *pixsrc* to 80 µJy beam⁻¹ provides the best balance between goodness of fit and regularization. For the beam, we use an elliptical Gaussian with an FWHM of 0.12 and 0.10 along the major and minor axis, respectively, and a PA of 74°.

Figures 20 and 21 show results of pixelated source reconstructions for the models with mild and tight shape priors, respectively. The inability to match the observed positions is apparent. Moreover, the models fail to match the orientation of the extended structure in the SW image, as well as general features of the E, C, and W images.
<figure><img src="image_24.png" /><figcaption>**Figure 21.** Similar to Figure 20, but for the lens model only including the four main galaxies with strong priors (as shown in Figure 19(b)).</figcaption></figure>

### B.3. Models with an Additional Mass Component

In many lens systems with multiple galaxies in close proximity, the galaxies are part of a group or cluster (see, e.g., P. Natarajan et al. 2024; Appendix A). We therefore consider a new set of models with a fifth mass component that could represent a common dark matter halo. We treat the additional mass component as a cored isothermal sphere or a spherical NFW halo. In both cases the halo is described by four parameters (position, Einstein radius or density parameter, and scale radius), bringing the total number of mass model parameters to 27.

We introduce an approximate method to impose constraints from the extended image shapes without requiring a full source reconstruction. In the absence of beam effects, lensing conserves surface brightness, so image shapes are determined by spatial distortions. In the vicinity of an image *i*, distortions are governed by the local lensing magnification tensor $\mu_i$ such that small displacements in the image plane, $\Delta\mathbf{x}_i$, are related to

small displacements in the source plane, $\Delta\mathbf{u}$, by $\Delta\mathbf{x}_i = \mu_i\Delta\mathbf{u}$. Repeating the same analysis for a different image *j*, we see that we can describe the transformation from image *i* to image *j* as $\Delta\mathbf{x}_j = (\mu_j \mu_i^{-1})\Delta\mathbf{x}_i$. We can apply the idea to lens modeling as shown in Figure 22. Specifically, given a model we can compute the tensor $\mu_j \mu_{NE}^{-1}$ that would transform the NE image to each of the other image locations (middle row). We can then compare the transformed images to the observed images (top row) and use the residuals (bottom row) to compute a shape penalty for the lens models. This method is less precise than a full pixelated source reconstruction, in particular because it does not account for beam smearing, but it is much faster so it can be easily incorporated into the parameter space exploration.

Models with the additional mass component are quite successful, and are discussed fully in Section 4. For completeness, we list all of the parameter values here in Table 11.
<figure><img src="image_25.png" /><figcaption>**Figure 22.** Illustration of the lensed image shape analysis using linear transformations. The top row shows cutouts around the five lensed images in the ALMA 292 GHz map. The middle row shows transformed versions of the NE image computed using the best lens model with a cored halo. The bottom row shows residuals; such residuals are used to compute an image shape penalty when exploring the model parameter space. Note that this analysis does not properly account for beam effects, but it is still useful for assessing how well different models can match the image shapes.</figcaption></figure>

<table><caption>**Table 11**
Parameters for Lens Models with Four Galaxies and an Additional Dark Matter Halo</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="2">Cored Halo Model</th><th colspan="2">NFW Halo Model</th></tr><tr><th>Best</th><th>MCMC</th><th>Best</th><th>MCMC</th></tr></thead><tbody><tr><td>log b<sub>1</sub></td><td>-0.27</td><td>-0.13<sup>+0.08</sup><sub>-0.15</sub></td><td>-0.11</td><td>-0.09<sup>+0.06</sup><sub>-0.06</sub></td></tr><tr><td>x<sub>1</sub></td><td>0.00</td><td>0.00<sup>+0.03</sup><sub>-0.03</sub></td><td>0.00</td><td>0.01<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>y<sub>1</sub></td><td>0.02</td><td>0.02<sup>+0.03</sup><sub>-0.03</sub></td><td>0.00</td><td>0.01<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>e<sub>c1</sub></td><td>0.03</td><td>0.06<sup>+0.10</sup><sub>-0.10</sub></td><td>0.06</td><td>0.07<sup>+0.11</sup><sub>-0.10</sub></td></tr><tr><td>e<sub>s1</sub></td><td>-0.23</td><td>-0.19<sup>+0.15</sup><sub>-0.13</sub></td><td>-0.03</td><td>-0.03<sup>+0.16</sup><sub>-0.17</sub></td></tr><tr><td>log b<sub>2</sub></td><td>-0.58</td><td>-0.46<sup>+0.11</sup><sub>-0.13</sub></td><td>-0.44</td><td>-0.43<sup>+0.08</sup><sub>-0.08</sub></td></tr><tr><td>x<sub>2</sub></td><td>2.13</td><td>2.14<sup>+0.03</sup><sub>-0.03</sub></td><td>2.12</td><td>2.13<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>y<sub>2</sub></td><td>1.79</td><td>1.79<sup>+0.03</sup><sub>-0.03</sub></td><td>1.80</td><td>1.79<sup>+0.02</sup><sub>-0.03</sub></td></tr><tr><td>e<sub>c2</sub></td><td>0.03</td><td>0.02<sup>+0.13</sup><sub>-0.14</sub></td><td>-0.01</td><td>-0.01<sup>+0.13</sup><sub>-0.12</sub></td></tr><tr><td>e<sub>s2</sub></td><td>-0.36</td><td>-0.37<sup>+0.15</sup><sub>-0.13</sub></td><td>-0.37</td><td>-0.38<sup>+0.16</sup><sub>-0.13</sub></td></tr><tr><td>log b<sub>3</sub></td><td>-0.55</td><td>-0.47<sup>+0.09</sup><sub>-0.12</sub></td><td>-0.43</td><td>-0.44<sup>+0.08</sup><sub>-0.08</sub></td></tr><tr><td>x<sub>3</sub></td><td>2.52</td><td>2.51<sup>+0.03</sup><sub>-0.03</sub></td><td>2.52</td><td>2.51<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>y<sub>3</sub></td><td>1.13</td><td>1.13<sup>+0.03</sup><sub>-0.03</sub></td><td>1.13</td><td>1.14<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>e<sub>c3</sub></td><td>-0.14</td><td>-0.20<sup>+0.16</sup><sub>-0.13</sub></td><td>-0.24</td><td>-0.23<sup>+0.15</sup><sub>-0.12</sub></td></tr><tr><td>e<sub>s3</sub></td><td>0.60</td><td>0.55<sup>+0.09</sup><sub>-0.14</sub></td><td>0.49</td><td>0.52<sup>+0.10</sup><sub>-0.13</sub></td></tr><tr><td>log b<sub>4</sub></td><td>-0.25</td><td>-0.13<sup>+0.13</sup><sub>-0.15</sub></td><td>-0.09</td><td>-0.07<sup>+0.12</sup><sub>-0.11</sub></td></tr><tr><td>x<sub>4</sub></td><td>-6.33</td><td>-6.33<sup>+0.03</sup><sub>-0.03</sub></td><td>-6.33</td><td>-6.33<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>y<sub>4</sub></td><td>-12.33</td><td>-12.33<sup>+0.03</sup><sub>-0.03</sub></td><td>-12.34</td><td>-12.34<sup>+0.03</sup><sub>-0.03</sub></td></tr></tbody></table>
<table><caption>**Table 11**
(Continued)</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="2">Cored Halo Model</th><th colspan="2">NFW Halo Model</th></tr><tr><th>Best</th><th>MCMC</th><th>Best</th><th>MCMC</th></tr></thead><tbody><tr><td>log b<sub>h</sub></td><td>0.89</td><td>0.75<sup>+0.22</sup><sub>-0.20</sub></td><td>-0.91</td><td>-0.87<sup>+0.14</sup><sub>-0.09</sub></td></tr><tr><td>x<sub>h</sub></td><td>-2.36</td><td>-3.05<sup>+0.73</sup><sub>-1.06</sub></td><td>-2.92</td><td>-3.42<sup>+0.87</sup><sub>-1.87</sub></td></tr><tr><td>y<sub>h</sub></td><td>-1.97</td><td>-2.56<sup>+0.67</sup><sub>+1.00</sub></td><td>-2.38</td><td>-2.79<sup>+0.83</sup><sub>-1.73</sub></td></tr><tr><td>log s<sub>h</sub></td><td>0.70</td><td>0.57<sup>+0.21</sup><sub>-0.23</sub></td><td>1.52</td><td>1.47<sup>+0.28</sup><sub>-0.31</sub></td></tr><tr><td>γ<sub>c</sub></td><td>0.02</td><td>0.02<sup>+0.01</sup><sub>-0.01</sub></td><td>0.03</td><td>0.03<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td>γ<sub>s</sub></td><td>-0.06</td><td>-0.06<sup>+0.02</sup><sub>-0.02</sub></td><td>-0.07</td><td>-0.07<sup>+0.03</sup><sub>-0.03</sub></td></tr></tbody></table>

**Note.** We list the best value along with the median and 68% confidence interval for each parameter in each of the two types of lens models with a cored isothermal sphere halo and a spherical NFW halo. Note: in a high-dimensional parameter space with an asymmetric posterior distribution, the best value may not fall within the 68% confidence interval of the marginalized distribution for some parameters. All x- and y-positions are given in arcseconds relative to the observed position of galaxy G1. For galaxies G1–G4, *b* is the Einstein radius parameter in arcseconds. Galaxies G1–G3 are allowed to be elliptical, while G4 is assumed to be spherical. For the cored halo, *b*<sub>h</sub> is the Einstein radius parameter and *s*<sub>h</sub> is the core radius, both in arcseconds. For the NFW halo, *b*<sub>h</sub> is the NFW lensing strength $\kappa_s = \rho_s r_s / \Sigma_{crit}$ (which is dimensionless), and *s*<sub>h</sub> is the scale radius *r*<sub>s</sub> in arcseconds. *e*<sub>c</sub> and *e*<sub>s</sub> are the quasi-Cartesian components of ellipticity, while $\gamma_c$ and $\gamma_s$ are the analogous parameters for external shear (which are all dimensionless).

## Appendix C
Radiation Transfer Large Velocity Gradient Model

To investigate the CO line excitation of HerS-3 and estimate the physical conditions of the molecular gas, we derived an LVG radiative transfer model for the CO SLED of HerS-3 using the MCMC implementation of RADEX (F. F. S. van der Tak et al. 2007), as presented in C. Yang et al. (2017). This 1D non-LTE radiative transfer code assumes spherical symmetry and describes the velocity gradients, $dv/dr$, in the expanding sphere approximation. It performs an MCMC sampling of the parameter space that includes the molecular hydrogen density ($n_{H_2}$), the gas kinetic temperature ($T_k$), the CO column density per unit velocity gradient ($N_{CO}/dv$), and the solid angle of the source ($\Omega_{app}$), and it samples the posterior probability

distributions functions from the CO line fluxes modeled by RADEX. A Bayesian approach is used to fit the line fluxes from the model, and the code `emcee` is adopted to perform the MCMC calculation.

As the CO SLED of HerS-3 displays a strong $^{12}$CO(9–8) emission line (Figure 14), we used a two-excitation component model, assigning different parameters to each of the components and assuming two additional priors on the relative sizes and temperatures within the Bayesian analysis. The results of the two-component model are shown in Figure 14. Figure 23 displays the posterior distributions of the molecular gas density $n_{H_2}$, the gas temperature $T_{kin}$, and the CO column density per velocity $N_{CO}/dv$ of the source for the cold and warm components.
<figure><img src="image_26.png" /><figcaption>**Figure 23.** The posterior probability distributions of the molecular gas density $n_H$, gas temperature $T_\text{kin}$, and CO column density per velocity $N_{\text{CO}}/\text{dv}$ of the source for the cold (left) and warm components (right), with the maximum posterior possibility point in the parameter space shown in orange lines and points, while dotted lines mark the ±1σ range in the distributions. The contours are in steps of 0.5σ starting from the center. The resulting parameter values are listed above the corresponding histograms.</figcaption></figure>

## ORCID iDs

P. Cox https://orcid.org/0000-0003-2027-8221

K. M. Butler https://orcid.org/0000-0001-7387-0558

C. R. Keeton https://orcid.org/0000-0001-6812-2467

L. Eid https://orcid.org/0009-0003-5696-9355

E. Borsato https://orcid.org/0009-0009-9483-8763

T. J. L. C. Bakx https://orcid.org/0000-0002-5268-2221

R. Neri https://orcid.org/0000-0002-7176-4046

B. M. Jones https://orcid.org/0000-0002-0675-0078

P. Prajapati https://orcid.org/0000-0002-3094-1077

A. J. Baker https://orcid.org/0000-0002-7892-396X

S. Berta https://orcid.org/0000-0002-0320-1532

A. Cooray https://orcid.org/0000-0002-3892-0190

E. M. Corsini https://orcid.org/0000-0003-3460-5633

L. Marchetti https://orcid.org/0000-0003-3948-7621

A. Omont https://orcid.org/0000-0002-4721-3922

A. Beelen https://orcid.org/0000-0003-3201-0185

R. Gavazzi https://orcid.org/0000-0002-5540-6935

D. Ismail https://orcid.org/0009-0007-2281-4944

R. J. Ivison https://orcid.org/0000-0001-5118-1313

M. Krips https://orcid.org/0000-0001-6971-4851

M. D. Lehnert https://orcid.org/0000-0003-1939-5885

H. Messias https://orcid.org/0000-0002-2985-7994

D. Riechers https://orcid.org/0000-0001-9585-1462

C. Vlahakis https://orcid.org/0000-0003-3745-4228

A. Weiss https://orcid.org/0000-0003-4678-3939

P. van der Werf https://orcid.org/0000-0001-5434-5942

C. Yang https://orcid.org/0000-0002-8117-9991

## References

Adam, G., Bacon, R., Courtes, G., et al. 1989, A&A, 208, L15

Ahumada, R., Allende Prieto, C., Almeida, A., et al. 2020, ApJS, 249, 3

ALMA Partnership, Vlahakis, C., Hunter, T. R., et al. 2015, ApJL, 808, L4

Bartelmann, M. 1996, A&A, 313, 697

Barvainis, R., Tacconi, L., Antonucci, R., Alloin, D., & Coleman, P. 1994, Natur, 371, 586

Berta, S., Young, A. J., Cox, P., et al. 2021, A&A, 646, A122

Berta, S., Stanley, F., Ismail, D., et al. 2023, A&A, 678, A28

Bettoni, D., Falomo, R., Scarpa, R., et al. 2019, ApJL, 873, L14

Bolatto, A. D., Warren, S. R., Leroy, A. K., et al. 2013, Natur, 499, 450

Bolton, A. S., Moustakas, L. A., Stern, D., et al. 2006, ApJL, 646, L45

Borsato, E., Marchetti, L., Negrello, M., et al. 2024, MNRAS, 528, 6222

Bosch, J., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S5

Bradford, C. M., Aguirre, J. E., Aikin, R., et al. 2009, ApJ, 705, 112

Butler, K. M., van der Werf, P. P., Rybak, M., et al. 2021, ApJ, 919, 5

Calzetti, D. 2001, PASP, 113, 1449

Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105

Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé, R. 2018, MNRAS, 480, 4379

Carnall, A. C., McLure, R. J., Dunlop, J. S., et al. 2022, ApJ, 929, 131

CASA Team, Bean, B., Bhatnagar, S., et al. 2022, PASP, 134, 114501

Chabrier, G. 2003, PASP, 115, 763

Chantry, V., & Magain, P. 2007, A&A, 470, 467

Chartas, G., Eracleous, M., Agol, E., & Gallagher, S. C. 2004, ApJ, 606, 78

Chen, G. C. F., Fassnacht, C. D., Suyu, S. H., et al. 2019, MNRAS, 490, 1743

Chwolson, O. 1924, AN, 221, 329

Cikota, A., Bertolla, I. T., Huang, X., et al. 2023, ApJL, 953, L5

Conley, A., Cooray, A., Vieira, J. D., et al. 2011, ApJL, 732, L35

Cortés, P., Vlahakis, C., Hales, A., et al. 2024, ALMA Technical Handbook, ALMA Doc. 11.3, v1.4, Zenodo, doi: 10.5281/zenodo.4511521

Cox, P., Neri, R., Berta, S., et al. 2023, A&A, 678, A26

da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13

Delhaize, J., Smolčić, V., Delvecchio, I., et al. 2017, A&A, 602, A4

Draine, B. T., Aniano, G., Krause, O., et al. 2014, ApJ, 780, 172

Dunne, L., Maddox, S. J., Papadopoulos, P. P., Ivison, R. J., & Gomez, H. L. 2022, MNRAS, 517, 962

Dunne, L., Maddox, S. J., Vlahakis, C., & Gomez, H. L. 2021, MNRAS, 501, 2573

Eales, S., Dunne, L., Clements, D., et al. 2010, PASP, 122, 499

Einasto, M., Einasto, J., Tenjes, P., et al. 2024, A&A, 681, A91

Einstein, A. 1936, Sci, 84, 506

Fox, D. C., & Pen, U.-L. 2001, ApJ, 546, 35

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1

Gilman, D., Birrer, S., Treu, T., Nierenberg, A., & Benson, A. 2019, MNRAS, 487, 5721
Girelli, G., Pozzetti, L., Bolzonella, M., et al. 2020, A&A, 634, A135

Hege, E. K., Hubbard, E. N., Strittmatter, P. A., & Worden, S. P. 1981, ApJL, 248, L1

Herrera-Camus, R., Förster Schreiber, N., Genzel, R., et al. 2021, A&A, 649, A31

Hezaveh, Y. D., Dalal, N., Marrone, D. P., et al. 2016, ApJ, 823, 37

Huchra, J., Gorenstein, M., Kent, S., et al. 1985, AJ, 90, 691

Inada, N., Oguri, M., Pindor, B., et al. 2003, Natur, 426, 810

Inoue, K. T., Minezaki, T., Matsushita, S., & Nakanishi, K. 2023, ApJ, 954, 197

Ismail, D., Beelen, A., Buat, V., et al. 2023, A&A, 678, A27

Jiang, L., Fan, X., Bian, F., et al. 2014, ApJS, 213, 12

Keeton, C. R. 2001, arXiv:astro-ph/0102340

Kelly, P. L., Rodney, S., Treu, T., et al. 2023, Sci, 380, abh1322

Kennicutt, R. C. J. 1989, ApJ, 344, 685

Kneib, J. P., Alloin, D., & Pello, R. 1998, A&A, 339, L65

Kochanek, C. S., & Schechter, P. L. 2004, in Measuring and Modeling the Universe, Carnegie Observatories Centennial Symp., ed. W. L. Freedman (Cambridge: Cambridge Univ. Press), 117

Lucey, J. R., Schechter, P. L., Smith, R. J., & Anguita, T. 2018, MNRAS, 476, 927

Magnelli, B., Ivison, R. J., Lutz, D., et al. 2015, A&A, 573, A45

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127

Meijerink, R., Kristensen, L. E., Weiß, A., et al. 2013, ApJL, 762, L16

Napolitano, N. R., Li, R., Spiniello, C., et al. 2020, ApJL, 904, L31

Natarajan, P., Williams, L. L. R., Bradač, M., et al. 2024, SSRv, 220, 19

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Nayyeri, H., Keele, M., Cooray, A., et al. 2016, ApJ, 823, 17

Negrello, M., Hopwood, R., De Zotti, G., et al. 2010, Sci, 330, 800

Neri, R., Cox, P., Omont, A., et al. 2020, A&A, 635, A7

Nierenberg, A. M., Keeley, R. E., Sluse, D., et al. 2024, MNRAS, 530, 2960

Nightingale, J., Amvrosiadis, A., Hayes, R., et al. 2023, JOSS, 8, 4475

Nightingale, J. W., Dye, S., & Massey, R. J. 2018, MNRAS, 478, 4738

Oliver, S. J., Bock, J., Altieri, B., et al. 2012, MNRAS, 424, 1614

Omont, A., Yang, C., Cox, P., et al. 2013, A&A, 551, A115

Pérez-Beaupuits, J. P., Güsten, R., Harris, A., et al. 2018, ApJ, 860, 23

Pettini, M., Christensen, L., D'Odorico, S., et al. 2010, MNRAS, 402, 2335

Planck Collaboration, Aghanim, N., Altieri, B., et al. 2015, A&A, 582, 29

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6

Ratnatunga, K. U., Ostrander, E. J., Griffiths, R. E., & Im, M. 1995, ApJL, 453, L5

Refsdal, S. 1964, MNRAS, 128, 307

Riechers, D. A. 2025, ApJ, 980, 59

Riechers, D. A., Walter, F., Brewer, B. J., et al. 2008, ApJ, 686, 851

Riechers, D. A., Carilli, C. L., Maddalena, R. J., et al. 2011, ApJL, 739, L32

Riechers, D. A., Bradford, C. M., Clements, D. L., et al. 2013, Natur, 496, 329

Schmidt, T., Treu, T., Birrer, S., et al. 2023, MNRAS, 518, 1260

Sengül, A. C., & Dvorkin, C. 2022, MNRAS, 516, 336

Sharon, C. E., Tagore, A. S., Baker, A. J., et al. 2019, ApJ, 879, 52

Sheu, W., Cikota, A., Huang, X., et al. 2024, ApJ, 973, 3

Spilker, J. S., Phadke, K. A., Aravena, M., et al. 2020a, ApJ, 905, 85

Spilker, J. S., Aravena, M., Phadke, K. A., et al. 2020b, ApJ, 905, 86

Stanley, F., Jones, B. M., Riechers, D. A., et al. 2023, ApJ, 945, 24

Stern, D., Djorgovski, S. G., Krone-Martins, A., et al. 2021, ApJ, 921, 42

Stoughton, D., Adelman, S. G., Annis, A., et al. 2002, Proc. SPIE, 4836, 339

Tagore, A. S., & Keeton, C. R. 2014, MNRAS, 445, 694

Tagore, A. S., & Jackson, N. 2016, MNRAS, 457, 3066

van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627

van der Werf, P. P., Isaak, K. G., Meijerink, R., et al. 2010, A&A, 518, L42

Vieira, J. D., Crawford, T. M., Switzer, E. R., et al. 2010, ApJ, 719, 763

Walsh, D., Carswell, R. F., & Weymann, R. J. 1979, Natur, 279, 381

Weiß, A., Henkel, C., Downes, D., & Walter, F. 2003, A&A, 409, L41

Weymann, R. J., Latham, D., Angel, J. R. P., et al. 1980, Natur, 285, 641

Wisotzki, L., Schechter, P. L., Bradt, H. V., Heinmüller, J., & Reimers, D. 2002, A&A, 395, 17

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

Wong, K. C., Suyu, S. H., Chen, G. C. F., et al. 2020, MNRAS, 498, 1420

Yang, C., Omont, A., Beelen, A., et al. 2017, A&A, 608, A144

Yang, C., Gavazzi, R., Beelen, A., et al. 2019, A&A, 624, A138

Yee, H. K. C. 1988, AJ, 95, 1331

Zhang, L., Zhang, Z.-Y., Nightingale, J. W., et al. 2023, MNRAS, 524, 3671

Zwicky, F. 1937, PhRv, 51, 679