# Close-up view of a luminous star-forming galaxy at z = 2.95

S. Berta, A. J. Young, P. Cox, R. Neri, B. M. Jones, A. J. Baker, A. Omont, L. Dunne, A. Carnero Rosell, L. Marchetti, et al.

► To cite this version:

S. Berta, A. J. Young, P. Cox, R. Neri, B. M. Jones, et al. Close-up view of a luminous star-forming galaxy at z = 2.95. Astronomy and Astrophysics - A&A, 2021, 646, pp.A122. 10.1051/0004-6361/202039743. hal-03147428

HAL Id: hal-03147428

https://hal.science/hal-03147428

Submitted on 19 Feb 2021

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.
Close-up view of a luminous star-forming galaxy at z = 2.95*

S. Berta¹, A. J. Young², P. Cox³, R. Neri¹, B. M. Jones⁴, A. J. Baker², A. Omont³, L. Dunne⁵, A. Carnero Rosell⁶,⁷,
L. Marchetti⁸,⁹,¹⁰, M. Negrello⁵, C. Yang¹¹, D. A. Riechers¹²,¹³, H. Dannerbauer⁶,⁷, I. Perez-Fournon⁶,⁷,
P. van der Werf¹⁴, T. Bakx¹⁵,¹⁶, R. J. Ivison¹⁷, A. Beelen¹⁸, V. Buat¹⁸, A. Cooray¹⁹, I. Cortzen¹, S. Dye²⁰, S. Eales⁵,
R. Gavazzi³, A. I. Harris²¹, C. N. Herrera¹, D. Hughes²², S. Jin⁶,⁷, M. Krips¹, G. Lagache¹⁸, M. Lehnert³,
H. Messias²³, S. Serjeant²⁴, F. Stanley³, S. Urquhart²⁴, C. Vlahakis²⁵, and A. Weiß²⁶

(Affiliations can be found after the references)

Received 22 October 2020 / Accepted 1 December 2020

ABSTRACT

Exploiting the sensitivity of the IRAM NOrthern Extended Millimeter Array (NOEMA) and its ability to process large instantaneous bandwidths, we have studied the morphology and other properties of the molecular gas and dust in the star forming galaxy, H-ATLAS J131611.5+281219 (HerBS-89a), at z = 2.95. High angular resolution (0.′3) images reveal a partial 1.′0 diameter Einstein ring in the dust continuum emission and the molecular emission lines of ¹²CO(9–8) and H₂O(2₀−1₁). Together with lower angular resolution (0.′6) images, we report the detection of a series of molecular lines including the three fundamental transitions of the molecular ion OH⁺, namely (1₀−0₁), (1₂−0₁), and (1₀−0₁), seen in absorption; the molecular ion CH⁺(1−0) seen in absorption, and tentatively in emission; two transitions of amidogen (NH₂), namely (2₀₂−1₁₁) and (2₂₀−2₁₁) seen in emission; and HCN(11−10) and/or NH(1₂−0₁) seen in absorption. The NOEMA data are complemented with Very Large Array data tracing the ¹²CO(1−0) emission line, which provides a measurement of the total mass of molecular gas and an anchor for a CO excitation analysis. In addition, we present Hubble Space Telescope imaging that reveals the foreground lensing galaxy in the near-infrared (1.15 µm). Together with photometric data from the Gran Telescopio Canarias, we derive a photometric redshift of zphot = 0.9+0.3-0.5 for the foreground lensing galaxy. Modeling the lensing of HerBS-89a, we reconstruct the dust continuum (magnified by a factor μ ≈ 5.0) and molecular emission lines (magnified by μ ~ 4–5) in the source plane, which probe scales of ~0.′1 (or 800 pc). The ¹²CO(9–8) and H₂O(2₀₂−1₁₁) emission lines have comparable spatial and kinematic distributions; the source-plane reconstructions do not clearly distinguish between a one-component and a two-component scenario, but the latter, which reveals two compact rotating components with sizes of ≈1 kpc that are likely merging, more naturally accounts for the broad line widths observed in HerBS-89a. In the core of HerBS-89a, very dense gas with nH₂ ~ 10⁷⁻⁹ cm⁻³ is revealed by the NH₂ emission lines and the possible HCN(11−10) absorption line. HerBS-89a is a powerful star forming galaxy with a molecular gas mass of Mmol = (2.1 ± 0.4) × 10¹¹ M☉, an infrared luminosity of LIR = (4.6 ± 0.4) × 10¹² L☉, and a dust mass of MDust = (2.6 ± 0.2) × 10⁹ M☉, yielding a dust-to-gas ratio δGDR ≈ 80. We derive a star formation rate SFR = 614 ± 59 M☉ yr⁻¹ and a depletion timescale τdepl = (3.4 ± 1.0) × 10⁸ years. The OH⁺ and CH⁺ absorption lines, which trace low (~100 cm⁻³) density molecular gas, all have their main velocity component red-shifted by ΔV ~ 100 km s⁻¹ relative to the global CO reservoir. We argue that these absorption lines trace a rare example of gas inflow toward the center of a galaxy, indicating that HerBS-89a is accreting gas from its surroundings.

**Key words.** galaxies: high-redshift – galaxies: ISM – galaxies: star formation – gravitational lensing: strong – submillimeter: galaxies – radio lines: ISM

**1. Introduction**

In the last two decades, surveys in the far-infrared and submillimeter wavebands have opened up a new window for our understanding of the formation and evolution of galaxies, revealing a population of massive, dust-enshrouded galaxies forming stars at enormous rates in the early Universe (see, e.g., Blain et al. 2002; Carilli & Walter 2013; Casey et al. 2014; Hodge & da Cunha 2020). In particular, the extragalactic imaging surveys done with the Herschel Space Observatory (Pilbratt et al. 2010), such as Herschel-ATLAS (Eales et al. 2010a), HerMES (Oliver et al. 2012), and PEP (Lutz et al. 2011), have increased the number of dust-obscured star-forming galaxies from hundreds to several hundred thousand. Together with other large-area surveys, such as the all-sky Planck-HFI (Planck Collaboration Int. XXVII 2015), the South Pole Telescope (SPT; Carlstrom et al. 2011) cosmological survey (Staniszewski et al. 2009; Vieira et al. 2010), and the Atacama Cosmology Telescope (ACT; Marsden et al.

2014; Gralla et al. 2020), today we have vast samples of luminous dusty star-forming galaxies (DSFGs) that are amongst the brightest galaxies in the Universe, including numerous examples of strongly lensed systems (e.g., Negrello et al. 2010, 2017; Cox et al. 2011; Bussmann et al. 2013; Spilker et al. 2014; Cañamas et al. 2015; Reuter et al. 2020) and rare cases of galaxies with intrinsic infrared luminosities, L<sub>FIR</sub> ≳ 10<sup>13</sup> L<sub>☉</sub>, and star formation rates (SFRs) in excess of 1000 M<sub>☉</sub> yr<sup>-1</sup>, known as Hyper-Luminous Infrared Galaxies (HyLIRGs, see, e.g., Ivison et al. 1998, 2013, 2019; Fu et al. 2013; Oteo et al. 2016; Riechers et al. 2013, 2017).

Detailed follow-up studies of these galaxies to investigate their nature and physical properties require robust estimates of their distances. Due to the dust obscuration in these objects, searching for CO emission lines via millimeter and submillimeter spectroscopy has proved to be the most reliable method for measuring accurate redshifts. This approach has become more and more efficient thanks to the increased bandwidths of the receivers and backends, most notably at the NOrthern Extended Millimeter Array (NOEMA) and the Atacama Large Millimeter/submillimeter Array (ALMA; e.g., Weiß et al. 2013;

<footnote>* The reduced datacubes are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/646/A122</footnote>
Fudamoto et al. 2017; Neri et al. 2020; Reuter et al. 2020, and references therein).

Using NOEMA, Neri et al. (2020) reported the results of a project whose aim was to measure robust spectroscopic redshifts for 13 bright Herschel-selected galaxies with $S_{500\mu m} > 80$ mJy, preferentially selecting lensed systems (Negrello et al. 2010). Reliable spectroscopic redshifts were derived for 12 individual sources, demonstrating the efficiency of deriving redshifts of high-z galaxies using the new correlator and broadband receivers on NOEMA. Based on the success of this Pilot Program, we started a comprehensive redshift survey of a sample of 125 of the brightest Herschel-selected galaxies (using the same selection criteria as above), the NOEMA Large Program z-GAL, whose main scientific goal is to further characterize the properties of luminous DSFGs in the early Universe. Interestingly, half of the sources in the Pilot Program sample display CO emission line widths in excess of 800 km s⁻<fnref n="1" />. Based on their estimated locations relative to the $L'_{CO(1-0)}$ versus ΔV(CO) relationship of Harris et al. (2012), several of the sources are inferred to be gravitationally amplified, while a number of them appear to belong to the rare class of hyper-luminous infrared galaxies (Neri et al. 2020).

One of these galaxies, H-ATLAS J131611.5+281219 (hereafter HerBS-89a), at z = 2.9497, displays a very strong 2 mm continuum (with a flux density $S_{159\text{GHz}} = 4.56 \pm 0.05$ mJy) and CO emission lines that are the broadest of the entire sample with a line width (FWHM) of ΔV ~ 1100 km s⁻<fnref n="1" />. Both the 2 mm continuum and ¹²CO(5–4) line emission were resolved by the 1.′2 imaging, with an extension of 0.′9 ± 0.′1 in the east-west direction. The corresponding infrared luminosity of HerBS-89a (between 8 and 1000 µm in the rest-frame; uncorrected for amplification) is estimated to be (2.9 ± 0.2) × 10¹³ $L_\odot$, which suggested that HerBS-89a could be a HyLIRG.

In order to further explore the properties of HerBS-89a, we used NOEMA to perform high angular resolution observations at 1 mm. We observed, in addition to the underlying dust continuum, the emission lines of ¹²CO(9–8) and H₂O(2₀₂–1₁₁); the three lines of the ground state of the molecular ion OH⁺, (1₁–0₁), (1₂–0₁) and (1₀–0₁), which are for the first time reported together in a high-z galaxy; the CH⁺(1–0) line seen in absorption (and tentatively in emission); NH₂(2₀₂–1₁₁)(5/2–3/2) seen in emission; and HCN(11–10) and/or NH(1₂–0₁) seen in absorption, here also reported for the first time in a high-z galaxy.

In HerBS-89a, the 0.′3 images of the molecular emission lines and the dust continuum reveal a partial 1.′0 diameter Einstein ring. All of the observed transitions of the molecular ions OH⁺ and CH⁺ are seen in absorption against the strong dust continuum. Together, the molecular emission and absorption lines allow us to probe simultaneously in HerBS-89a regions that have very different properties, ranging from the (very) dense molecular gas traced by ¹²CO(9–8), H₂O(2₀₂–1₁₁), NH₂ and HCN(11–10) (overlapping with NH(1₂–0₁)), to the low H₂ fraction, diffuse gas traced in OH⁺ and the reservoirs of turbulent, cold and low-density molecular gas traced by CH⁺, thereby providing a unique view at sub-kpc spatial resolution (in the source plane) of the physical properties and feedback activity in this high-z system.

The NOEMA data are complemented by Hubble Space Telescope (HST) imaging that traces the foreground massive lensing galaxy in the near-infrared (using the F110W filter around 1.1 µm), and by optical/near-infrared data obtained with the Gran Telescopio Canarias (GTC) that provide constraints on the redshift of the lensing galaxy. In addition, we also present data obtained with the Karl G. Jansky Very Large Array (VLA) on

the ¹²CO(1–0) emission line, which allow us to derive the mass of molecular gas in HerBS-89a and give an anchoring point for an analysis of its CO spectral line energy distribution.

The structure of the paper is as follows. Section 2 describes the NOEMA, VLA, HST, and GTC observations, and the reduction of the respective data sets; Sect. 3 presents the main results, including the morphology of the source, the properties of the molecular emission and absorption lines and the underlying continuum; Sect. 4 describes the characteristics of the foreground lensing galaxy; Sect. 5 presents the lensing model and the morphology of HerBS-89a in the source plane; and Sect. 6 outlines the global intrinsic properties of the dust and molecular gas (derived from the CO and water emission lines) corrected for amplification, including the CO excitation. The molecular gas kinematics and the dynamical mass are reported in Sect. 7; the properties of the molecular absorption and emission lines other than CO and H₂O are outlined in Sect. 8; and a discussion of the gas inflow suggested by the red-shifted molecular absorption lines of OH⁺ and CH⁺ is presented in Sect. 9. Finally, Sect. 10 summarizes the main conclusions of this paper and outlines future prospects.

Throughout this paper we adopt a spatially flat $\Lambda$CDM cosmology with $H_0 = 67.4$ km s⁻<fnref n="1" /> Mpc⁻<fnref n="1" /> and $\Omega_M = 0.315$ (Planck Collaboration VI 2020) and assume a Chabrier (2003) initial mass function (IMF). At the redshift z = 2.9497 of HerBS-89a, one arc-second corresponds to 7.9 kpc and the luminosity distance to the source is $D_L = 2.5 \times 10^4$ Mpc.

## 2. Observations and data reduction

### 2.1. NOEMA

We used NOEMA to target high-frequency molecular lines in HerBS-89a, redshifted into the 1 mm band. The observations were carried out under two separate projects.

The first project was a Discretionary Directorial Time (DDT) project, labeled E18AE (P.I.: R. Neri), observed on February 5, 2019 with ten antennas using the extended A-configuration, yielding an angular resolution of 0.′3, for a total observing time of 4.3 h on-source time. This project, which was specifically tailored to measure the ¹²CO(9–8) and H₂O(2₀₂–1₁₁) emission lines, also enabled the detection of two strong absorption lines of the molecular ion OH⁺ corresponding to the red-shifted ground state transitions of OH⁺(1₁–0₁) and OH⁺(1₂–0₁) ($v_{\text{rest}} = 1033.118$ GHz and 971.803 GHz, respectively)<fnref n="1" />.

The second project, labeled W19DE (P.I.: S. Berta), was completed on March 30th, 2020 with ten antennas using the intermediate C-configuration, yielding an angular resolution of 1.′0 × 0.′6, for a total on-source observing time of 4.2 h. This project was a follow-up of the DDT to measure the third ground state transition of OH⁺(1₀–0₁) ($v_{\text{rest}} = 909.158$ GHz)<fnref n="1" /> and the ground state transition of the molecular ion CH⁺(1–0) ($v_{\text{rest}} = 835.137$ GHz)<fnref n="2" />.

Observing conditions were excellent for both projects with an atmospheric phase stability of 20° rms and 0.8 mm of precipitable water vapor. The correlator was operated in the low resolution mode to provide spectral channels with nominal resolution of 2 MHz. The NOEMA antennas were equipped with 2SB receivers that cover a spectral window of 7.744 GHz in each

<footnote n="1">¹ The adopted frequencies correspond to the strongest of the hyperfine transitions for the specific $\Delta J$ (Indriolo et al. 2015, and references therein).</footnote>

<footnote n="2">² The rest frequency of the CH⁺(1–0) transition is taken from Müller (2010).</footnote>
sideband and polarization. For the first series of observations, we covered the frequency range from 244.9 to 252.6 GHz and 260.4 to 268.1 GHz; for the second series, the frequency ranges were from 209.4 to 217.1 GHz and 224.9 to 232.6 GHz.

The strength of the continuum in HerBS-89a at 1 mm (~20 mJy – see Table 1) enabled phase self-calibration, which significantly improved the image fidelity and dynamic range both in the continuum and molecular line emission (see Sects. 3.1.1 and 3.1.2). The flux calibrator used in both projects was Lk Hα101. The phase calibrator for the high-angular resolution project (E18AE) was 1308+326. The data were calibrated, averaged in polarization, mapped, and analyzed in the GILDAS software package. The absolute flux calibration is accurate to within 10% and the $1\sigma$ error on the absolute position of HerBS-89a is estimated to be <0.2.

## 2.2. VLA

The National Radio Astronomy Observatory's (NRAO) VLA was used to observe the $^{12}$CO(1–0) emission line in HerBS-89a when the array was in the C configuration. The observations are part of a larger project (program I.D.:VLA/20A-083 – P.I.: D. Riechers) whose goal was to measure the $^{12}$CO(1–0) emission lines in the sample of Herschel-selected galaxies studied in Neri et al. (2020); the complete results of that project will be presented in a forthcoming paper (Stanley et al., in prep.). The data were acquired on March 22, 2020 under stable atmospheric conditions. The 0.9 cm Ka band receivers were tuned (1 GHz bandwidth per IF) to the expected frequency of the CO emission line, that is 29.1848 GHz based on the redshift determined by Neri et al. (2020), and to 38.499 GHz for the second baseband. In total, we observed for 1.6 h (with 62 min on-source) recording 1024×2 MHz dual-polarization channels across a total bandwidth of 2 GHz, which was chosen to maximize line sensitivity while retaining as much bandwidth as possible. The 2 GHz bandwidth setup was used to maximize the potential for stacking of faint lines (for the entire VLA project), while at same time retaining sufficient spectral resolution (2 MHz, i.e., 16–21 km s⁻¹) to finely sample the broad CO emission line of HerBS-89a. The 8-bit samplers were selected to maximize sensitivity. The source 3C286 was observed to determine accurate complex gain and bandpass correction solutions and was also observed to set the absolute scale flux density based on the Perley & Butler (2017) models; the pointing accuracy was checked regularly. The data were calibrated, averaged in polarization, mapped using natural baseline weighting and analyzed in the CASA (Common Astronomy Software Applications) package. The resulting line map has a spatial resolution of 1.24 × 0.79 (PA –62.3°) and a rms noise level of 32.1 µJy beam⁻¹ over a band width of 0.13 GHz. The absolute flux scale is estimated to be accurate to within 10%.

## 2.3. HST

HerBS-89a was observed with the HST in March 2012 as part of a cycle-19 SNAPshot proposal (program I.D.: 12488 – P.I.: M. Negrello), which aimed, amongst other goals, at characterizing the nature of the putative lenses in a large sample of candidate lensing systems selected at 500 µm in the Herschel extragalactic surveys. Observations were obtained with the Wide Field Camera 3 (WFC3) using the wide-J filter F110W (peak wavelength 1.15 µm). The total exposure time is 252 s. Data were reduced using the Python ASTRODRIZZLE package, with parameters optimized to improve the final image

<table><caption>**Table 1.** Far-infrared and submillimeter continuum flux densities of HerBS-89a (and HerBS-89b).</caption><thead><tr><th rowspan="2">λ [mm]</th><th rowspan="2">ν [GHz]</th><th>HerBS-89a</th><th>HerBS-89b</th><th rowspan="2">Ref.</th></tr><tr><th>S<sub>ν</sub> [mJy]</th><th></th></tr></thead><tbody><tr><td></td><td></td><td colspan="2"><em>Herschel</em></td><td></td></tr><tr><td>0.10</td><td></td><td>&lt;130.8</td><td>—</td><td>(1,2)</td></tr><tr><td>0.16</td><td></td><td>&lt;136.2</td><td>—</td><td>"</td></tr><tr><td>0.25</td><td></td><td>71.8 ± 5.7</td><td>—</td><td>"</td></tr><tr><td>0.35</td><td></td><td>103.4 ± 5.7</td><td>—</td><td>"</td></tr><tr><td>0.50</td><td></td><td>95.7 ± 7.0</td><td>—</td><td>"</td></tr><tr><td></td><td></td><td colspan="2"><em>SCUBA-2</em></td><td></td></tr><tr><td>0.85</td><td></td><td>52.8 ± 4.3</td><td>—</td><td>(3)</td></tr><tr><td></td><td></td><td colspan="2"><em>NOEMA</em></td><td></td></tr><tr><td>1.14</td><td>263.3</td><td>22.33 ± 0.04</td><td>0.66 ± 0.04</td><td>(4)</td></tr><tr><td>1.21</td><td>248.3</td><td>19.22 ± 0.04</td><td>0.55 ± 0.04</td><td>"</td></tr><tr><td>1.31</td><td>229.3</td><td>13.62 ± 0.05</td><td>0.34 ± 0.05</td><td>"</td></tr><tr><td>1.41</td><td>212.4</td><td>11.31 ± 0.04</td><td>0.25 ± 0.04</td><td>"</td></tr><tr><td>1.89</td><td>158.6</td><td>4.56 ± 0.05</td><td>0.24 ± 0.05</td><td>(5)</td></tr><tr><td>2.01</td><td>149.0</td><td>3.40 ± 0.30</td><td>&lt;0.1</td><td>"</td></tr><tr><td>2.09</td><td>143.2</td><td>3.02 ± 0.04</td><td>—</td><td>"</td></tr><tr><td>2.25</td><td>133.5</td><td>2.20 ± 0.20</td><td>—</td><td>"</td></tr><tr><td>2.69</td><td>111.5</td><td>1.10 ± 0.08</td><td>—</td><td>(5)</td></tr><tr><td>2.89</td><td>103.7</td><td>0.83 ± 0.05</td><td>—</td><td>"</td></tr><tr><td>3.12</td><td>96.0</td><td>0.56 ± 0.06</td><td>—</td><td>"</td></tr><tr><td>3.40</td><td>88.3</td><td>0.44 ± 0.06</td><td>—</td><td>"</td></tr><tr><td></td><td></td><td colspan="2"><em>VLA</em></td><td></td></tr><tr><td>7.89</td><td>38.0</td><td>&lt;0.147</td><td>—</td><td>(4)</td></tr><tr><td>10.34</td><td>29.0</td><td>&lt;0.065</td><td>—</td><td>"</td></tr></tbody></table>

**Notes.** The widths of each of the NOEMA sidebands is 7.744 GHz and their frequency ranges are provided for the 1 mm data in Sect. 2.1 and in Table 2 of Neri et al. (2020) for the 2 and 3 mm data. The flux densities of HerBS-89b have been corrected for the primary beam. The quoted statistical uncertainties do not account for the absolute flux calibration uncertainty. All the upper limits are $3\sigma$.

**References.** (1) Bakx et al. (2018); (2) Valiante et al. (2016); (3) Bakx et al. (2020); (4) this paper; (5) Neri et al. (2020).

quality. The pixel scale of the Infrared-Camera is 0.128, but the image was resampled to a finer pixel scale of 0.064 by exploiting the adopted sub-pixel dither pattern. The astrometry of the resulting image was calibrated by matching the positions of 11 stars from SDSS DR12 and 2MASS and is accurate to within 0.1.

## 2.4. GTC

HerBS-89a was observed with the GTC 10.4 m telescope in two observing runs (program ID: GTC03-19ADDT and GTC09-19ADDT; PI: H. Dannerbauer). First, optical imaging was obtained using the instrument OSIRIS on April 6, 2019. The observations were conducted in service mode under clear skies (although with non-photometric conditions), integrating for 10 min in the Sloan r-band filter with a seeing of 0.8. The field of view is 7.5 × 6.0 and the pixel size 0.254. Standard procedures for data reduction and calibration of the raw images were performed with the IRAF data reduction package (Tody 1986).

The second series of observations was performed on June 7, 2020 using the visitor instrument HiPERCAM for a total observing time of 1.25 h under average seeing conditions (0.9). HiPERCAM is a quintuple-beam CCD imager enabling the five
<figure><img src="image_1.png" /></figure>

**Fig. 1.** Continuum images of HerBS-89a at 1.2 mm, resulting from merging the lower and upper side bands of the two sets of NOEMA observations obtained with natural weighting of the visibilities (see Sect. 3.1.1 for further details). Contours are at 5, 10, 20, 40, 80$\sigma$ levels. North is up, east is left. *Left*: high-angular resolution continuum image at 254.6 GHz with a beam size of 0.43 $\times$ 0.22. *Right*: lower-angular resolution image at 220.8 GHz with a beam size of 1.1 $\times$ 0.66. The synthesized beam is shown in the lower left corner of each image, and the cross indicates the position of the phase center (RA 13:16:11.52, Dec +28:12:17.7). The weak source HerBS-89b, about 6.0 to the east, is detected at both frequencies; the position of the source is indicated with an arrow and further details are provided in Sect. 3.1.1.

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

**Fig. 2.** Zoom-in of the continuum image of HerBS-89a at 254.6 GHz, obtained with uniform weighting, i.e., enhanced resolution, showing the details of the dust-continuum emission, in particular the faint extended emission of the ~1.0 diameter partial Einstein ring. See caption of Fig. 1 for further details.

Sloan filter *ugriz* to be observed simultaneously over a field of view of 2.8 $\times$ 1.4 with a pixel size of 0.081 (Dhillon et al. 2018, and references therein). The HiPERCAM team designed a dedicated data reduction tool<fnref n="3" /> that applies standard procedures for the reduction and calibration of the raw images.

## 3. NOEMA and VLA results

In this section, we describe the new data obtained on HerBS-89a, outline the properties derived from the NOEMA 1 mm high-angular resolution observations for both the continuum and the molecular emission and absorption lines (Sect. 3.1), and present the results on the $^{12}$CO(1–0) emission line measured with the VLA (Sect. 3.2).

### 3.1. NOEMA results

#### 3.1.1. Continuum emission

Figure 1 presents the two continuum maps of HerBS-89a obtained by merging the lower and upper side-bands of each set of observations, resulting in a high-angular resolution image centered at ~254.6 GHz and a lower angular resolution image at ~220.8 GHz. We reach a sensitivity of $\sigma = 29$ and $33~\mu$Jy beam$^{-1}$ in the two bands, respectively. Merging the two side-bands of the high-resolution continuum images improves the S/Ns and resulted, by applying a uniform weight, in a final image with beam size of 0.43 $\times$ 0.22, to be compared to the 0.55 $\times$ 0.33 achieved by the A-configuration that was used for these observations with natural weighting.

Figure 2 presents a zoom-in on this high-resolution continuum image of HerBS-89a. It displays an Einstein ring-like morphology, with a double source linked by weak arc structures, indicating that HerBS-89a is lensed. The high quality of the image reveals the details of the morphology of the dust continuum, such as the differences between the northern and southern continuum peaks as well as the weaker emission extending between them.

For the lower frequency data obtained using the medium-compact C-configuration of NOEMA, the source is only marginally resolved. Combining the data from the upper and lower side-bands resulted in a beam of 1.1 $\times$ 0.66.

<footnote n="3">³ http://deneb.astro.warwick.ac.uk/phsaap/hipercam/docs/html/</footnote>
To the east of HerBS-89a, the weak unresolved source (labeled HerBS-89b), which was already detected in 2 mm continuum by Neri et al. (2020), is detected in the 1.2 mm continuum emission in both images (Fig. 1), about 6" east of the phase center of our observations. The detection reported here confirms the authenticity of this source. However, there is no corresponding source in the SDSS catalog at that position. The faintness of HerBS-89b precludes the extraction of spectroscopic information from the available NOEMA data.

Continuum aperture flux densities have been measured for each side-band separately. Table 1 summarizes the values, both for HerBS-89a and the serendipitous source HerBS-89b, which is seen to be ≈40× fainter than HerBS-89a. The effective frequencies of the adopted continuum bands and the statistical uncertainties on the flux densities are also provided. The latter should be added in quadrature to the 10% absolute flux calibration uncertainty (see Sect. 2.1).

### 3.1.2. Molecular emission and absorption lines

The large bandwidth of the NOEMA receivers allows us to search in HerBS-89a for redshifted emission and absorption lines over a wide range in frequency. The complete spectra of HerBS-89a, integrated over the areas subtended by the lines (see below) and normalized by the continuum, are displayed in Fig. 3. Each of the panels covers the frequency range of one sideband (LSB or USB) of one of the two observational projects discussed above, namely between 245 and 268 GHz (two upper panels) and between 209 and 233 GHz (two lower panels). Two strong molecular emission lines are present, $^{12}$CO(9–8) and H$_2$O(2$_{02}$−1$_{11}$), which both display wide profiles similar to those seen in the $^{12}$CO(3–2) and $^{12}$CO(5–4) emission lines (Neri et al. 2020). In addition, the three lines of the ground state of the molecular ion OH$^+$ are all detected in absorption: the (1$_1$−0$_1$) line, which is adjacent in frequency to the $^{12}$CO(9–8) emission line, and the (1$_2$−0$_1$) and (1$_0$−0$_1$) lines, which are detected for the first time in a high-$z$ galaxy. Blue-shifted from the strong OH$^+(1_2-0_1)$ absorption line, another weaker absorption line is detected that is due to HCN(11–10) ($v_{\text{rest}}$ = 974.487 GHz) and/or NH(1$_2$−0$_1$) ($v_{\text{rest}}$ = 974.471 GHz) (see Sect. 8.1 for a detailed discussion). Next to the OH$^+(1_0-0_1)$ absorption line, another line is seen in emission, which corresponds to o−NH$_2$(2$_{02}$−1$_{11}$)(5/2−3/2) ($v_{\text{rest}}$ = 907.433 GHz). A second NH$_2$ line is detected with lower signal-to-noise ratio (S/N) at higher frequency, the NH$_2$(2$_{20}$−1$_{11}$)(5/2−3/2) transition at $v_{\text{rest}}$ = 993.3226 GHz. Finally, the molecular ion CH$^+$ is detected in the ground transition, CH$^+(1-0); its profile is dominated by an absorption line that has a width similar to those of the OH$^+$ absorption lines and at the same red-shifted velocity; in addition, a weak and broad (~850 km s$^{-1}$) emission component is likely present, of which the extreme red and blue wings are detected at low S/N at either side of the absorption line.

Also identified in Fig. 3 are the positions of red-shifted molecular lines that fall within the observed frequency range but remained undetected in HerBS-89a, including: SH$^+(2-1)$ ($v_{\text{rest}}$ = 893.152 GHz), o−NH$_2$(2$_{02}$−1$_{11}$)(3/2−1/2) ($v_{\text{rest}}$ = 902.210 GHz), H$_2$O$^+(3_{12}-3_{03})$ ($v_{\text{rest}}$ = 982.955 GHz), and o−H$_3$O$^+(0_0^-1_0^-)$ ($v_{\text{rest}}$ = 984.708 GHz), most of which are seen in the spectrum of the local merger Arp 220 (Rangwala et al. 2011).

Figure 4 shows the angular extents of the eight molecular emission and absorption lines detected with NOEMA in HerBS-89a, each compared to the continuum extracted in the respective sideband to which it belongs. The $^{12}$CO(9–8) and H$_2$O(2$_{02}$−1$_{11}$)

emission lines are detected in both lensed components of HerBS-89a, but are slightly shifted in position with respect to each other, with the water line being well centered on the continuum emission peaks. In order to recover more information on the distribution of the molecular emission lines, we combined the $^{12}$CO(9–8) and H$_2$O(2$_{02}$−1$_{11}$) velocity averaged images, under the assumption that they are probing roughly the same reservoir of high density and temperature molecular gas in HerBS-89a. The increased S/N of the combined image (see upper right panel in Fig. 4) reveals the molecular emission to the west, joining the southern and northern peaks, and resembling the image of the dust continuum emission.

The two strong OH$^+$ absorption lines detected in our higher angular resolution data, namely the (1$_1$−0$_1$) and (1$_2$−0$_1$) transitions, also coincide with the southern and northern peaks of the continuum emission. The combined OH$^+(1_1-0_1)$ and (1$_2$−0$_1$) map, integrated between −400 and +400 km s$^{-1}$ from the line center in order to avoid contamination by neighboring lines (right panel in the second row of Fig. 4), further highlights this spatial coincidence.

The faint absorption lines of our lower angular resolution data are detected mainly toward the southern component of HerBS-89a; the CH$^+$ absorption seems to be more extended than the OH$^+$ absorption. The image of the HCN(11–10) and/or NH(1$_2$−0$_1$) absorption line shows two peaks exactly centered on the dust continuum emission peaks with no indication of emission extending further out. Finally, the NH$_2$(2$_{02}$−1$_{11}$)(5/2−3/2) emission line peaks westward of the continuum peak, as expected for red-shifted gas (see below).

In Fig. 5, we show the integrated spectra of all the individual molecular emission and absorption lines detected in HerBS-89a, including the $^{12}$CO(3–2) and $^{12}$CO(5–4) emission lines reported in Neri et al. (2020). In the following, we discuss the properties of each of these molecular lines.

For the two low-J CO emission lines, Neri et al. (2020) fitted each of the double-peaked line profiles with two Gaussian components, with identical widths ($FWHM$ = 631 km s$^{-1}$) and central velocities shifted by ±289 km s$^{-1}$ from the nominal line frequency at the redshift of HerBS-89a ($z$ = 2.9497) (Table 2). We adopt here the same approach for the $^{12}$CO(9–8), H$_2$O(2$_{02}$−1$_{11}$) and NH$_2$(2$_{02}$−1$_{11}$)(5/2−3/2) emission lines, fixing their component velocities to the values found for $^{12}$CO(5–4) and (3–2). We also force the two Gaussian components (of each line) to have the same width (see Table 2). For the NH$_2$(2$_{20}$−1$_{11}$)(5/2−3/2) emission line, the detection has too low S/N to enable a reasonable fit and we therefore here report only upper limits to its line flux.

As mentioned above, the CH$^+(1-0)$ line seems to display weak emission features on both sides of the absorption line. The fit to the profile shows the overlapping deep absorption feature and the less-well defined, broad component seen in emission.

The OH$^+(1_2-0_1)$ absorption line displays strong, red-shifted absorption, whose profile is well reproduced by a narrow and deep absorption line, red-shifted by ~100 km s$^{-1}$. Next to this OH$^+$ absorption line, the broad weak absorption line, which corresponds to HCN(11–10) and/or NH(1$_2$−0$_1$), is centered at zero velocity and has a width of ~800 km s$^{-1}$. As shown in Fig. 4, this broad absorption line is very concentrated spatially, peaking on the northern and southern peaks of the dust continuum emission.

The OH$^+(1_0-0_1)$ and NH$_2$(2$_{02}$−1$_{11}$)(5/2−3/2) lines are close to each other in frequency (i.e., 230.18 and 229.74 GHz, corresponding to a velocity difference of $\Delta V$ ~ 150 km s$^{-1}$). In the observed spectrum (Fig. 5) there is an apparent shift for both species to the red by ~250 km s$^{-1}$ with respect to their nominal
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Spectra of HerBS-89a in the frequency ranges between 245 and 268 GHz and between 209 and 233 GHz. The spectra have been normalized by the continuum, which was modeled with a linear function (see text). The velocity channels were binned to 40 km s⁻¹ for the 245–268 GHz data and 60 km s⁻¹ for the 209–233 GHz data, reaching rms of 0.6 and 0.4 mJy per channel, respectively. The rest-frame and observed frequencies are given on the upper and lower horizontal axes, respectively. All the detected molecular emission and absorption lines are identified (with the arrows positioned at the line frequencies). In addition, the redshifted positions of molecular lines (at z = 2.9497) that fall within the observed frequency range, but were not detected, are indicated (in italics): o – H₃O⁺(0₀⁻¹₀), H₂O⁺(3₁₂⁻3₀₃), SH⁺(2⁻¹), and o – NH₂(2₀,₂⁻¹₁,₁)(3/2⁻¹/₂).</figcaption></figure>

frequencies. For the OH⁺(1₀–0₁) absorption, this shift is larger than for the two other OH⁺ transitions, and in the case of NH₂, the profile seems to be inconsistent with those seen for other molecular emission lines tracing dense gas (CO and H₂O). Since

NH₂ (amidogen) is a photodissociation product of NH₃ (ammo-
nia), it should probe the dense gas traced in the high-*J* CO or
water lines, and therefore display a similar double Gaussian pro-
file. However, we only see the red-shifted component of the
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Images of the molecular emission and absorption lines detected in HerBS-89a, compared to the underlying continuum (natural weighting; shown as contours with $10\sigma$ spacing) extracted in the respective side-band where each molecular line was detected (for each image the continuum has been subtracted from the line emission). The *three upper rows* display the high-angular resolution data, whereas the *bottom row* shows the lower-angular resolution images. **Upper right panel:** combined $^{12}$CO(9–8) and H$_2$O(2$_{02}-1_{11}$) image that enhances the low-level line emission, in particular along the western arc of the partial Einstein ring. **Right panel in the second row:** combined map of OH$^+$ ((1$_1-0_1$) and (1$_2-0_1$)) absorption lines. The synthesized beam is shown in the lower left corner of each of the images.</figcaption></figure>

double-peaked line profile, most likely because the blue-shifted part of the NH$_2$ emission line has been entirely suppressed by the OH$^+$ absorption line. We have therefore fitted simultaneously the overlapping OH$^+$ absorption and NH$_2$ emission lines using three Gaussian profiles, namely: (i) two components for NH$_2$, where the velocities are derived from the $^{12}$CO(5–4) and $^{12}$CO(3–2) emission lines and the ratio of the blue and red peaks similar to that of the double Gaussian profile of the $^{12}$CO(5–4) emis-

sion line; (ii) one Gaussian for the OH$^+$ absorption line whose width and the center velocity are derived from the narrow and deep OH$^+(1_2-0_1)$ absorption line. The fit shown in Fig. 5 reproduces very well the observed spectrum including the velocity shift described above.

Table 2 summarizes the results of the spectral line fitting. For each molecular line, the main properties of the best Gaussian fits are listed, including the central velocity and frequency, the full
<table><caption>**Table 2.** Summary of molecular emission and absorption line properties in HerBS-89a, obtained by the line fitting.</caption><thead><tr><th>Molecule</th><th>Transition</th><th>ν<sub>rest</sub><br>[GHz]</th><th>ν<sub>obs</sub><br>[GHz]</th><th>V<br>[km s<sup>-1</sup>]</th><th>FWHM<br>[km s<sup>-1</sup>]</th><th>I<sub>line</sub><br>[Jy km s<sup>-1</sup>]</th><th>μL'<sub>line</sub><br>[10<sup>10</sup> K km s<sup>-1</sup>pc<sup>-2</sup>]</th></tr></thead><tbody><tr><td colspan="8" style="text-align: center; font-weight: bold;">Emission lines</td></tr><tr><td rowspan="4"><sup>12</sup>CO</td><td>1 – 0</td><td>115.271</td><td>29.186</td><td>-57.6 ± 119.5</td><td>1443 ± 293</td><td>0.64 ± 0.13</td><td>25.76 ± 5.23</td></tr><tr><td>3 – 2</td><td>345.796</td><td>87.549</td><td>0.0</td><td>1102 ± 83</td><td>4.00 ± 0.60</td><td>17.75 ± 2.66</td></tr><tr><td>5 – 4</td><td>576.268</td><td>145.902</td><td>0.0</td><td>1028 ± 107</td><td>8.40 ± 0.80</td><td>13.42 ± 1.28</td></tr><tr><td>9 – 8</td><td>1036.912</td><td>262.530</td><td>0.0</td><td>1168 ± 94</td><td>1.78 ± 0.21</td><td>0.88 ± 0.10</td></tr><tr><td>H<sub>2</sub>O</td><td>2<sub>02</sub> – 1<sub>11</sub></td><td>987.927</td><td>250.127</td><td>0.0</td><td>1128 ± 174</td><td>1.59 ± 0.37</td><td>0.86 ± 0.20</td></tr><tr><td rowspan="2">NH<sub>2</sub></td><td>2<sub>02</sub> – 1<sub>11</sub> (<sup>†</sup>)</td><td>907.433</td><td>229.747</td><td>0.0</td><td>972 ± 237<sup>(*)</sup></td><td>0.64 ± 0.27<sup>(*)</sup></td><td>0.41 ± 0.17</td></tr><tr><td>2<sub>20</sub> – 2<sub>11</sub> (<sup>†</sup>)</td><td>993.322</td><td>251.493</td><td>0.0</td><td></td><td>&lt;0.33</td><td></td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Absorption lines</td></tr><tr><td rowspan="3">OH<sup>+</sup></td><td>1<sub>1</sub> – 0<sub>1</sub></td><td>1033.118</td><td>261.546</td><td>+102.2 ± 13.0</td><td>447 ± 32</td><td>-.96 ± 0.14</td><td></td></tr><tr><td>1<sub>2</sub> – 0<sub>1</sub></td><td>971.804</td><td>246.026</td><td>+91.6 ± 20.1</td><td>558 ± 51</td><td>-3.08 ± 0.28</td><td></td></tr><tr><td>1<sub>0</sub> – 0<sub>1</sub></td><td>909.159</td><td>230.132</td><td>+267.6 ± 32.2</td><td>349 ± 121</td><td>-0.28 ± 0.15<sup>(**)</sup></td><td></td></tr><tr><td>CH<sup>+</sup></td><td>1 – 0</td><td>835.137</td><td>211.419</td><td>+137.6 ± 12.6</td><td>354 ± 108<sup>(***)</sup></td><td>-0.39 ± 0.09<sup>(***)</sup></td><td></td></tr><tr><td>HCN</td><td>11 – 10<sup>(‡)</sup></td><td>974.487</td><td>246.724</td><td>0.0</td><td>523 ± 239</td><td>-0.63 ± 0.28</td><td></td></tr><tr><td>NH</td><td>1<sub>2</sub> – 0<sub>1</sub><sup>(‡)</sup></td><td>974.471</td><td>246.720</td><td>0.0</td><td>523 ± 239</td><td>-0.63 ± 0.28</td><td></td></tr></tbody></table>

**Notes.** The data for the <sup>12</sup>CO(3–2) and <sup>12</sup>CO(5–4) emission lines are from Neri et al. (2020). The fitting of each line is shown in Fig. 5. When one parameter is quoted without uncertainty, its value was fixed during the fitting procedure. The observed line luminosities μL'<sub>line</sub> have been derived using the standard relation given by Solomon & Vanden Bout (2005) – see Sect. 6.2 and Eq. (6). <sup>(*)</sup>Combination of the two Gaussian components. <sup>(**)</sup>Only the visible red component. <sup>(***)</sup>Net absorption only. <sup>(†)</sup>The two NH<sub>2</sub> transitions are (2<sub>02</sub>–1<sub>11</sub>) (5/2–3/2) and (2<sub>20</sub>–2<sub>11</sub>) (5/2–3/2), respectively. In the case of the latter transition, we provide a 3σ upper limit. <sup>(‡)</sup>In the case of the absorption line observed at 246.72 GHz, we here provide the derived properties for HCN and NH separately, noting that both molecules likely do contribute to the observed absorption (see Sect. 8.1).

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

**Fig. 5.** Spectra of the molecular emission and absorption lines detected in HerBS-89a reported in this study, together with the <sup>12</sup>CO(5–4) and <sup>12</sup>CO(3–2) spectra from Neri et al. (2020). Each molecular line is identified in the lower corner of each panel (labeled in black) and, in the case of a line overlap, as <sup>12</sup>CO(9–8)/OH<sup>+</sup>(1<sub>1</sub>–0<sub>1</sub>), OH<sup>+</sup>(1<sub>2</sub>–0<sub>1</sub>)/HCN(11–10) and OH<sup>+</sup>(1<sub>0</sub>–0<sub>1</sub>)/NH<sub>2</sub>(2<sub>02</sub>–1<sub>11</sub>)(5/2–3/2), the second molecular line is labeled in red. The spectra are displayed with the continuum subtracted and, in each panel, the molecular line (labeled in black) is plotted relative to the zero velocity corresponding to its rest frequency. Fits to the emission and absorption molecular lines are shown as solid lines, and, when multi-component and/or multiple lines have been fitted, they are displayed individually as dashed lines (see text for details).
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Maps of the optical depth of the OH⁺(1₁⁻₀₁) and OH⁺(1₂⁻₀₁) absorption lines, evaluated as $F_v/F_{cont} = -e^{-\tau_v}$ at the central frequency $v_0$ of each line.</figcaption></figure>

width at half maximum (FWHM), the integrated intensity, and the intrinsic frequency. All the uncertainties are propagated into the derived quantities.

### 3.1.3. Optical depths and column densities of lines seen in absorption

The optical depth of an absorption line ($\tau_v$) can be evaluated from the observed spectrum normalized by the continuum as follows

$$1 - e^{-\tau_v} = \frac{S_v^{\text{cont}} - S_v}{S_v^{\text{cont}}}, \qquad (1)$$

where $S_v$ is the flux density of the spectrum and $S_v^{\text{cont}}$ that of the continuum only. The optical depth is evaluated at the central frequency $v_0$ of the line. Table 3 reports the values of $\tau(v_0)$ obtained from the integrated spectra shown in Fig. 3, that is, averaged over the angles covered by the lines on the sky. In the case of the two strongest absorption lines (the OH⁺(1₁⁻₀₁) and (1₂⁻₀₁) transitions), maps of the opacity are displayed in Fig. 6.

Using Eq. (1) and the continuum-normalized spectra, we compute “optical depth spectra” that describe the variation of $\tau_v$ along the absorption lines, as a function of frequency or velocity. The column density $N_l$ of the absorption line is derived by integrating over the velocity range of the line along the line of sight (e.g., Comito et al. 2003; see also Indriolo et al. 2015, 2018):

$$N_l = \frac{8\pi v^3}{c^3} \frac{g_l}{A_{ul}g_u} \int \tau_v dV. \qquad (2)$$

The quantity $\int \tau_v dV$ and the speed of light $c$ can be expressed in units of [cm s⁻¹], and the transition frequency $v$ in [Hz] (see Table 2 for the corresponding frequencies of each line), so that $N_l$ is directly computed in units of [cm⁻²]. The quantities $g_l$ and $g_u$ are the lower and upper state statistical weights, and $A_{ul}$ is the spontaneous emission coefficient. We compute the column density for the strongest of the hyperfine transitions for the specific rotational $\Delta J$, and we assume that nearly all molecules are in the ground rotational state. The values of the statistical weights and the emission coefficients are taken from the Cologne Database for Molecular Spectroscopy (CDMS; Müller et al. 2005).

Table 3 summarizes the results for the absorption lines observed in HerBS-89a, i.e, the three OH⁺ transitions,

<table><caption>**Table 3.** Optical depths and column densities for the molecular absorption lines.</caption><thead><tr><th>Line</th><th>τ(v<sub>0</sub>)</th><th>ΔV [km s<sup>-1</sup>]</th><th>∫ τ<sub>v</sub>dV [km s<sup>-1</sup>]</th><th>N<sub>l</sub> [10<sup>14</sup> cm<sup>-2</sup>]</th></tr></thead><tbody><tr><td>OH<sup>+</sup>(1<sub>1</sub>-0<sub>1</sub>)</td><td>0.62 ± 0.14</td><td>-280 to +530</td><td>191.4 ± 17.2</td><td>11.2 ± 1.0</td></tr><tr><td>OH<sup>+</sup>(1<sub>2</sub>-0<sub>1</sub>)</td><td>0.61 ± 0.16</td><td>-400 to +600</td><td>306.0 ± 23.9</td><td>9.6 ± 0.8</td></tr><tr><td>OH<sup>+</sup>(1<sub>0</sub>-0<sub>1</sub>)</td><td>0.24 ± 0.09</td><td>-33 to +560</td><td>53.8 ± 15.1</td><td>7.2 ± 2.0</td></tr><tr><td>CH<sup>+</sup>(1-0)</td><td>0.36 ± 0.13</td><td>-145 to +500</td><td>91.4 ± 22.8</td><td>2.6 ± 0.7</td></tr><tr><td>HCN(11-10)</td><td>0.24 ± 0.12</td><td>-570 to +430</td><td>66.1 ± 18.3</td><td>0.38 ± 0.1</td></tr><tr><td>NH(1<sub>2</sub>-0<sub>1</sub>)</td><td>0.24 ± 0.12</td><td>-570 to +430</td><td>66.1 ± 18.3</td><td>6.7 ± 1.9</td></tr></tbody></table>

Notes. The frequency $v_0$ corresponds to the deepest dip of each lines. In the case of the absorption line observed at 246.72 GHz, we provide two possible column densities obtained assuming that the line is a transition of either HCN or NH separately, although we note that both lines likely contribute at the same time (see Sect. 8.1).

CH⁺(1-0) and the absorption line observed at 246.72 GHz that is due to HCN(11-10) and/or NH(1₂⁻₀₁), listing the velocity integration ranges, the integrated optical depths, and the column densities. Uncertainties are computed from the dispersions of the spectra, via standard error propagation.

### 3.1.4. Velocity fields

Detailed views on the dynamics of HerBS-89a are presented in Fig. 7, which displays the velocity fields of all the molecular emission and absorption lines detected in this source. The emission lines of ¹²CO(9-8) and H₂O(2₀₂-1₁₁) show similar, regular east-west velocity gradients, from blue to red along the southern arc, and reversed in the northern peak, as expected from gravitational lensing (see Sect. 5). The combined image of the ¹²CO(9-8) and H₂O(2₀₂-1₁₁) emission lines reveals further details, in particular the velocity distribution along the western arc. The velocity field of these dense gas tracers is consistent with the presence of kinematically distinct components that could indicate a merger system, or with a single rotating disk-like structure (see Sect. 7). The main peaks of the three OH⁺ and one CH⁺ absorption lines are seen to be red-shifted across most of the continuum (Fig. 7). Finally, the red-shifted component of the NH₂(2₀₂-1₁₁)(5/2-3/2) emission line is clearly seen
<figure><img src="image_8.png" /></figure>

**Fig. 7.** Velocity field maps (shown in colors) in each of the molecular emission and absorption lines detected in HerBS-89a superimposed on the continuum at the individual frequencies of each molecular line (in contours of $10\sigma$ spacing). The velocity maps were obtained above $3\sigma$ thresholds in the zero-th moment map of the molecular lines. As in Fig. 4, the *three top rows* display the high-angular resolution data, whereas the *lower row* shows the lower-angular resolution images. *Upper right panel:* combined $^{12}$CO(9–8) and H$_2$O(2$_{02}-1_{11}$) image that reveals the kinematics along the western arm of the partial Einstein ring. *Right panel in the second row:* combined velocity field of the two OH+ absorption lines ((1$_1-0_1$) and (1$_2-0_1$)). The synthesized beam is shown in the lower left corner of each of the images.

west of the continuum peak, as expected from the velocity field revealed in the higher angular resolution data for the $^{12}$CO(9–8) and H$_2$O(2$_{02}-1_{11}$) emission lines.

## 3.2. VLA results: The $^{12}$CO(1–0) emission line

Figure 5 shows the VLA $^{12}$CO(1–0) emission line spectrum, resampled to a velocity resolution of ~120 km s$^{-1}$ to highlight the key properties of the line profile. Because of the limited S/Ns,

we fit the observed line with a single Gaussian profile. The position of the $^{12}$CO(1–0) emission line is consistent with the nominal redshift of HerBS-89a. Its FWHM is $\Delta V \sim 1400$ km s$^{-1}$, comparable to those of the three high-$J$ CO emission lines, and the observed total line intensity is $I_v = 0.64 \pm 0.13$ Jy km s$^{-1}$ (Table 2).

The emission line map, which is shown in Fig. 8, was obtained by integrating the VLA data cube between –800 and +500 km s$^{-1}$ relative to the frequency of $^{12}$CO(1–0) at z =
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** $^{12}$CO(1–0) emission in HerBS-89a. The $1\sigma$ sensitivity is $32.1~\mu\text{Jy}\text{beam}^{-1}$ and the beam size (shown in the lower left corner) is $1.24 \times 0.79$ (PA $-62.3^\circ$, E of N). The emission is resolved in the east-west direction. Contours are in $1\sigma$ steps starting from the $3\sigma$ level.</figcaption></figure>

2.9497. The $^{12}$CO(1–0) line emission is slightly resolved and displays an east-west extension, similar to what is observed in the other CO lines for which the emission has been resolved (i.e., the $^{12}$CO(9–8) and (5–4) lines). It is noteworthy that, contrary to what is seen in the $^{12}$CO(9–8) and dust continuum emission, the northern component in $^{12}$CO(1–0) appears slightly stronger than the southern component; however, the difference is at a $<1\sigma$ level. If confirmed, this inconsistency may be the result of differential magnification and indicate that the molecular gas reservoir traced in CO (1–0) is more extended than that traced by the high-J CO and water transitions (see Sect. 5 for a detailed discussion).

Due to the limited sensitivity of the VLA data, the underlying continuum of HerBS-89a was not detected. We measure an rms of ~15.4 and ~20.1 $\mu$Jy beam$^{-1}$ in the continuum maps at 29 GHz (band width 0.74 GHz) and 38 GHz (band width 0.93 GHz) in the observed frame, with beam sizes of 1.23 × 0.79 and 0.96 × 0.59, respectively.

# 4. The foreground lensing galaxy

In this section, we derive the characteristics of the foreground lensing galaxy from the near-infrared HST image and the GTC photometry, supplemented by NOEMA data.

## 4.1. Imaging results (HST)

Figure 9 displays the HST wide-J filter F110W image (at $1.15~\mu\text{m}$) toward HerBS-89a compared to the 1.2 mm dust continuum emission. The HST image reveals a galaxy along the line of sight to HerBS-89a, which lies close to its northern component but is offset by ~0.2 southwest from the northern dust continuum peak. The lensed galaxy itself remains undetected in the HST image. The near-infrared emission of the foreground lensing galaxy is resolved in the HST image. It is dominated by a bright bulge-like structure. It also displays a weak emission that extends in the east-west direction and that could resemble a disk. A bright spot to the northwest of the spheroidal component seems to wrap around the northern peak of the dust continuum, and could be a small satellite galaxy. From the HST

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** View of the foreground lensing galaxy in the HST WFC3 F110W filter ($1.15~\mu\text{m}$) image, with the 254.6 GHz dust continuum measured with NOEMA superposed as red contours.</figcaption></figure>

data, the foreground galaxy has an estimated flux density of ~4.36 $\pm$ 0.88 $\mu$Jy at $1.15~\mu\text{m}$.

## 4.2. Photometric redshift (GTC)

The foreground lensing galaxy is also detected in our GTC OSIRIS r-band imaging with a magnitude $r = 24.5 \pm 0.3$ AB mag. The deeper imaging with HiPERCAM detects the foreground galaxy in the griz bands, while in the u band the object remains undetected. The target lies very close to a bright ($r = 19.5$) point-like object, 3.3 to the east, and the photometric measurement needs to take properly into account the wings of the bright neighbor.

The HiPERCAM maps were calibrated on SDSS DR12 stars present in the field. The code GALFIT (Peng et al. 2002) was used to subtract the profile of the bright point-like nearby object using a model PSF. Aperture photometry was performed on the residual image to measure the emission of the target lensing galaxy. Uncertainties were derived measuring the noise level in an annulus around the object, thus taking into account the presence of the nearby source. Finally, systematic zero-point uncertainties were propagated in the estimate of the magnitude uncertainties. Table 4 lists the derived magnitudes and the $3\sigma$ u band upper limit in the AB system (Oke & Gunn 1983).

The foreground galaxy is also detected by WISE (Wright et al. 2010) at $3.4$ and $4.6~\mu\text{m}$. However, the WISE photometry likely includes contributions from the lensing galaxy, the lensed HerBS-89a galaxy, and the nearby bright point-like source. Because of the broad WISE beam ($6.1$ at $3.4~\mu\text{m}$), it is impossible to disentangle the different components. Therefore only the GTC and HST data points have been retained in this analysis.

At $r = 24.5$ it is far from easy to obtain an optical spectroscopic redshift, but thanks to the sensitivity of HiPERCAM and the capability to observe simultaneously the five ugriz bands, together with the HST F110W imaging results, it is possible to derive a photometric redshift.
<table><caption>**Table 4.** Optical photometry of the foreground lensing galaxy.</caption><thead><tr><th>Band</th><th>λ<br>[µm]</th><th>mag<br>[AB]</th><th>Flux density<br>[µJy]</th><th>Instrument</th></tr></thead><tbody><tr><td>u</td><td>0.36</td><td>&gt;25.0</td><td>&lt;0.36</td><td>GTC HiPERCAM</td></tr><tr><td>g</td><td>0.48</td><td>26.0 ± 0.3</td><td>0.15 ± 0.05</td><td>"</td></tr><tr><td>r</td><td>0.64</td><td>24.5 ± 0.3</td><td>0.63 ± 0.06</td><td>GTC OSIRIS</td></tr><tr><td>r</td><td>0.63</td><td>24.6 ± 0.2</td><td>0.53 ± 0.11</td><td>GTC HiPERCAM</td></tr><tr><td>i</td><td>0.77</td><td>23.8 ± 0.2</td><td>1.10 ± 0.22</td><td>"</td></tr><tr><td>z</td><td>0.91</td><td>23.1 ± 0.2</td><td>2.09 ± 0.42</td><td>"</td></tr><tr><td>F110W</td><td>1.15</td><td>22.3 ± 0.2</td><td>4.36 ± 0.88</td><td>HST WFC3</td></tr><tr><td>W1</td><td>3.35</td><td>19.4 ± 0.1</td><td>63.1 ± 6.1</td><td>WISE</td></tr><tr><td>W2</td><td>4.60</td><td>19.2 ± 0.2</td><td>75.9 ± 15.3</td><td>"</td></tr></tbody></table>

Notes. Due to the contamination of a nearby point-like source, the WISE data were not considered in the derivation of the photometric redshift of the foreground lensing galaxy (see Sect. 4.2). The upper limit to the u-band flux density is given at a $3\sigma$ level.

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Best derivation of the photometric redshift of the foreground lensing galaxy, using the Le PHARE code (Ilbert et al. 2006; Arnouts et al. 1999). The filled triangle shows the HST F110W flux density and the filled circles represent the GTC measurements (see text). The WISE data (open circles) have not been included in the fit. The curve shows the best fit model, a SB spiral template (Polletta et al. 2007; Ilbert et al. 2009). The inset displays the distribution of $\chi^2$ as a function of redshift; the dotted horizontal lines correspond to the 1, 2, $3\sigma$ levels, defined as $\chi^2_{min} + 1, 2.71, 6.63$. The best photometric redshift for the foreground galaxy is $z_{phot} = 0.9^{+0.3}_{-0.5}$ where the errors are $3\sigma$.</figcaption></figure>

To get the best possible result, and test systematics due to the chosen method, we applied various photometric redshift codes EAZY (Brammer et al. 2008), Hyper-z (Bolzonella et al. 2000), Le PHARE (Arnouts et al. 1999; Ilbert et al. 2006), and MAGPHYS-z (da Cunha et al. 2008; Battisti et al. 2019) to the ugriz and F110W data, adopting different template libraries (Coleman et al. 1980; Bolzonella et al. 2000; Berta et al. 2003, 2013; Babbedge et al. 2004; Polletta et al. 2007; Brammer et al. 2008; Ilbert et al. 2009; da Cunha et al. 2015). Since it best accounts for the upper limits, the final choice is the Le PHARE code.

The best photometric redshift of the foreground lensing galaxy thus obtained is $z_{phot} = 0.9$, with a 99% confidence interval of $0.4 < z_{phot,3\sigma} < 1.2$. The best fit solution is obtained with a SB spiral template from the COSMOS library (Ilbert et al. 2009), originally belonging to the SWIRE library (Polletta et al. 2007). Figure 10 shows the best fit model together with the

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Spectrum between 139.29 and 147.04 GHz toward HerBS-89a, from Neri et al. (2020), showing the $^{12}$CO(5–4) emission line of HerBS-89a (at 145.902 GHz) and the weaker, narrow (~200 km s$^{-1}$) emission line at 143.74 GHz, possibly belonging to the foreground lensing galaxy (see Sect. 4.3 for further details).</figcaption></figure>

$\chi^2$ distribution. The photometric redshift is best constrained by the r-i color, tracing the D4000 break$^4$. The lack of rest-frame near-infrared coverage and the upper limit in the u band restrict the accuracy of the photo-z measurement. Deeper data and a wider wavelength coverage would be necessary to improve its precision.

In addition, various attempts were made using templates that include emission lines (e.g., Babbedge et al. 2004), with the aim to verify if the r and F110W bands could be contaminated by lines at $z \sim 0.7$, namely [OII] ($\lambda3727$) and H$\alpha$ ($\lambda6563$). However, none produced a better fit to the data.

### 4.3. Spectroscopic redshift (NOEMA)

An independent indication of the redshift of the lensing foreground can possibly be derived from the 2 mm spectral scan that was obtained to derive the spectroscopic redshift of HerBS-89a in Neri et al. (2020). In the frequency range 139.3–147.0 GHz, where the $^{12}$CO(5–4) emission line of HerBS-89a was detected, a weaker narrow (~200 km s$^{-1}$) emission line is seen at 143.74 GHz, close to the center position of HerBS-89a (see Fig. 11). The narrowness of this emission line excludes the possibility that it could be related to HerBS-89a. If it traces instead the foreground lensing galaxy, the line could correspond to a $^{12}$CO emission line, either detected in the $J = 2-1$ transition, in which case $z_{lens} = 0.604$, or in the $J = 3-2$ transition for which $z_{lens} = 1.406$. The former is consistent within the $2\sigma$ uncertainties with the photometric redshift derived from the GTC and HST data. A search for the counterpart of this emission line (at 1 mm) could allow us to verify the above identification, and eventually derive a precise spectroscopic redshift of the foreground galaxy gravitationally amplifying HerBS-89a.

## 5. Lens modeling

In order to recover the source-plane morphology and the intrinsic properties of HerBS-89a, we derive a lens model from the high-angular resolution NOEMA continuum data as shown in Fig. 2 (with initial lens parameters informed by HST

<footnote>$^4$ The 4000 Å discontinuity is defined as the ratio between the average flux density between 4050 and 4250 Å and that between 3750 and 3950 Å (Bruzual 1983).</footnote>
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Lens modeling results for the 1.2 mm dust continuum emission and spectral lines in HerBS-89a (identified in the right most panels). In each row, from left to right: observed image (with the synthesized beam shown in the lower left corner) and critical curve (red line); reconstructed model image convolved with the PSF; residuals; reconstructed model image at full angular resolution; and (separated and with a different angular scale) reconstructed source-plane image and caustic. The residuals are within ±2.5σ, showing that the model agrees well with the data. The cross shows the coordinates of the phase center of NOEMA observations (RA 13:16:11.52 and Dec +28:12:17.7 in J2000). The OH+ absorption lines have been inverted to allow for a positive normalization of the Sérsic profile reconstruction.</figcaption></figure>

observations) and proceed to apply it to the molecular line data. We perform the lens modeling and source reconstruction using the *lensmodel* package (Keeton 2001) with the *pixsrc* extension (Tagore & Keeton 2014; Tagore & Jackson 2016). For the sake of computational efficiency, all modeling has been done in the sky plane, although we recognize that fitting directly to

visibility data can prevent specific choices in the interferometric imaging process (e.g., related to uv weighting and deconvolution) from having undue impact on a model. Although the quality of our NOEMA observations is excellent, the limited number of resolution elements across the source prevents us from leveraging the full power of the nonparametric approach of
<table><caption>**Table 5.** Results of lens modeling: parameters of the best fit model for the foreground lens.</caption><thead><tr><th rowspan="2">Parameter</th><th rowspan="2">Value</th><th colspan="2">Uncertainties</th></tr><tr><th>1σ</th><th>Systematic</th></tr></thead><tbody><tr><td>Ref. point (RA)</td><td>13:16:09.79</td><td>—</td><td>—</td></tr><tr><td>Ref. point (Dec)</td><td>28:12:18.17</td><td>—</td><td>—</td></tr><tr><td>Position ΔRA</td><td>-0′.0612</td><td>0.0005</td><td></td></tr><tr><td>Position ΔDec</td><td>-0′.0989</td><td>0.0008</td><td></td></tr><tr><td>Einstein radius <i>b</i></td><td>0′.4832</td><td>0.0006</td><td>0.0483</td></tr><tr><td>Ellipticity <i>e</i></td><td>0.177</td><td>0.002</td><td>0.018</td></tr><tr><td>Position angle θ<sub>e</sub></td><td>32°4</td><td>0.4</td><td>3.2</td></tr></tbody></table>

**Notes.** The position angle's orientation is measured east of north.

Tagore & Keeton (2014), instead limiting us to the use of parametric (Sérsic profile) components in reconstructing the sources.

Because we are not using regularization to guide the reconstruction of nonparametric sources, correlations in the noise introduced by the finite uv sampling do not require special handling à la Riechers et al. (2008). In all cases, we convolve the lensed model images (fourth column of Fig. 12) with the appropriate (two-dimensional Gaussian) synthesized beams to generate the reconstructed model images (second column of Fig. 12) that are directly compared to the data.

The foreground deflector is assumed to have a singular isothermal elliptical mass distribution (SIEMD) described by five parameters: Einstein radius (b), RA, Dec, ellipticity (e)<fnref n="5" />, and ellipticity position angle (θe; measured east of north). Table 5 lists the best fit values of these parameters and their uncertainties as determined from a Markov chain Monte Carlo (MCMC) exploration of parameter space. The position of the lens is given relative to the reference position RA = 13:16:09.79 and Dec = 28:12:18.17. Because of the very high S/N of the continuum data, we use the lens model derived from the continuum data when reconstructing the spectral line data.

The left panels of Fig. 12 include the critical curve (i.e., the locus of infinite magnification in the image plane; red lines). In the case of a SIEMD with no external shear, this curve is also an iso-density contour of the deflector's mass distribution whose position and geometry can be compared to the HST image of the lensing galaxy (Fig. 9). The corresponding source-plane caustic is plotted in the right panels of Fig. 12.

The continuum and (almost) all of the spectral line zeroth moment maps are fit well by single Sérsic components in the source plane (see Table 6 and Fig. 12 for fitting results, and below for more discussion of the line data). The profile is parametrized in terms of its position, normalization $I_0$, ellipticity $e$, position angle $θ_e$, scale radius $R_s$, and Sérsic index $n$:

$$I(r) = I_0 \exp \left[ - \left( \frac{r}{R_s} \right)^{1/n} \right]. \qquad (3)$$

Using the estimated photometric redshift of the foreground galaxy and its uncertainty (see Sect. 4.2), we estimate that the deflector mass enclosed within the Einstein radius is $M_{\text{lens}} = 9.8^{+3.2}_{-4.9} \times 10^{10} M_{\odot}$ (where the upper and lower bounds are at $3\sigma$).

Statistical uncertainties in the lens and continuum source model were evaluated with two MCMC cycles. The first fixed the source model and determined the uncertainties of the lens parameters; the second fixed the lens model to its optimal parameters

and let the source parameters vary. This approach was chosen for computing speed and for ease of comparison between continuum and line results (because the latter assume the lens model derived from the continuum data, i.e., the lens model is fixed in this case). However, it should be noted that the MCMC approach underestimates the lens+source errors. Appendix A presents the conditional posterior distributions of the lens and source model parameters.

To these statistical uncertainties, we add fiducial systematic uncertainties of 10% associated with the choices of lensing potential, source brightness profile, interpolation errors, etc. The only exceptions are the source model positional errors, for which we adopt systematic uncertainties given by

$$\frac{\sqrt{\text{HPBW}_{\text{min}} \text{HPBW}_{\text{maj}}}}{2} \frac{\sigma_I}{I_{\text{line}}} \frac{1}{\sqrt{\mu}}, \qquad (4)$$

in terms of the synthesized beam dimensions $\text{HPBW}_{\text{min}}$ and $\text{HPBW}_{\text{max}}$, the S/N of the observed detection $\sigma_I/I_{\text{line}}$, and the magnification $\mu$ (see, e.g., Downes et al. 1999). For the lens model, we estimate systematic positional errors as the astrometric errors from the continuum imaging. Because we expect systematic errors to dominate and because of computational constraints, we do not perform MCMC on the line spectral source parameters.

As mentioned above, nearly all of the spectral line data are fit well by single Sérsic components; the exception is the $H_2O(2_{02}-1_{11})$ zeroth-moment map, which shows faint (~5σ, compared to the ~10σ peak) emission to the west and northwest that motivated some exploration of a two-component model. We adopt here a single-component fit for this line, due to numerical instability in the two-component fit, but there are some systematic features in the residuals associated with the peak of the emission. The results are listed in Table 6, along with the corresponding magnification factors $\mu$, for the high resolution data and for $^{12}\text{CO}(1-0)$. Because of the lower resolutions and S/Ns, we do not model the spectral lines and continuum that were measured at lower frequencies (from 209 to 232 GHz).

The magnification factor $\mu$ for the continuum and line data was computed using samples from the lens and source parameter MCMC runs. We estimate a magnification of $5.05 \pm 0.03(\pm 0.51)$ for the 254.6 GHz dust continuum emission, where the term in parentheses represents the 10% systematic uncertainty.

In the case of the $^{12}\text{CO}(1-0)$ line, we infer a lensing magnification $\mu = 5.0$, with a large systematic uncertainty mainly due to possible different choices in how the zeroth-moment map is derived (also reflecting the limited S/N of the current data set). Deeper and higher angular resolution observations of the $^{12}\text{CO}(1-0)$ emission would be helpful to better constrain the amplification factor as well as to study in greater detail the morphology (extent and kinematics) of the cold gas reservoir of HerBS-89a.

## 6. Global properties of HerBS-89a: Dust and molecular gas

In this section, we derive the intrinsic properties of HerBS-89a from the new millimeter, centimeter, and optical/near-infrared data presented in Sect. 3, and the lensing model described in Sect. 5. The discussion is organized as follows: the analysis of the continuum spectral energy distribution and the derivation of the molecular gas mass and dust-to-gas ratio are described in Sects. 6.1 and 6.2; the CO line excitation is investigated in

<footnote n="5">⁵ Ellipticity is defined as $e = 1 - q$, where $q$ is the axis ratio.</footnote>
<table><caption>**Table 6.** Results of lens modeling: parameters of reconstructed source-plane emission (and absorption) in HerBS-89a.</caption><thead><tr><th rowspan="2">Tracer</th><th rowspan="2">I<sub>0</sub><br>[Jy beam<sup>-1</sup>]</th><th colspan="2">ΔRA</th><th rowspan="2">e</th><th rowspan="2">θ<sub>e</sub><br>[deg]</th><th colspan="2">R<sub>s</sub></th><th rowspan="2">n</th><th rowspan="2">μ</th></tr><tr><th colspan="2">[mas]</th><th>[mas]</th><th>[pc]</th></tr></thead><tbody><tr><td>Continuum <sup>(a)</sup><br>syst. unc.</td><td>0.16 ± 0.01<br>(0.02)</td><td>-1.9 ± 0.2<br>(9.5)</td><td>-277.8 ± 0.4<br>(9.5)</td><td>0.522 ± 0.007<br>(0.052)</td><td>-97.4 ± 0.4<br>(9.7)</td><td>12 ± 1<br>(1.0)</td><td>95 ± 8<br>(10.0)</td><td>1.82 ± 0.04<br>(0.18)</td><td>5.05 ± 0.03<br>(0.51)</td></tr><tr><td>CO(9–8)<br>syst. unc.</td><td>0.30<br>(0.03)</td><td>-9.6<br>(11.2)</td><td>-254.5<br>(11.2)</td><td>0.408<br>(0.041)</td><td>-92.8<br>(9.3)</td><td>0.4<br>(0.1)</td><td>3.0<br>(1.0)</td><td>3.21<br>(0.32)</td><td>5.23<br>(0.52)</td></tr><tr><td>CO(1–0)<br>syst. unc.</td><td>0.010<br>(0.001)</td><td>15.35<br>(45.0)</td><td>-248.3<br>(45.0)</td><td>0.189<br>(0.019)</td><td>87.8<br>(8.8)</td><td>0.8<br>(0.1)</td><td>5.9<br>(0.6)</td><td>4.08<br>(0.41)</td><td>5.0<br>(1.0)</td></tr><tr><td>H<sub>2</sub>O(2<sub>02</sub>–1<sub>11</sub>)<br>syst. unc.</td><td>0.008<br>(0.001)</td><td>-19.1<br>(22.5)</td><td>-312.4<br>(22.5)</td><td>0.559<br>(0.056)</td><td>-79.7<br>(8.0)</td><td>50.0<br>(5.0)</td><td>396.0<br>(40.0)</td><td>1.50<br>(0.15)</td><td>3.87<br>(0.39)</td></tr><tr><td>OH<sup>+</sup>(1<sub>1</sub>–0<sub>1</sub>)<br>syst. unc.</td><td>0.044<br>(0.004)</td><td>-3.3<br>(7.0)</td><td>-287.1<br>(7.0)</td><td>0.045<br>(0.005)</td><td>-23.9<br>(2.4)</td><td>26.0<br>(3.0)</td><td>206.0<br>(24.0)</td><td>0.861<br>(0.086)</td><td>4.69<br>(0.47)</td></tr><tr><td>OH<sup>+</sup>(1<sub>2</sub>–0<sub>1</sub>)<br>syst. unc.</td><td>0.042<br>(0.004)</td><td>-2.7<br>(9.4)</td><td>-314.1<br>(9.4)</td><td>0.271<br>(0.027)</td><td>-68.6<br>(6.9)</td><td>29.0<br>(3.0)</td><td>230.0<br>(24.0)</td><td>1.113<br>(0.111)</td><td>4.24<br>(0.42)</td></tr></tbody></table>

**Notes.** <sup>(a)</sup>The dust continuum emission is measured at 254.5 GHz, and the uncertainties are ±1σ. For each parameter, the systematic uncertainties are tabulated separately in parentheses on a second line (see Sect. 5). In the heading of the table, *e* stands for the ellipticity, *θ*<sub>e</sub> for the position angle (measured east of north), *R*<sub>s</sub> for the scale radius, *n* for the Sérsic index, and *μ* for the amplification factor. Surface brightness units of Jy beam<sup>-1</sup> are inherited from the surface brightness units of the data. Strong covariance between *R*<sub>s</sub> and *n* (see, e.g., Fig. A.1) means that the <sup>12</sup>CO(9–8) and <sup>12</sup>CO(1–0) data are formally compatible with (larger) *R*<sub>s</sub> and (smaller) *n* like those recovered for the other reconstructions.

Sect. 6.3, and the properties of the water emission line are examined in Sect. 6.4.

## 6.1. Continuum spectral energy distribution

The new NOEMA continuum measurements (see Table 1), combined with the data from Neri et al. (2020) and Herschel and SCUBA-2 photometry (Bakx et al. 2018, 2020), provide an exquisite sampling of the far-infrared (FIR) and sub-mm spectral energy distribution (SED) of HerBS-89a. A total of 16 broad band measurements cover the range from 250 µm to 3 mm in the observed frame, 12 of which are spread along the Rayleigh-Jeans tail.

Figure 13 shows the continuum data available for HerBS-89a as well as the continuum measurements extracted for the nearby source HerBS-89b. The latter source shows an increase of flux at wavelengths above ~2 mm with respect to the dust emission, which could be due to non-thermal processes and rules out the possibility that HerBS-89b is a third lensed image of HerBS-89a.

Following and extending the work by Neri et al. (2020) and Berta et al. (2016), we model the SED of HerBS-89a using three different approaches: (i) the Draine & Li (2007, hereafter DL07) dust model; (ii) a single-temperature modified black body (MBB) in the optically-thin approximation; and (iii) an MBB model in its general form.

In the case of the DL07 modeling, interstellar dust is described as a mixture of carbonaceous and amorphous silicate grains, whose size distributions are chosen to reproduce different observed extinction laws. We defer to Draine & Li (2001, 2007), Li & Draine (2001), and Berta et al. (2016) for a detailed description of the model and its implementation. Despite the richness of the long-wavelength data, the SED is still lacking measurements in the mid-infrared, so that the fine details of the DL07 model, for example, the possible contribution of the poly-cyclic aromatics hydrocarbons (PAHs), are not constrained. The dust mass estimate is nevertheless robust, because it is dominated by the colder dust components.

For the optically thin MBB case, the emergent luminosity from a given dust mass $M_{\text{dust}}$ is described as:

$$L_{\nu} \sim M_{\text{dust}}\kappa_{\nu}B_{\nu}(T_{\text{dust}}), \quad (5)$$

where $B_{\nu}(T_{\text{dust}})$ is the Planck function, $T_{\text{dust}}$ is the dust temperature, and $\kappa_{\nu} = \kappa_{0}(\nu/v_{0})^{\beta}$ is the mass absorption coefficient of dust at rest frequency $\nu$.

The general form of the MBB differs from the optically thin approximation by the factor $(1 - e^{-\tau_{\nu}})$, with $\tau_{\nu} = \tau_{0}(\nu/v_{0})^{\beta}$, instead of the $\nu^{\beta}$ term. Here $v_{0}$ is the frequency at which $\tau = \tau_{0} = 1.0$ (see Berta et al. 2016). Therefore, the general form of the MBB tends to produce best fit solutions with higher dust temperature than the optically thin case. In the optically thin approximation, three free parameters are at play: $T$, $\beta$ and the model normalization. In the general form, $v_{0}$ is also a free parameter.

The MBB fits are limited to the data with rest-frame wavelength $\lambda_{\text{rest}} > 50$ µm, thus avoiding biases toward warmer temperatures. The effects of the cosmic microwave background (CMB) have been taken into account, following the prescription of da Cunha et al. (2013). A correction similar to that derived for the MBB has also been applied to the DL07 model.

Berta et al. (2016) and Bianchi (2013) present thorough discussions of the proper use of $\kappa_{\nu}$ and $\beta$. Appendix B summarizes the current $\kappa_{\nu}$ panorama. Following the choice often found in the literature, we adopt the Draine (2003) $\kappa_{\nu}$ and apply the correction prescribed by Draine et al. (2014): $\kappa_{\nu} = 0.047 \text{ m}^2 \text{ kg}^{-1}$ at 850 µm. The same correction is also applied to the dust masses derived with the DL07 approach.

The best fit models are shown in Fig. 13, and the results are listed in Table 7. We include all the values obtained with our fits using the continuum flux densities not corrected for the lens magnification, as well as our preferred source-plane values for the infrared luminosity and the dust mass based on the DL07 fit corrected for gravitational magnification (see Sect. 5). The CMB effect produces a steepening of the millimetric slope by ~+0.15, comparable to the uncertainty on $\beta$. At the redshift $z = 2.95$ and for the dust temperature of HerBS-89a, the effect is milder than what found at $z = 3.6–5.8$ by Jin et al. (2019).

Finally, the VLA $3\sigma$ upper limits for HerBS-89a that are shown in Fig. 13 have been computed over an area equivalent to the NOEMA 263.3 GHz continuum (see Table 1). For comparison, a power-law spectrum $S_{\nu} \propto \nu^{-0.8}$ is plotted,
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Spectral energy distribution. Top panel: observed SED of HerBS-89a (filled symbols) and 89b (open symbols). The HerBS-89a data are fit with three different dust models: DL07 (blue lines), optically-thin MBB (green lines); MBB in general form (red lines). The dotted lines represent the models as they are, and the dashed lines include the correction for the effects of CMB (see text for details). The light-blue continuous lines represent the synchrotron emission with a spectral index $\alpha = -0.8$, normalized on the basis of the radio-FIR correlation subject to the redshift evolution found by Magnelli et al. (2015, upper line) and Delhaize et al. (2017, lower line). The data include SPIRE flux densities (black dots; from Bakx et al. 2018); SCUBA-2 photometry (blue dot; from Bakx et al. 2020); the NOEMA derived continuum flux densities at 1 mm continuum (green dots) from this work and 2 and 3 mm (red dots) from Neri et al. (2020); and 3$\sigma$ upper limits at 29 and 38 GHz (10 and 7.8 mm, respectively) are derived from the VLA data (arrows). The data for HerBS-89b are taken from this paper for 1 mm and Neri et al. (2020) for 2 mm flux densities. Bottom panel: HerBS-89a NOEMA spectra (Neri et al. 2020, and this work) compared to the general form MBB model with CMB correction.</figcaption></figure>

representing a synchrotron emission (e.g., Ibar et al. 2010; Thomson et al. 2014), normalized on the basis of the radio-far-infrared correlation (Magnelli et al. 2015; Delhaize et al. 2017).

## 6.2. Molecular gas mass and gas-to-dust ratio

The observed $^{12}$CO(1–0) emission line intensity can be translated into a luminosity (in K km s$^{-1}$ pc$^2$), from which an estimate

<table><caption>**Table 7.** Properties of HerBS-89a derived from the dust continuum spectrum.</caption><thead><tr><th>Quantity</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td colspan="3">DL07 model</td></tr><tr><td>μL<sub>DL07</sub>(8 – 1000 μm)</td><td>(2.89 ± 0.25) 10<sup>13</sup></td><td>L<sub>☉</sub></td></tr><tr><td>μM<sub>DL07</sub> (dust)</td><td>(1.30 ± 0.07) 10<sup>10</sup></td><td>M<sub>☉</sub></td></tr><tr><td colspan="3">Optically thin MBB</td></tr><tr><td>T<sub>MBB,thin</sub> (dust)</td><td>28.0 ± 1.3</td><td>K</td></tr><tr><td>β<sub>MBB,thin</sub></td><td>2.01 ± 0.11</td><td>–</td></tr><tr><td>μL<sub>MBB,thin</sub>(50 – 1000 μm)</td><td>(1.94 ± 0.08) 10<sup>13</sup></td><td>L<sub>☉</sub></td></tr><tr><td>μM<sub>MBB,thin</sub> (dust)</td><td>(1.11 ± 0.06) 10<sup>10</sup></td><td>M<sub>☉</sub></td></tr><tr><td colspan="3">MBB in general form</td></tr><tr><td>T<sub>MBB,gen</sub> (dust)</td><td>49.2 ± 2.7</td><td>K</td></tr><tr><td>β<sub>MBB,gen</sub></td><td>2.24 ± 0.13</td><td>–</td></tr><tr><td>ν<sub>0</sub>, MBB, gen</td><td>1214 ± 121</td><td>GHz</td></tr><tr><td>λ<sub>0</sub>, MBB, gen</td><td>247 ± 24</td><td>μm</td></tr><tr><td>μL<sub>MBB,gen</sub>(50 – 1000 μm)</td><td>(2.00 ± 0.08) 10<sup>13</sup></td><td>L<sub>☉</sub></td></tr><tr><td>μM<sub>MBB,gen</sub> (dust)</td><td>(6.06 ± 0.71) 10<sup>9</sup></td><td>M<sub>☉</sub></td></tr><tr><td colspan="3">De-magnified values from DL07</td></tr><tr><td>L<sub>IR</sub>(8 – 1000 μm)</td><td>(4.6 ± 0.4) 10<sup>12</sup></td><td>L<sub>☉</sub></td></tr><tr><td>M<sub>dust</sub></td><td>(2.6 ± 0.2) 10<sup>9</sup></td><td>M<sub>☉</sub></td></tr></tbody></table>

**Notes.** The MBB results include the correction for the CMB effect on T and β. The last two rows are based on the DL07 model and have been corrected for the gravitational lens magnification.

of the molecular gas mass can be made, using the following relation (Solomon & Vanden Bout 2005):

$$L'_{\text{Line}} = 3.25 \times 10^7 S_{\text{Line}} \Delta V \times \frac{D_{\text{L}}^2}{(1+z)} \frac{1}{v_{\text{rest}}^2}, \quad (6)$$

where $S_{\text{Line}}\Delta V = I_{\text{Line}}$ is the integrated line flux in units of Jy km s$^{-1}$; $v_{\text{rest}}$ is the line rest frequency in GHz; and $D_{\text{L}}$ is the luminosity distance in Mpc.

The line luminosity, expressed in units of $L_{\odot}$, can be derived as:

$$L_{\text{Line}} = 3.2 \times 10^{-11} v_{\text{rest}}^3 L'_{\text{Line}}. \quad (7)$$

We define the molecular and total gas masses as

$$M_{\text{mol}} = M_{\text{H}_2} + M_{\text{He}}, \quad (8)$$

$$M_{\text{gas}} = M_{\text{mol}} + M_{\text{HI}}, \quad (9)$$

where $M_{\text{He}}$ is the mass of helium, $M_{\text{H}_2}$ is the mass of molecular hydrogen and $M_{\text{HI}}$ the mass of atomic hydrogen. Hereafter, we assume that $M_{\text{mol}} \gg M_{\text{HI}}$ and thus that $M_{\text{gas}} \approx M_{\text{mol}}$.

The molecular hydrogen gas mass is computed from the $^{12}$CO(1–0) luminosity applying a conversion factor $\alpha_{\text{CO}}$ in units of $M_{\odot}$ (K km s$^{-1}$ pc$^2$)$^{-1}$:

$$M_{\text{mol}} = \alpha_{\text{CO}} L'_{\text{CO}(1-0)}. \quad (10)$$

In the Milky Way and nearby star-forming galaxies with near-solar metallicity, the empirical conversion factor $\alpha_{\text{CO}}$ is $\alpha_{\text{MW}} = 4.4 \pm 0.9$, already including a contribution due to helium (see, e.g., Magnelli et al. 2012; Bolatto et al. 2013; Carilli & Walter 2013; Tacconi et al. 2020). For extreme local starbursts, there is good evidence from dynamical arguments that $\alpha_{\text{CO}}$ is 0.8–1.5 (e.g., Downes & Solomon 1998; Daddi et al. 2010; Genzel et al. 2010), as has also been suggested that for DSFGs for which $\alpha_{\text{CO}} = 0.8$ (Carilli & Walter 2013). Applying a
1.36× correction factor to account for the helium contribution, the latter value becomes $\alpha_{CO} = 1.09$. Dunne et al. (2021) used the far-infrared continuum of a small but statistically complete sample of 12 low redshift galaxies together with CO and CI measurements to determine the CO-H$_2$ conversion factor $\alpha_{CO}$. In another paper, Dunne et al. (in prep.) combined this sample with others measurements found in the literature for a total of more than 70 far-infrared and sub-mm selected galaxies with CO and CI detection and 230 with CO only. In these studies, a value of $\kappa_{850} = 0.065 \text{ m}^2 \text{ kg}^{-1}$ is adopted (in agreement with the results by Planck Collaboration XXI 2011, see Appendix B for details) and it is found that a value $\alpha_{CO} \approx 3$ is consistent with the behavior of the three gas tracers in a diverse population of galaxies, ranging from high-z powerful star-forming galaxies to local disks.

For the computation of the molecular gas mass of HerBS-89a, we here adopt $\alpha_{CO} = 3.0 M_{\odot} (\text{K km s}^{-1} \text{pc}^{-1})^{-1}$, and we multiply it by the usual factor 1.36 to account for the helium contribution. Possible consequences of using a different value are discussed at the end of this section.

The resulting line luminosity and molecular gas mass, corrected for lensing magnification, are thus

$$L'_{\text{CO(1-0)}} = (5.15 \pm 1.05) \times 10^{10} \text{ K km s}^{-1} \text{ pc}^2, \quad (11)$$

$$M_{\text{mol}} = (2.1 \pm 0.4) \times 10^{11} M_{\odot}. \quad (12)$$

The molecular gas mass derived from the $^{12}$CO(1-0) emission line and the SED-based estimate of the dust mass (see Sect. 6.1) provide the basis for deriving the average gas-to-dust ratio ($\delta_{GDR}$) in HerBS-89a:

$$\delta_{\text{GDR}} = \frac{M_{\text{gas}}}{M_{\text{dust}}}, \quad (13)$$

where the assumption is made that the two estimates are representative for the whole galaxy, and that $M_{mol} \gg M_{HI}$.

The data lead to a value of $\delta_{GDR}$ between 80 and 174, depending on whether the DL07 or the MBB (general form) dust mass estimate is used, with an uncertainty of $\pm 21$ based on error propagation. Interestingly, this value is consistent within the errors with those seen in typical star-forming galaxies of solar metallicity, near the so-called “Main Sequence of star formation” (MS), for which $\delta_{GDR} \sim 100$ (e.g., Leroy et al. 2011; Magdis et al. 2012; Rémy-Ruyer et al. 2014). This result is nevertheless subject to the adopted value for $\alpha_{CO}$.

The total molecular gas mass ($M_{mol}$) and SFR relate to each other via the so-called depletion timescale (in years):

$$\tau_{\text{depl}} = \frac{M_{\text{mol}}}{\text{SFR}}, \quad (14)$$

where the rate of star formation is computed as SFR = $1.09 \times 10^{-10} L_{8-1000 \mu\text{m}}$ with the infrared luminosity expressed in units of $L_{\odot}$ (Kennicutt 1998a, after modification for a Chabrier 2003 IMF).

Assuming that the observed infrared luminosity is entirely associated with the ongoing star formation activity of HerBS-89a, we derive an intrinsic star formation rate SFR = 614 $\pm 59 M_{\odot} \text{ yr}^{-1}$ and a depletion timescale $\tau_{depl} = (3.4 \pm 1.0) \times 10^8$ years. For its current formation rate, HerBS-89a would exhaust its molecular gas reservoir in only 340 million years, not considering any mass return to the interstellar medium or further gas inflow (see Sect. 9).

The value for $M_{mol}$ was derived from the $^{12}$CO(1-0) luminosity using $\alpha_{CO} = 3.0 \times 1.36$, a value that is close to the Milky Way value $\alpha_{MW} = 4.4 M_{\odot} (\text{K km s}^{-1} \text{pc}^{-1})^{-1}$, commonly used for

MS galaxies, rather than $\alpha_{SB} = 1.09$, often associated to starbursts or outliers in the main sequence. This choice leads to $\delta_{GDR}$ and $\tau_{depl}$ values that would position HerBS-89a close to or on the MS star-forming galaxies on the integrated Kennicutt-Schmidt relation of star formation (e.g., Lada 2015; Liu et al. 2019), which suggests that HerBS-89a is dominated by a “secular” star-forming mode.

However, despite the evolution of the MS normalization (increasing as a function of redshift, e.g., Elbaz et al. 2011; Whitaker et al. 2014; Schreiber et al. 2015), and because of the bending of the MS at the high stellar mass end, even at $z = 3$, a SFR $\sim 600 M_{\odot} \text{ yr}^{-1}$ would hardly be associated with a galaxy on the main sequence. It is therefore more likely that HerBS-89a is located in between the MS and the starbursts region. Using $\alpha_{SB}$, the gas-to-dust ratio and the depletion time scale of HerBS-89a would both decrease by a factor of 3.74, yielding $\delta_{GDR} \approx 22-47$ and $\tau_{depl} = 9.1 \pm 2.7 \times 10^7 \text{ yr}$. Realistically, the actual value of $\alpha_{CO}$, on which these results depend, is likely between the two above mentioned extremes.

Further insights into the position of HerBS-89a with respect to the MS would require the knowledge of its stellar mass $M^*$. Detecting HerBS-89a and distinguishing it from its deflector seem out of reach for currently existing optical/NIR facilities, but the James Webb Space Telescope (JWST) could provide the needed sensitivity.

### 6.3. $^{12}$CO spectral line energy distribution

Together with the measurements reported in Neri et al. (2020), we have detected in total four $^{12}$CO emission lines in HerBS-89a, namely the $J = 1-0, 3-2, 5-4$ and $9-8$ transitions. The $^{12}$CO Spectral Line Energy Distribution (SLED) of HerBS-89a, normalized to the $J = 1-0$ transition after correcting for the lens magnification (see Table 6), is shown in Fig. 14. The velocity-integrated fluxes of the $^{12}$CO lines increase monotonically with rotational quantum number up to $J = 5-4$, with a turnover around that transition, before decreasing at higher $J$ values, with the $J = 9-8$ emission line about a factor of 8 times weaker than the $J = 5-4$ line. The peak of the CO SLED for HerBS-89a should therefore occur around the $J = 5-4$ and $J = 6-5$ transitions. This result makes HerBS-89a comparable to other starburst galaxies as shown in Fig. 14, where various SLEDs of several individual DSFGs and one quasar (the Cloverleaf) are displayed together with the center of M 82 and the Milky Way disk. In addition, the average SLEDs of three samples are shown: The (mostly) unlensed sample of 32 DSFGs selected at $850 \mu\text{m}$, at $1.2 < z < 4.1$ from Bothwell et al. (2013), the sample of 22 (mostly lensed) high-z ($2.0 < z < 5.7$) DSFGs selected from the SPT survey (Spilker et al. 2014), and the sample of 78 “main sequence” star-forming galaxies from Valentino et al. (2020). The CO SLED of HerBS-89a most closely resembles that of the Cosmic Eyelash, which has a similar peak around $J_{up} \sim 5-6$ and a falloff toward higher $J_{up}$. It is also remarkably similar to the CO SLED of the starburst GN20 (Carilli et al. 2010), up to $J_{up} = 6$; however, recent observations by Cortzen et al. (2020, and in prep.) report an increase of the GN20 CO SLED at $J_{up} = 7$.

To investigate the CO line excitation further and constrain the physical conditions of the molecular gas in HerBS-89a, we model the CO line fluxes using a Large Velocity Gradient (LVG) statistical equilibrium method (following the approach described in Yang et al. 2017, and references therein). The free parameters are the kinetic temperature of the molecular gas ($T_{kin}$), the volume density ($n_{H_2}$), the column density of CO per unit velocity
<figure><img src="image_15.png" /><figcaption>Fig. 14. $^{12}$CO Spectral Line Energy Distribution (SLED) for HerBS-89a (red dots), compared to other well-studied local and high-z galaxies (or samples of galaxies). The CO SLEDs are normalized to the $J=1-0$ transition, after correcting for the magnification, and the line fluxes of HerBS-89a have been corrected for amplification (see Sect. 5 and Table 6) before normalizing. The three samples of galaxies are from Bothwell et al. (2013) including 32 DSFGs at $1.2 < z < 4.1$, Spilker et al. (2014), including 22 SPT-selected DFSGs at $2.0 < z < 5.7$, and Valentino et al. (2020), including 78 “main sequence” star-forming galaxies. References to individual sources are: Cosmic Eyelash (Swinbank et al. 2011; Danielson et al. 2011); Cloverleaf (Barvainis et al. 1997; Weiß et al. 2003; Bradford et al. 2009; Riechers et al. 2011); GN20 (Carilli et al. 2010; Cortzen et al. 2020, and in prep.); M82 (Weiß et al. 2005); and the Milky Way disk (Fixsen et al. 1999).</figcaption></figure>

gradient ($N_{CO}/dV$), and the solid angle ($\Omega_{app}$) of the source. The CO SLED only depends on $T_{kin}$, $n_{H_2}$ and $N_{CO}/dV$. A Bayesian approach is used to fit the line fluxes generated from the model and the code *emcee* is adopted to perform a MCMC calculation. The resulting single component fit to the CO SLED of HerBS-89a is shown in Fig. 15. The SLED peaks around $J_{up} = 6-5$, and, from the single excitation component fitting, the molecular gas density is found to be $\log_{10}(n_{H_2}/\text{cm}^{-3}) = 2.79_{-0.54}^{+0.72}$, the kinetic temperature $T_{kin} = 174_{-92}^{+157}$ K and the CO column density $\log_{10}(N_{CO}/dV/(\text{cm}^{-2}\text{km s}^{-1})) = 17.59_{-0.59}^{+0.55}$, where the errors are $3\sigma$.

A second fit including an additional, colder and more extended component is shown in the lower panel of Fig. 15. Because of the paucity of available data points, the cold component is poorly constrained, with $\log_{10}(n_{H_2}) = 2.30_{-0.57}^{+1.06} \text{cm}^{-3}$, $T_{kin} = 20_{-7}^{+20}$ K, and $\log_{10}(N_{CO}/dV) = 15.63_{-0.79}^{+1.14} \text{cm}^{-2} \text{km s}^{-1}$, where the errors are $3\sigma$. The warm component is in broad agreement with the previous single component fitting results.

### 6.4. The para-$H_2O(2_{02}-1_{11})$ emission line

Water is one of the most abundant molecules after $H_2$ and CO in the gaseous interstellar medium (ISM). It serves as a unique diagnostic for probing the physical conditions of the ISM in both local (e.g., González-Alfonso et al. 2014) and high-redshift galaxies (e.g., Omont et al. 2013; Riechers et al. 2013; Yang et al. 2016, 2020, and references therein). In particular, since the $E_{upper} > 200$ K levels of $H_2O$ are primarily excited

through absorption of far-infrared photons from dust emission in warm dense regions, it is a useful diagnostic of the far-infrared radiation field independent of gas conditions. In addition, sub-millimeter water lines are the second strongest molecular emitter after the CO lines (e.g., Yang et al. 2013). The detection of the para-$H_2O(2_{02}-1_{11})$ emission line ($v_{rest} = 987.9268$ GHz) in HerBS-89a is an illustration thereof: the line width is broad ($\Delta V \sim 1100 \text{ km s}^{-1}$) and comparable to those of the mid/high-$J$ CO lines, consistent with the finding that the $J=2 H_2O$ lines are spatially co-located with the CO lines (e.g., Omont et al. 2013; Yang et al. 2016). The observed line flux is bright, with $I_{H_2O} = 1.59 \pm 0.37 \text{ Jy km s}^{-1}$, and comparable to the line flux of the $^{12}$CO(9–8) emission line (Table 2). The CO and $H_2O$ lines have morphologies and velocity fields that are broadly comparable (see Figs. 4 and 7) indicating that they both trace the same warm high-density gas in HerBS-89a. Using Eqs. (6) and (7), we derive an intrinsic line luminosity (corrected for lensing magnification) of $L_{H_2O(2_{02}-1_{11})} = (6.4 \pm 1.5) \times 10^7 L_{\odot}$, yielding a ratio $L_{H_2O(2_{02}-1_{11})}/L_{IR} = (1.12 \pm 0.36) \times 10^{-5}$, in terms of the 8 to 1000 $\mu$m $L_{IR}$ (Table 7). HerBS-89a lies on the correlation between $L_{IR}$ and $L_{H_2O}$, that is $L_{H_2O(2_{02}-1_{11})} \propto L_{IR}^{1.06 \pm 0.19}$ derived by Yang et al. (2016), as shown in Fig. 16. This result is in line with the fact that far-infrared pumping is the dominant mechanism for the excitation of the submillimeter $H_2O$ emission that originates in very dense, heavily obscured star-formation dominated regions (see, e.g., González-Alfonso et al. 2014; Yang et al. 2016). The steeper-than-linear growth is likely the result of an increase of the optical depth in the dust continuum (at $\sim 100 \mu$m) with increasing infrared luminosity. Finally,
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Large Velocity Gradient (LVG) model fit to the ¹²CO SLED of HerBS-89a (after correcting for lens magnification). Top panel: data are shown as black dots and the solid curve shows the best fit from the single component (MCMC) model corresponding to the maximum posterior possibility (see Sect. 6.3 for details). Bottom panel: two-components fit, including an additional cold and extended component (darker, dot-dashed line).</figcaption></figure>

the ratio of $L_{H_2O(2\text{O}_2-1\text{I}_1)}/L_{\text{IR}}$ also indicates that an AGN component, if present in HerBS-89a, is not contributing much to the infrared luminosity, in contrast to a system like Mrk 231 or APM 08279+5455 where the AGN is dominant (and references therein Yang et al. 2016).

## 7. Molecular gas kinematics

In an ideal situation, analysis of spectral line data cubes like those presented here for HerBS-89a can enable detailed analyses of source-plane kinematics. The lens model derived in Sect. 4 can in principle be applied to the emission in a particular line across many independent velocity channels, allowing them to be reconstructed in the source plane, used to compute moment maps, and fits to dynamical models (e.g., Sharon et al. 2019). In practice, although the quality of our NOEMA observations is excellent, the S/N and the number of resolution elements (i.e., synthesized beams) across the observed size of HerBS-89a are not high enough to fully achieve this goal. Several factors are at

<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Correlation between $L_{\text{IR}}$ and $L_{H_2O(2\text{O}_2-1\text{I}_1)}$ for local (black symbols) and high-redshift (blue symbols) galaxies. HerBS-89a is labeled with a red star and is located along the tight relationship between the infrared and water luminosity (shown as a light blue line). Figure is adapted from Yang et al. (2016) and references therein.</figcaption></figure>

play here. First, as noted in Sect. 5 above, the combination of S/N and effective resolution limits us to the use of parametric (Sérsic profile) sources in reconstructing individual channels. Second, reconstruction of a given source-plane channel map tends to (a) favor sources that are very (sometimes implausibly) compact as a result of the use of a constrained emission profile, and (b) couple with the noise fluctuations in that channel in such a way that introduces uncertainties on the centroids of the reconstructed sources. Finally, comparison of a delensed source-plane reconstruction with an observed channel map necessarily involves the convolution of the former to the (coarser) resolution of the latter, meaning that many possible source-plane reconstructions are formally compatible with the observed data.

To provide a transparent (and possibly cautionary) picture of what properties can be derived from the current data, we show in Figs. 17 and 18 (as well as in Figs. A.3 and A.4) some of the results of our lens modeling efforts. Our main focus is the ¹²CO(9–8) line, which has higher S/N than the $H_2O(2\text{O}_2-1\text{I}_1)$ line and is therefore more promising as a kinematic probe. In order to improve the S/N, we rebin the 40 km s⁻¹ data to a resolution of 120 km s⁻¹, yielding a total of seven independent velocity channels that are suitable for modeling, although they do not span the line’s full velocity width. Each of these channels is independently reconstructed with a Sérsic profile in the source plane; these reconstructions are shown in the rightmost column of Fig. A.3, and are represented with the colored ellipses centered on dots in the left panel of Fig. 17 superposed on the dust continuum Sérsic profile in the source plane. The modest (generally <2.5σ in the vicinity of the line emission) level of the data–model residuals in the third column of Fig. A.3 indicates that these reconstructions are formally acceptable, although the five channels blue-ward of 168 km s⁻¹ have scale radii that are notably smaller than those of the two most highly redshifted channels. The centroids of all seven channels fall within the envelope defined by the reconstructed dust continuum emission, and in turn by the reconstruction of the ¹²CO(9–8) zeroth moment map (second row of Fig. 12), giving confidence that the overall spatial extent of the warmest and densest material in HerBS-89a is understood.

The center panel of Fig. 17 presents a source-plane zeroth moment map constructed by summing the source-plane
<figure><img src="image_18.png" /><figcaption>Fig. 17. Source plane reconstruction split by spectral channels. Left panel: positions (colored thick dots) and shapes (ellipses) of the reconstructed profiles of the individual $^{12}$CO(9–8) channels, over-plotted on the map of the dust continuum emission in the source plane. The ellipses are plotted at 1/20th of the peak flux of each reconstructed channel map. Central panel: zeroth moment map of the $^{12}$CO(9–8) emission line derived from the reconstruction of the seven individual channels shown in the left panel. The contours show the dust continuum emission, starting at 1/20 of the peak flux in steps of 1/20th (see also Fig. A.3). Right panel: zeroth moment map obtained splitting the $^{12}$CO(9–8) in two halves, corresponding to the velocity ranges from −312 to 168 km s$^{-1}$ and from 168 to 528 km s$^{-1}$, and reconstructing them separately. In all panels, the cross at the lower right corner marks the phase center of the NOEMA observations.</figcaption></figure>

reconstructions of the individual channel maps. This map seems to indicate an arc-like shaped morphology and the possibility of two separate components. However, we have here entered a regime in which the three limiting factors noted above are highly salient: because of the compactness of reconstructed Sérsic-profile sources in most of the velocity channels, the moment map derived from them suggests that we have reliable information about the internal source structure on much finer scales than is actually the case. It is also worth noting that when the moment map derived from the reconstructed channel maps is lensed forward into the image plane, it provides a poorer match to the observed $^{12}$CO(9–8) moment map than does the single Sérsic-component reconstruction in the second row of Fig. 12. This apparent paradox occurs because small noise fluctuations in the observed channel maps have unavoidably affected the source centroids, shapes, and orientations of the corresponding source-plane reconstructions. Similar inconsistencies have been seen by other authors who have separately delensed both channel maps and integrated moment maps (e.g., Dong et al. 2019). Use of prior information on the expected relationships between emitting regions in adjacent velocity channels (e.g., Rizzo et al. 2018; Young et al., in prep.) would be one way to mitigate this tendency, but we have not used such an approach here.

Figure 17 provides a further illustration of the degree to which small noise fluctuations in the observed channel maps can perturb the centroids of the source-plane reconstructions. Initially inspired by the apparent bent morphology (and nominal associated velocity gradients; see below) in the central panel of Fig. 17, we constructed new, roughly half-line moment maps from the bluest four and reddest three of our seven 120 km s$^{-1}$ channels, and reconstructed these using the lens model optimized in Sect. 5. The sum of these two reconstructions (i.e., another zeroth moment map) is presented in the right panel of Fig. 17. The apparent double morphology might at first right suggest a pre-coalescence merging system; however, the ‘doubleness’ of the morphology stems at least in part from the general compactness of our Sérsic profile reconstructions. Moreover, notable aspects of the internal source structure (e.g., the distinct southeastern offset of the reconstructed −312 to −192 km s$^{-1}$ channel map) in the center panel do not translate to

the reconstructions used to derive the right panel. This level of inconsistency raises the question of whether detailed information on the intrinsic or the internal source structure (as distinct from the overall extent of emission) can be robustly recovered from this dataset, and suggests that the apparent small-scale structure in the resulting zeroth-moment maps has to be interpreted with caution.

With the above discussion as preamble, we turn now to Fig. 18, which shows source-plane velocity fields derived from the seven-channel (120 km s$^{-1}$ resolution) reconstructions of Fig. A.3 in two different ways. The left panel presents an intensity-weighted first moment map calculated in the usual way. We find that the dominant velocity gradient is in an east-west direction, roughly consistent with the east-west elongation of the reconstructed dust continuum emission. The comparison of the center and right panels in Fig. 17 discussed above suggests that the perturbation from an “exactly” aligned east-west velocity gradient due to the most highly blue-shifted emission may not be robust.

Assuming a single rotationally supported disk, the radial profile of the dynamical mass of the system is given by

$$M_{\text{dyn}}(R) = f \cdot \frac{R V_{\text{obs}}^2(R)}{G \sin^2 i}, \quad (15)$$

where i is the inclination angle of the rotating disk with respect to the line of sight, and $V_{\text{obs}}$ is the observed velocity as a function of radius. Although the inclination could be derived from the ellipticity of the light profile, as $i = \arccos(1 - e)$, since the velocity gradient does not fully correspond to the continuum Sérsic profile in the source plane, we opt to keep the $\sin(i)$ term explicit in our computations. The constant f is a dimensionless scale factor that depends on the structure of the galaxy and has value $f \approx 1.0$ for a thin disk embedded in a massive spheroid (Lequeux 1983).

Along the major axis of the continuum profile (red line in the top-left panel of Fig. 18), we extract the velocity curve $V(R)$ shown in the middle-left panel. The corresponding bottom-left panel presents the radial dynamical mass profile, which leads to an estimate of $M_{\text{dyn}} \sin^2 i \approx 5.1 \times 10^{10} M_{\odot}$ within $R = 1.25$ kpc.

In contrast, the top-right panel of Fig. 18 presents a first moment map in which masking has been used to highlight the
<figure><img src="image_19.png" /><figcaption>**Fig. 18.** Kinematic and dynamical analysis of the $^{12}$CO(9–8) emission line in the source plane of HerBS-89a. **Top panels:** left: first moment map, compared to the Sérsic profile of the dust continuum emission. The contours start from 1/20 of the peak intensity with steps of 1/20th. The red line marks the continuum major axis. The reference position is defined where the rest-velocity ($V = 0$ km s$^{-1}$) and is denoted with a star. From this position we calculate the radial distance. The cross at the lower right corner marks the phase center of the NOEMA observations. **Right:** based on Figs. A.3 and 17, two possible components are identified; the first moment maps of the two components are here shown, after shifting the velocity scale by ±289 km s$^{-1}$ for clarity. The axes of the two components are shown as green and orange lines. The dashed line marks the contour of the velocity field shown in the top-left panel. **Middle panels:** velocity curves extracted along the axes marked in the **top panels**. Colors match those of the respective axes, with magenta points extracted from the top left moment map and the green and orange points extracted from the top right moment map. In the two right-most panels, the data are shown both before (small symbols) and after (large symbols) shifting by ±289 km s$^{-1}$. The amplitude of the error bars is shown in the left-most panel. **Bottom panels:** dynamical mass profile as a function of radius, as derived from the velocity curves.</figcaption></figure>

two regions of strongest integrated emission (compare also to Figs. 17 and A.3). These two regions velocity scales have been artificially offset by ±289 km s⁻¹, corresponding to the separation between the two peaks in the integrated CO spectrum (see Fig. 5) in order to better present them as “candidate” rotating disks. Velocity curves extracted along the axes of these two possible structures (green and orange lines) are shown in the middle row of Fig. 18, both before (small symbols) and after (large symbols) applying the ±289 km s⁻¹ shift.

The individual dynamical masses of these components, assuming rotational support as above, would be $M_{dyn} \sin^2 i \approx 2.6 \times 10^{10} M_{\odot}$ within $R = 0.85$ kpc for the eastern component and $M_{dyn} \sin^2 i \approx 2.2 \times 10^{10} M_{\odot}$ within $R = 0.40$ kpc for the western component, again under the assumption that $f = 1.0$ in Eq. (15). The bottom panels of Fig. 18 show the dynamical mass radial distributions, limited to the extension of the first-moment maps.

Assuming a baryon fraction $f_{baryon}$ (e.g., Tacconi et al. 2020), and supposing that at high redshift the gas content of a galaxy
is dominated by its molecular component and the atomic gas mass $M_{\text{HI}}$ can be neglected (see, e.g., Tacconi et al. 2018), the dynamical mass could yield an estimate of the stellar mass $M^* = f_{\text{baryon}}M_{\text{dyn}}-M_{\text{mol}}-M_{\text{HI}}$ and of the molecular gas fraction $f_{\text{mol}} = M_{\text{mol}}/(f_{\text{baryon}}M_{\text{dyn}}-M_{\text{mol}})$, thereby solving the question of whether or not the galaxy lies on the MS. However, the dynamical mass profiles derived for HerBS-89a do not extend to large enough radii to allow for a determination of the total mass.

For the sake of completeness, we also performed the reconstructions of multiple velocity channels across the H₂O(2₀₂−1₁₁) emission line, adopting the same approach as above, using a coarser velocity re-binning to boost the S/Ns (see Appendix A). As already suggested by the slightly different observed and reconstructed source-plane zeroth moment maps (see Fig. 12), the individual ¹²CO(9–8) and H₂O(2₀₂−1₁₁) velocity channels reconstruct to slightly different locations in the source plane (see Figs. A.3 and A.4). However, a first moment map derived from the reconstructed H₂O(2₀₂−1₁₁) channel maps shows a velocity gradient similar to the one shown in the top left panel of Fig. A.3 for the ¹²CO(9–8) emission line, albeit at a slightly different position angle, suggesting once more that this aspect of the system’s kinematics has been robustly recovered. Due to the higher S/N of the ¹²CO(9–8) emission line, we have limited the kinematic modeling to this line, rather than adventuring further in the realm of water.

To summarize, given the uncertainties discussed above, ambiguity remains as to whether the broad, double-peaked spectral lines in HerBS-89a reflect a single rotating disk or a pair of galaxies on their way to merging, as in the case of previous studies of high-z lensed starburst galaxies (e.g., Genzel et al. 2003; Ivison et al. 2010; Sharon et al. 2015). On balance, we favor the latter scenario, which would provide a natural explanation of the very large widths of the ¹²CO emission lines in this system. Observations with higher S/N and angular resolution that can support a more detailed lens model will be needed to resolve the one vs. two sources question in a definitive way for this system, shed further light on the properties of this galaxy, unveil the finer details of its structure and kinematics, and fully recover the secrets of its starburst nature.

# 8. Molecular lines other than CO and water

In the following subsections, we present and discuss the various molecular absorption and emission lines other than CO and H₂O that have been detected in the spectrum of HerBS-89a. The individual molecular lines are described starting with the very dense molecular gas tracers, namely: HCN and/or NH, and NH₂ (Sect. 8.1), followed by the two molecular ions, OH⁺ and CH⁺, which probe the low density molecular gas with low H₂ fractions (Sect. 8.2).

## 8.1. Very dense gas tracers: NH₂, and HCN and/or NH

NH₂. Amidogen (NH₂) is a molecule that is an important reactant intermediate in the production and destruction of N-bearing molecules, in particular as a photodissociation product of ammonia (NH₃). It is an asymmetrical molecule with a ²B₁ ground electronic state characterized by a complex rotational spectrum (similar to H₂O⁺). It was first detected in the interstellar medium (in SgrB2) by van Dishoeck et al. (1993) and Goicoechea et al. (2004). NH₂ was later detected in the Herschel SPIRE-FTS spectra of Mrk 231 and Arp 220 (González-Alfonso et al. 2018; Rangwala et al. 2011), where the NH₂(2₀₂−1₁₁)(5/2−3/2) emis-

sion line is seen next to the OH⁺(1₀−0₁) absorption line, and, possibly, the NH₂(2₀₂−1₁₁)(3/2−1/2) emission line as well.

In the case of HerBS-89a, two emission lines of NH₂ are detected. The first is the NH₂(2₀₂−1₁₁)(5/2−3/2) emission line, which overlaps with the OH⁺(1₀−0₁) absorption line, such that only the red-shifted part of the NH₂ emission line is apparent in the spectrum. However, the nearby NH₂(2₀₂−1₁₁)(3/2−1/2) remains undetected, which is in line with the expected ratio of ≈2 between the (5/3−3/2) and (3/2−1/2) transitions. The second detected transition NH₂(2₀₂−2₁₁)(5/2−3/2) is weaker but uncontaminated by nearby lines.

As the data are of low S/N, a detailed study of NH₂, first reported here in a high-z galaxy, remains out of reach. However, from the velocity-integrated emission map of the high-angular resolution data of the NH₂(2₀₂−2₁₁)(5/2−3/2) emission line (Fig. 4), it is clear that the NH₂ emission closely resembles the distributions of the ¹²CO(9–8) and water emission lines, with the southern peak being stronger, in accordance with the fact that this molecule also traces the high-density gas. Dedicated future observations can further explore the properties of this molecule in HerBS-89a or in other high-z galaxies to study the properties of their dense gas reservoirs.

**HCN(11–10) and/or NH.** The absorption feature observed at 246.72 GHz in HerBS-89a may be due to redshifted HCN(11–10) and/or NH(1₂–0₁). With the current data, there is no way to distinguish between these molecular species, hence we briefly describe each molecule separately. The absorption feature is relatively weak and detected at a low S/Ns. As shown in Fig. 4, it is spatially well centered on both dust continuum emission peaks. This is the first detection of a high-J HCN line and/or NH in a high-z galaxy. Clearly, detecting other transitions of HCN and NH in HerBS-89a will help to determine the relative contributions of each of these species to the observed 246.72 GHz absorption feature.

– HCN is known to be a reliable tracer of dense molecular gas ($n_{\text{H}_2} \gtrsim 3 \times 10^4 \text{ cm}^{-3}$), with critical densities 100 to 1000 times higher than of CO. Although HCN emission lines are typically one to two orders of magnitude fainter than CO emission lines, HCN has been observed in its low-J transitions in local starburst galaxies, revealing a tight correlation between the HCN luminosity and the star formation rate (e.g., Gao & Solomon 2004; Zhang et al. 2014). Emission lines of low-J HCN have also been detected in a few high-z galaxies (Weiß et al. 2007; Riechers et al. 2010; Oteo et al. 2017; Cañamas et al. 2021, and references therein), and HCN is seen in the stacked rest-frame 220–770 GHz spectrum of 22 high-z DSFGs selected from the SPT survey (Spilker et al. 2014).

Between J = 5–4 and 11–10, HCN flips from emission to absorption as shown in the case of Arp 220 where Rangwala et al. (2011) detect the J = 12–11 to 17–16 transitions in absorption<fnref n="6" />. These transitions of HCN have critical densities around 10⁹–10¹⁰ cm⁻³ and probe very dense ($n_{\text{H}_2} \approx 3 \times 10^6 \text{ cm}^{-3}$) gas. The high-J lines of HCN are thought to be populated by radiative pumping of infrared photons (rather than by collisional excitation) in an intense high-temperature (>350 K) radiation field Rangwala et al. (2011).

In HerBS-89a, the optical depth of the 246.72 GHz line is estimated to be $\tau(\nu_0) \sim 0.24 \pm 0.12$ yielding, in the case of HCN, a column density of $(0.38 \pm 0.1) \times 10^{14} \text{ cm}^{-2}$ (Table 3). Future observations of both high-J HCN lines (in absorption) and low-J lines (in emission) in HerBS-89a would allow us to confirm if

<footnote n="6">⁶ The HCN(11–10) transition was not discussed in Rangwala et al. (2011).</footnote>
this absorption line is indeed due (or partly due) to HCN, and, once verified, allow us to probe and constrain the properties of the reservoir of extremely dense gas in this high-$z$ starburst galaxy.

– NH (imidogen) is a molecular radical that is linked (as for NH₂) to the formation and destruction of NH₃. It was first detected in the diffuse interstellar medium (Meyer & Roth 1991) and later in the dense molecular gas of SgrB2 (Cernicharo et al. 2000; Goicoechea et al. 2004). The Herschel SPIRE-FTS spectra of Mrk 231 and Arp 220 (González-Alfonso et al. 2018; Rangwala et al. 2011) show the presence of NH with various transitions detected. NH has a spectrum similar to that of OH⁺, since these molecules are isoelectronic. Many of the transitions of NH overlap with OH⁺ transitions, in particular NH(1₂–1₀) with OH⁺(1₀–0₁), with similar column densities derived for OH⁺ and NH in Mrk 231 (González-Alfonso et al. 2018). In the case of HerBS-89a, attributing all of the 246.72 GHz absorption line to NH implies a column density of (6.7 ± 1.9) × 10¹⁴ cm⁻², which is also comparable to the column density of OH⁺ (Table 3). Together with NH₂, the likely detection of NH in HerBS-89a opens up the possibility of exploring the chemistry of N-bearing molecules in a high-$z$ galaxy by measuring other transitions of NH and searching for NH₃, and exploring whether (as is the case in SgrB2) the derived column densities are compatible with grain-surface chemistry and sputtering by shocks (see, e.g., Goicoechea et al. 2004).

## 8.2. Low density gas tracers: OH⁺ and CH⁺

OH⁺. The molecular ion OH⁺ is a reliable and powerful tracer of molecular outflows in DSFGs as has been demonstrated with the Herschel SPIRE-FTS studies in Arp 220 (Rangwala et al. 2011) and Mrk 231 (González-Alfonso et al. 2018), and in high-$z$ galaxies with NOEMA in HLFS 3 (Riechers et al. 2013) and ALMA (Indriolo et al. 2018, who reported OH⁺(1₁–0₁) absorption line in two lensed sources, i.e., the Cosmic Eyelash and SDP.17). As discussed in Indriolo et al. (2018), OH⁺ absorption is thought to arise in the extended cool, diffuse, low H₂ fraction gas that surrounds galaxies. OH⁺ has been found to be useful in constraining the cosmic-ray ionization rate of atomic hydrogen, ζ$_{H_2}$, particularly in gas with a low H₂ fraction. The three ground state lines of OH⁺ (all are split into hyperfine components, which are blended due to the broad line widths measured in galaxies) are generally observed in absorption. In a few cases, they can be seen in emission as a red wing, as observed in Mrk 231 or Arp 220, or as recently reported for the OH⁺(1₁–0₁) line in the $z = 6.03$ quasar SDSS231038+1855 by Li et al. (2020) and the hot dust-obscured galaxy W0410–0913 at $z = 3.63$ by Stanley et al. (2021), where they can provide a measurement of the local density (e.g., Rangwala et al. 2011).

In HerBS-89a, all three OH⁺ ground-state $1_J-0_I$ lines ($J = 0, 1, 2$) are seen in absorption (Fig. 5 middle panel). The ($1_I-0_I$) and ($1_2-0_I$) transitions are the strongest, while the ($1_0-0_I$) transition is nearly five times weaker (Table 2). None of the OH⁺ lines in HerBS-89a show any indication of emission above the continuum as in Mrk 231 or Arp 220 (González-Alfonso et al. 2018; Rangwala et al. 2011). All three lines have comparable line widths with $\Delta V \approx 500 \text{ km s}^{-1}$ and the inferred OH⁺ column density is $N_{\text{OH}^+} \sim 10^{15} \text{ cm}^{-2}$ (Table 3), comparable to the value derived for the Cosmic Eyelash (Indriolo et al. 2018). The OH⁺ spectra show clear indications of red-shifted gas, with all three transitions peaking in opacity at $\approx+100 \text{ km s}^{-1}$ relative to the systemic velocity, including the ($1_0-0_I$) transition after the contamination of the NH₂ emission line is accounted for (see Sect. 3.1.2).

CH⁺. The methylidyne cation, CH⁺, has been shown to be a sensitive tracer of feedback mechanisms in high-$z$ DSFGs (Falgarone et al. 2017), who report the detection of CH⁺(1–0) in six lensed galaxies at redshifts $z \approx 2.5$, including the Cosmic Eyelash and SDP.17⁷. The CH⁺(1–0) line is found both in emission and absorption, providing critical information on the mechanisms of mechanical energy dissipation and/or strong UV irradiation in these high-$z$ galaxies. Due to its high critical density for excitation ($\sim 10^7 \text{ cm}^{-3}$ for the $J = 1-0$ transition), the CH⁺ line is seen in emission in dense ($>10^5 \text{ cm}^{-3}$) gas, probing shocked regions powered by galactic winds, and in absorption, in large ($>10 \text{ kpc}$) reservoirs of turbulent, cold ($\sim 100 \text{ K}$), low-density ($\sim 100 \text{ cm}^{-3}$) gas (see, e.g., Falgarone et al. 2017).

In HerBS-89a, the CH⁺(1–0) is detected in absorption with a width of $\approx 400 \text{ km s}^{-1}$, comparable to that of the OH⁺ lines, and is also red-shifted to $\approx+100 \text{ km s}^{-1}$ (Fig. 5). The inferred CH⁺ column density is estimated to be $N_{\text{CH}^+} \sim 10^{14} \text{ cm}^{-2}$ (Sect. 8), similar to what has been derived in the case of the Cosmic Eyelash (Falgarone et al. 2017).

The low S/N, tentative detection of CH⁺ emission could indicate that in HerBS-89a there is also dense gas ($>10^5 \text{ cm}^{-3}$) with a significant velocity dispersion (the line has a full width at zero intensity of $\gtrsim 1000 \text{ km s}^{-1}$ and is centered at the systemic velocity), as was reported for three of the sources studied by Falgarone et al. (2017). However, higher quality data are required to probe this dense gas component in greater detail.

## 9. Inflowing gas

In HerBS-89a, the three OH⁺ ground-state $1_J-0_I$ lines ($J = 0, 1, 2$) and the CH⁺(1–0) line are all seen in absorption and red-shifted with respect to the systemic velocity (defined by the CO emission lines) by $\approx+100 \text{ km s}^{-1}$. As shown in the velocity maps of these absorption lines, the red-shifted gas is distributed over the southern and northern dust continuum peaks (in the image plane), uniformly covering the distribution of the blue- and red-shifted gas traced in the $^{12}\text{CO}(9-8)$ and water emission lines (Fig. 4). The kinematics of the absorbing gas in HerBS-89a indicate that this red-shifted gas is therefore not kinematically related to the dense molecular gas. We argue here that we are in fact tracing, in both OH⁺ and CH⁺, low density gas that is flowing into the central regions of HerBS-89a.

The number of high-$z$ sources for which measurements of absorption by molecular ions such as CH⁺ or OH⁺ are available is still low; however, it is noteworthy that the majority of the sources for which observations of these species have been made do show blue-shifted absorption lines, indicating vigorous outflow activity. Out of the five high-$z$ lensed galaxies where CH⁺(1–0) has been found in absorption, only one source<fnref n="8" /> shows a red-shifted line, namely the Cosmic Eyelash (Falgarone et al. 2017); for this source, the CH⁺(1–0) absorption line is shifted by $\Delta V \sim +150 \text{ km s}^{-1}$ with respect to the systemic velocity as defined by the CO emission lines. A similar velocity shift in the Cosmic Eyelash is also observed in the OH⁺(1₁–0₁) and H₂O⁺(1₁₁–0₀₀) absorption lines

<footnote n="7">⁷ The first detection of CH⁺ was made in the far-infrared spectrum of the planetary nebula NGC 7027 through rotational lines from $J = 2-1$ to $6-5$ seen in emission Cernicharo et al. (1997).</footnote>

<footnote n="8">⁸ In the case of the galaxy G09v1.40 reported in Falgarone et al. (2017), adopting the precise redshift of $z = 2.0924 \pm 0.0001$ (Yang et al. 2017) shifts the CH⁺(1–0) absorption line to $\sim -200 \text{ km s}^{-1}$ with respect to the systemic velocity of the source as defined by the CO emission lines; the CH⁺(1–0) absorption line is therefore blue-shifted in this system – see also Butler et al. (in prep.).</footnote>
(Indriolo et al. 2018). Recently, Butler et al. (in prep.) report $OH^+(1_1-0_1)$ absorption lines toward ten gravitationally lensed Herschel-selected high-z galaxies; for all of them, the line is invariably seen in absorption and blue-shifted with respect to the systemic velocity, with the exception of one source where the absorption line peaks at the systemic velocity. Together, the currently available observations of $OH^+$ (and $CH^+$) have established that these molecular ions are robust tracers of molecular gas outflows in high-z galaxies, confirming the results of observations of $OH^+$ or OH in local ultra-luminous galaxies (e.g., Lu et al. 2017; González-Alfonso et al. 2017, and references therein). In addition, recent studies of absorption lines in high-z lensed DSFGs in the OH 119 µm doublet, in SPT-2319-55 at $z = 5.29$ (Spilker et al. 2018) and the quasar ULASJ131911+095051 at $z = 6.13$ (Herrera-Camus et al. 2020), and in the $H_2O(3_{30}-3_{21})$ and $(4_{23}-4_{14})$ transitions in SPT-0346-52 at $z = 5.656$ (Jones et al. 2019), have revealed blue-shifted lines tracing powerful outflows (with rates of ~200–500 $M_{\odot}$ yr$^{-1}$).

In contrast to the above trends, the detection of red-shifted molecular absorption lines in high-z systems seems much less common. Only two cases are known, namely the Cosmic Eye-lash (Falgarone et al. 2017; Indriolo et al. 2018) and HerBS-89a (here described), out of the total of 18 high-z sources observed to date for which $CH^+(1-0)$ and/or $OH^+(1_1-0_1)$ (or other tracers such as OH, $H_2O$ or $H_2O^+)$ have been seen in absorption.

Direct observational evidence of inflow is indeed scarce. Signatures of infalling gas are notoriously difficult to observe: the accreting material may have low metallicity, low density, and/or a small covering factor (less than 10% at $z \sim 1.5$; see, e.g., Fumagalli et al. 2011), and it might be ionized or obscured by outflows. The need for a good alignment of an accretion structure with a background continuum source along the line of sight lowers the chance of detection even further.

The findings here described for molecular absorption lines in high-z galaxies are in line with the scarcity of infall evidence gathered from previous studies. Using optical spectroscopy, Rubin et al. (2012) report evidence of cool, metal-enriched gas accretion onto galaxies at $z \sim 0.5$. Wiseman et al. (2017) found a rare case of inflow of metal-poor gas from the intergalactic medium onto a $z = 3.35$ galaxy. More recently, Ao et al. (2020) report evidence of infalling gas in a Lyman-$\alpha$ blob at $z = 2.3$, and Daddi et al. (2020) presented observations revealing cold gas filaments accreting toward the center of a massive group galaxy at $z = 2.91$. In the local Universe, a few examples of inflow activity have been found through Herschel spectroscopic observations of compact luminous infrared galaxies. These include NGC 4418 and Zw 049.057, where reversed P-Cygni profiles, a clear signature of infall, were revealed in the [OI] 63 µm fine-structure lines (González-Alfonso et al. 2012; Falstad et al. 2015); Arp 299 where the absorption of the ground state of OH at 119 µm is found red-shifted by ~+175 km s$^{-1}$, suggesting a low excitation inflow in this source (Falstad et al. 2017); Circinus, which shows an unambiguous OH 119 µm inverted P-Cygni profile (Stone et al. 2016); and IRAS 10091+4704, the starburst-dominated system at $z = 0.24$, were a clear inverted P-Cygni OH profile indicates the presence of a fast inflow of molecular gas with a rate of ~100 $M_{\odot}$ yr$^{-1}$ (Herrera-Camus et al. 2020).

Theoretically, gas flowing into less massive ($M \le 10^{12} M_{\odot}$) galaxies should be dynamically and thermally cold, while more massive halos receive most of their baryons as hotter ($T > 10^5$ K) gas (e.g., Dekel & Woo 2003; Kereš et al. 2005). Thus cold, dense, metal-poor circum-galactic gas is often interpreted as direct evidence of accretion. Metal-poor Lyman-limit systems have been used as tracers of accretion

<figure><img src="image_20.png" /><figcaption>Fig. 19. Schematic of the simplified geometry of the inflow.</figcaption></figure>

in observations (e.g., Lehner et al. 2013) and simulations (e.g., Fumagalli et al. 2011).

Inflow and outflow rates of galaxies are constrained by cosmological simulations (e.g., Kereš et al. 2005) and analytic models based on the mass assembly history of dark matter halos (Bouché et al. 2010; Lilly et al. 2013). Erb (2008) reproduces the gas mass fraction and metallicity distribution of galaxies with a simple chemical model involving gas inflow and outflow. Inverting their reasoning, Yabe et al. (2015) constrain gas inflow and outflow rates based on the observed stellar masses, gas mass fractions, and metallicities of star-forming galaxies at $z \sim 1.4$.

Taken together, the above evidence strongly indicates that galaxies possess large reservoirs of circum-galactic gas eligible for accretion (Tumlinson et al. 2017). Nevertheless, evidence for fuel does not necessarily imply the presence of fueling.

We can describe the gas flow as a cylindrical channel that brings the gas onto the galaxy (see Fig. 19 for a schematic representation). Given the distance $d_{\text{inflow}}$ of the infalling clouds from the galaxy (i.e., given the line-of-sight distance that the gas has yet to travel), and assuming that the accreting gas travels at a constant speed $V_{\text{inflow}}$, the mass inflow rate can be computed as follows:

$$ \frac{dM_X}{dt} = N_1(X)m_X \frac{A_{\text{obs}}}{\sin \phi} \frac{V_{\text{obs}}}{d_{\text{inflow}} \cos \phi} \quad (16) $$

Here $X$ denotes either the $OH^+$ or $CH^+$ molecular ion; $N_1$ is the column density of the inflowing gas; $m_X$ is the mass of the molecular ion; $A_{\text{obs}}$ is the de-magnified projected area covered by the inflowing gas in the sky; $\phi$ is the angle between the line of sight and the inflow direction; and $V_{\text{obs}}$ is the measured red-shifted velocity of the line. With these definitions, the term $V_{\text{obs}}/\cos \phi$ represents the de-projected inflow velocity $V_{\text{inflow}}$, and the cross-sectional area of the cylinder perpendicular to the flow is given by $A_{\text{inflow}} = A_{\text{obs}}/\sin \phi$.

The $N_1(\text{OH}^+)$ value measured from the high angular resolution spectra for $OH^+(1_1-0_1)$ and $(1_2-0_1)$ is $(9.6-11.2 \times 10^{14}) \text{cm}^{-2}$ (Table 3). The mass of the $OH^+$ molecular ion is $m_{\text{OH}^+} = 2.824 \times 10^{-23} \text{ g}$. The observed projected area can be computed from the Sérsic radius $R_s$ of the $OH^+$ source plane reconstruction ($0.21-0.23 \text{ kpc}$; Table 6) yielding $A_{\text{obs}} = 0.13-0.17 \text{ kpc}^2$. The observed velocity of the flow is 100 km s$^{-1}$. The distance that the inflowing gas travels from the observed position to the center of the galaxy is one of the most uncertain quantities in Eq. (16). Typical distances found in the literature are on the order of one to a few kpc. Using the positions of lines and continuum in the source plane (Table 6), we can only compute the distances projected on the sky as between ~0.1 and ~0.3 kpc, which are of the same order than the systematic positional uncertainty computed via Eq. (4) (~0.08 kpc). Assuming a geometry similar to that of $OH^+$, because no source plane reconstruction has been attempted for the low angular resolution
NOEMA data, a similar reasoning can be applied to the CH⁺ line. Based on these values and estimates and using Eq. (16), we find

$$ \frac{dM_{OH^+}}{dt} \sin \phi \cos \phi \sim 1.8-8.8 \times 10^{-6} M_{\odot} \text{yr}^{-1}, \quad (17) $$

and

$$ \frac{dM_{CH^+}}{dt} \sin \phi \cos \phi \sim 0.4-1.6 \times 10^{-6} M_{\odot} \text{yr}^{-1}, \quad (18) $$

modulo the angle $\phi$ between the flow direction and the line of sight, which is unknown.

In order to transform these results into the total inflow rate, we assume that the gas mass is dominated by hydrogen, and we adopt as abundances relative to hydrogen the quantities $f_{OH^+} = 5.7 \times 10^{-8}$ (computed as the Galactic average value from Table 5 of Indriolo et al. 2015) and $f_{CH^+} = 7.6 \times 10^{-9}$ (Godard et al. 2014; Indriolo et al. 2018). By dividing Eq. (16) by $f_X$ and multiplying by the ratio $m_H/m_X$ of the molecular ion's mass relative to H, we derive mass inflow rates of $\dot{M}_{\text{inflow}} \sin\phi \cos\phi \sim 3-14 M_{\odot} \text{yr}^{-1}$ from OH⁺, and $\sim 6-24 M_{\odot} \text{yr}^{-1}$ from CH⁺.

In the so-called “bathtub” or “gas regulator” model (Bouché et al. 2010; Lilly et al. 2013; Sánchez Almeida et al. 2014; Peng & Maiolino 2014; Dekel & Mandelker 2014; Somerville & Davé 2015, among many others), the evolution of the gas content of a galaxy can be formulated as follows:

$$ \dot{M}_{\text{gas}} = \Phi - (1 - R + (\lambda - \eta_{\text{rec}})) \epsilon M_{\text{gas}} + Y, \quad (19) $$

where $\Phi$ is the inflow rate; SFR = $\epsilon M_{\text{gas}}$, as per the Schmidt-Kennicutt star formation law (Schmidt 1959; Kennicutt 1998b); R is the fraction of mass of newly formed stars quickly returned to the ISM through stellar winds and supernovae; $\Psi = \lambda$SFR is the outflow rate, linked to the SFR by the mass loading factor $\lambda$; finally, $\eta_{\text{rec}}$ is the fraction of gas ejected by outflows falling back onto the galaxy. Additional factors, for example due to AGN feedback, possible mergers or other external causes, are here summarized by the term Y.

Different scenarios can be at play in the case of HerBS-89a. In the least likely scenario, we might be witnessing a flow channeling gas from the intergalactic medium onto the galaxy, described by the term $\Phi$ of Eq. (19). While such flow would in principle consist of pristine gas, the detected OH⁺ and CH⁺ transitions testify to the presence of enriched material. To produce these molecular ions, first of all, as the reaction rates go with $n^2$, pockets of dense gas should be present; second, the most important formation channel for OH⁺ (and CH⁺) requires high temperatures, which could be produced in an accretion shock in the infalling gas as it enters the halo of the galaxy, although this scenario remains speculative.

Yabe et al. (2015) use a simple analytic model of galaxy chemical evolution to constrain the inflow (and outflow) rates of three samples of star forming galaxies at $z = 0, 1.4$ and $2.2$ (Erb et al. 2006; Peeples & Shankar 2011, and their work). They show that the average gas inflow rate of the samples is 2 to 3x their SFRs. Turning into physical units, $\dot{M}_{\text{inflow}}$ is a few solar masses per year at $z = 0$ growing to $\sim 100 M_{\odot} \text{yr}^{-1}$ at $z = 2.2$. Model predictions reproduce these values for dark matter halos masses of $M_{\text{DH}} \sim 2 \times 10^{12} M_{\odot}$ (e.g., Bouché et al. 2010; Yabe et al. 2015). If the SFR is fed by the inflow, then it is naturally proportional to the inflow rate. The detailed process of gas accretion and its relation to the fuel of SFR is more complex and still poorly understood. The gas fraction and SFR efficiency (i.e., the reciprocal of $\tau_{\text{depl}}$) come into play. How they relate to

$\dot{M}_{\text{inflow}}$ is yet to be fully investigated by observations, because of the rare occurrence of direct inflow measurements.

In this frame, HerBS-89a would require a very high inflow rate to sustain its large measured star formation rate. Although the angle between the line of sight and the inflow is unknown, the values derived above seem too low to support this scenario.

Alternatively, the flow could be part of a galactic fountain, in which a fraction of the gas expelled by the starburst winds falls back onto the galaxy (the term $\eta \epsilon M_{\text{gas}}$ in Eq. (19)). In this case, an outflow should also be detected, as described by $-\lambda \epsilon M_{\text{gas}}$. No evidence of any outflow feature was found in the spectra of HerBS-89a, neither in the form of blue-shifted emission (e.g., inverted P-Cygni profile) nor through absorption components of the detected lines. Nevertheless, it is worth noting that the absorption lines of OH⁺ here reported are not isolated, and have emission or absorption lines of other molecular species on their blue sides, perhaps preventing the identification of a blue-shifted component. Moreover – despite the superb quality of the NOEMA spectra – the S/N reached in the spectra of the main emission lines is not high enough to detect any faint, very broad wings, that could indicate the presence of an outflow.

Finally, if HerBS-89a is in fact a merging pair or a multiple source (see Sect. 7), a third possibility is that the red-shifted OH⁺ and CH⁺ lines originate in a stream of matter between the two components. This case belongs to the term Y in Eq. (19), and would naturally explain the enriched composition of the flow. One possibility is that HerBS-89a is similar to a scaled-up version of the merger in the local Antennae galaxy (NGC 4038/9), in which an off-nuclear high concentration of dust and gas and a powerful starburst are found in the region between the two galaxies: a “Hyper-Antennae”.

With the data in hand, it is however not possible to distinguish between these different scenarios and draw a definitive conclusion. Also in this case, higher angular resolution and deeper sensitivity observations of HerBS-89a will be required to shed further light on the nature of the mass inflow observed in this system.

## 10. Summary and concluding remarks

We report new millimeter NOEMA and VLA observations of the lensed starburst galaxy HerBS-89a at $z = 2.95$ that are complemented with near-infrared and optical data (obtained with the HST and GTC) that reveal the foreground lensing galaxy. The ability to process large instantaneous bandwidths with NOEMA enables an in-depth exploration of a sizeable part of the 1 mm spectrum of HerBS-89a, covering in total 121 GHz in the rest-frame, leading to the detection of transitions of molecular ions and molecular species that have not been previously observed in a high-z galaxy. In addition to the strong emission lines of $^{12}$CO(9–8) and H$_2$O(2$_{02}$–1$_{11}$), we detect the three OH⁺ ground-state $1_{J-01}$ lines ($J = 0, 1, 2$) and the CH⁺(1–0) line, all four seen in absorption and red-shifted with respect to the systemic velocity of HerBS-89a, two transitions of NH$_2$ in emission, and HCN(11–10) and/or NH(1$_2-01$) lines, seen in absorption. Modeling the lensing of HerBS-89a, we reconstructed the dust continuum emission and the molecular emission lines in the source plane deriving the morphology and the kinematics of this high-z galaxy.

The main conclusions of this paper are as follows:

- The high-angular resolution (0.′3) images of HerBS-89a in dust continuum and the emission lines of $^{12}$CO(9–8) and H$_2$O(2$_{02}$–1$_{11}$) show two distinct components – one to the north and a brighter slightly more extended component to the
south, which are linked by an arc-like 1.′′0 diameter structure reminiscent of an Einstein ring. HerBS-89a is gravitationally amplified by a foreground massive galaxy at $z_{phot} = 0.9^{+0.3}_{-0.5}$, with an estimated deflector mass of $M_{lens} = 9.8^{+3.2}_{-4.9} \times 10^{10} M_{\odot}$. Using lens modeling, we have reconstructed the source-plane morphologies of HerBS-89a in dust continuum and molecular line emission. The magnification factor is estimated to be $\mu \sim 5$ for the dust continuum and $\mu \sim 4-5$ for the emission lines.

- HerBS-89a is a high-z dusty star forming galaxy. From the $^{12}$CO(1–0) measurements, we derive a mass for the molecular gas reservoir of $M_{mol} = (2.1 \pm 0.4) \times 10^{11} M_{\odot}$. modeling the available CO lines ($J = 1-0, 3-2, 5-4$ and $9-8$) with an LVG approach indicates that the molecular gas in HerBS-89a has an average density of $n_{H_2} \sim 500 \text{ cm}^{-3}$ and a kinetic temperature $T_{kin} \sim 200$ K. Using the fine sampling of the dust continuum flux densities observed from $250 \mu\text{m}$ to $3.4 \text{ mm}$ (or $63$ to $860 \mu\text{m}$ in the rest-frame), we derive an infrared luminosity of $L_{IR} = (4.6 \pm 0.4) \times 10^{12} L_{\odot}$, and a dust mass of $M_{dust} = (2.6 \pm 0.2) \times 10^9 M_{\odot}$, yielding a dust-to-gas ratio $\delta_{GDR} \approx 80$. The star formation rate is estimated to be SFR = $614 \pm 59 M_{\odot} \text{ yr}^{-1}$ and the depletion timescale $\tau_{depl} = (3.4 \pm 1.0) \times 10^8$ years.

- Lens modeling suggests that the warmest and densest material in HerBS-89a (i.e., as traced by dust, $^{12}$CO(9–8), and $H_2O(2\text{O}_2-1\text{I}_1)$ emission) is distributed across a region a few kpc across. We cannot distinguish definitively between a single-disk and a two-component description of this distribution on the basis of existing observations, but dynamical and molecular gas mass estimates are consistent in either scenario, and the two-component scenario in particular provides a natural explanation for the very large observed CO velocity widths.

- A series of molecular emission and absorption lines are detected in HerBS-89a that trace very dense gas (up to $n_{H_2} \sim 10^9-10^{10} \text{ cm}^{-3}$). $NH_2$ is seen in emission and detected in two transitions, while an absorption feature at $246.72 \text{ GHz}$ corresponds to HCN(11–10), $NH(1_2-1_0)$ or a combination of both. Future observations of additional transitions of these species (particularly HCN and NH) would be needed to further explore the densest regions of HerBS-89a.

- The three fundamental transitions of $OH^+$ and $CH^+(1-0)$ are all seen in absorption and all red-shifted by $\approx+100 \text{ km s}^{-1}$ with respect to the systemic velocity of HerBS-89a. The molecular ions $OH^+$ and $CH^+$ both trace low-density ($n_{H_2} \sim 100 \text{ cm}^{-3}$) gas. The fact that the velocity red-shifted absorption lines cover substantially the distribution of the CO and water emission lines in the image (and source) planes indicates that this low-density gas must be kinematically unrelated to the dense gas. We here argue that we are tracing a rare case of inflow of low-density gas, indicating that HerBS-89a is accreting matter from its surroundings. Based on the available data, we have estimated that the mass inflow rate is on the order of $\approx 10-20 M_{\odot} \text{ yr}^{-1}$. Various scenarios are possible to account for the observations including inflowing gas from the intergalactic medium, a galactic fountain or a stream of matter between the two merging components.

The results here described highlight the importance of processing large instantaneous bandwidth in the study of high-z galaxies, as such observations allow us to efficiently measure multiple lines and explore in an unbiased way the entire available frequency range for the presence of molecular species. In the case of HerBS-89a, we were able to probe this system from the very highest to the lowest molecular gas densities with high-quality

spectral data in the 830–1060 GHz rest-frame frequency range, which includes key molecular ionic transitions tracing feedback mechanisms. The findings here reported show that it is possible to unravel the complex physics at play within this starburst (likely to be a merger) system in the early Universe from the kinematics of its inner most regions to the low-density inflowing molecular gas. They provide the foundation for future follow-up observations to further study the nature of HerBS-89a and better constrain its morphology and kinematics, which will need higher-angular resolution and higher S/N data, as well as for future studies of larger samples of high-z galaxies, with properties akin to those of HerBS-89a, selected from on-going large redshift surveys.

**Acknowledgements.** This work is based on observations carried out under programs, E18AE (a Discretionary Directorial Time project) and W19DE, using the IRAM Northern Extended Millimeter Array (NOEMA), and, program VLA/20A-083 using the National Radio Astronomy Observatory's (NRAO) Karl G. Jansky Very Large Array (VLA). IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The authors are grateful to the IRAM Director, K.F. Schuster, for enabling the initial observations of this project. This paper makes also use of observations made with the Gran Telescopio Canarias (GTC), installed at the Spanish Observatorio del Roque de los Muchachos del Instituto de Astrofísica de Canarias, in the island of La Palma, under Director's Discretionary Time. The authors are grateful to the reference for constructive and insightful comments. We would like to thank J. Cernicharo, F. Cox, and E. Falgarone for useful discussions about some aspects of this study and E. Chapillon for her help with the JVLA data reduction. A.J.Y. and A.J.B. acknowledge enlightening discussions with C.R. Keeton about the lens modeling. The authors would also like to highlight interesting exchanges triggered by the thoughts of J. de Chabannes. This work benefited from the support of the project Z-GAL ANR-AAPG2019 of the French National Research Agency (ANR). A.J.B. and A.J.Y. acknowledge support from the National Science Foundation grant AST-1716585. B.M.J. acknowledges the support of a studentship grant from the UK Science and Technology Facilities Council (SFTC). L.D. acknowledges support from the ERC consolidator grant CosmicDust (PI: H. Gomez). T.B. acknowledges funding from NAOJ ALMA Scientific Research Grant 2018-09B and JSPS KAKENHI No. 17H06130. C.Y. acknowledges support from an ESO Fellowship. S.D. is supported by a UK STFC Rutherford Fellowship. H.D. acknowledges financial support from the Spanish Ministry of Science, Innovation and Universities (MICIU) under the 2014 Ramón y Cajal program RYC-2014-15686 and AYA2017-84061-P, the later one co-financed by FEDER (European Regional Development Funds). D.A.R. acknowledges support from the National Science Foundation under grant numbers AST-1614213 and AST-1910107 and from the Alexander von Humboldt Foundation through a Humboldt Research Fellowship for Experienced Researchers. S.J. acknowledges financial support from the Spanish Ministry of Science, Innovation and Universities (MICIU) under grant AYA2017-84061-P, co-financed by FEDER (European Regional Development Funds) and by the Spanish Space Research Program 'Participation in the NISP instrument and preparation for the science of EUCLID' (ESP2017-84272-C2-1-R).

## References

- Ao, Y., Zheng, Z., Henkel, C., et al. 2020, Nat. Astron., 4, 670

- Arnouts, S., Cristiani, S., Moscardini, L., et al. 1999, MNRAS, 310, 540

- Babbedge, T. S. R., Rowan-Robinson, M., Gonzalez-Solares, E., et al. 2004, MNRAS, 353, 654

- Bakx, T. J. L. C., Eales, S. A., Negrello, M., et al. 2018, MNRAS, 473, 1751

- Bakx, T. J. L. C., Eales, S. A., Negrello, M., et al. 2020, MNRAS, 494, 10

- Barvainis, R., Maloney, P., Antonucci, R., & Alloin, D. 1997, ApJ, 484, 695

- Battisti, A. J., da Cunha, E., Grasha, K., et al. 2019, ApJ, 882, 61

- Berta, S., Fritz, J., Franceschini, A., Bressan, A., & Pernechele, C. 2003, A&A, 403, 119

- Berta, S., Lutz, D., Santini, P., et al. 2013, A&A, 551, A100

- Berta, S., Lutz, D., Genzel, R., Förster-Schreiber, N. M., & Tacconi, L. J. 2016, A&A, 587, A73

- Bianchi, S. 2013, A&A, 552, A89

- Bianchi, S., Casasola, V., Baes, M., et al. 2019, A&A, 631, A102
Blain, A. W., Smail, I., Ivison, R. J., Kneib, J.-P., & Frayer, D. T. 2002, *Phys. Rep.*, **369**, 111

Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, *ARA&A*, **51**, 207

Bolzonella, M., Miralles, J. M., & Pelló, R. 2000, *A&A*, **363**, 476

Boselli, A., Eales, S., Cortese, L., et al. 2010, *PASP*, **122**, 261

Bothwell, M. S., Smail, I., Chapman, S. C., et al. 2013, *MNRAS*, **429**, 3047

Bouché, N., Dekel, A., Genzel, R., et al. 2010, *ApJ*, **718**, 1001

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

Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, *ApJ*, **686**, 1503

Bruzual, A. G. 1983, *ApJ*, **273**, 105

Bussmann, R. S., Pérez-Fournon, I., Amber, S., et al. 2013, *ApJ*, **779**, 25

Cañamas, R., Nesvadba, N. P. H., Guery, D., et al. 2015, *A&A*, **581**, A105

Cañamas, R., Nesvadba, N. P. H., Kneissl, R., et al. 2021, *A&A*, **645**, A45

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

Carilli, C. L., Daddi, E., Riechers, D., et al. 2010, *ApJ*, **714**, 1407

Carlstrom, J. E., Ade, P. A. R., Aird, K. A., et al. 2011, *PASP*, **123**, 568

Casey, C. M., Narayanan, D., & Cooray, A. 2014, *Phys. Rep.*, **541**, 45

Cernicharo, J., Liu, X. W., González-Alfonso, E., et al. 1997, *ApJ*, **483**, L65

Cernicharo, J., Goicoechea, J. R., & Caux, E. 2000, *ApJ*, **534**, L199

Chabrier, G. 2003, *ApJ*, **586**, L133

Clark, C. J. R., Schofield, S. P., Gomez, H. L., & Davies, J. I. 2016, *MNRAS*, **459**, 1646

Clark, C. J. R., De Vis, P., Baes, M., et al. 2019, *MNRAS*, **489**, 5256

Coleman, G. D., Wu, C. C., & Weedman, D. W. 1980, *ApJS*, **43**, 393

Comito, C., Schilke, P., Gerin, M., et al. 2003, *A&A*, **402**, 635

Cortzen, I., Magdis, G. E., Valentino, F., et al. 2020, *A&A*, **634**, L14

Cox, P., Krips, M., Neri, R., et al. 2011, *ApJ*, **740**, 63

da Cunha, E., Charlot, S., & Elbaz, D. 2008, *MNRAS*, **388**, 1595

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

da Cunha, E., Walter, F., Smail, I. R., et al. 2015, *ApJ*, **806**, 110

Daddi, E., Elbaz, D., Walter, F., et al. 2010, *ApJ*, **714**, L118

Daddi, E., Valentino, F., Rich, R. M., et al. 2020, *A&A*, submitted [arXiv:2006.11089]

Danielson, A. L. R., Swinbank, A. M., Smail, I., et al. 2011, *MNRAS*, **410**, 1687

Dekel, A., & Mandelker, N. 2014, *MNRAS*, **444**, 2071

Dekel, A., & Woo, J. 2003, *MNRAS*, **344**, 1131

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

Dhillon, V., Dixon, S., Gamble, T., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Proc. SPIE, **10702**, 107020L

Dong, C., Spilker, J. S., Gonzalez, A. H., et al. 2019, *ApJ*, **873**, 50

Downes, D., & Solomon, P. M. 1998, *ApJ*, **507**, 615

Downes, D., Neri, R., Greve, A., et al. 1999, *A&A*, **347**, 809

Draine, B. T. 2003, *ARA&A*, **41**, 241

Draine, B. T., & Li, A. 2001, *ApJ*, **551**, 807

Draine, B. T., & Li, A. 2007, *ApJ*, **657**, 810

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

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

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

Eales, S. A., Smith, M. W. L., Wilson, C. D., et al. 2010b, *A&A*, **518**, L62

Elbaz, D., Dickinson, M., Hwang, H. S., et al. 2011, *A&A*, **533**, A119

Erb, D. K. 2008, *ApJ*, **674**, 151

Erb, D. K., Shapley, A. E., Pettini, M., et al. 2006, *ApJ*, **644**, 813

Falgarone, E., Zwaan, M. A., Godard, B., et al. 2017, Nature, **548**, 430

Falstad, N., González-Alfonso, E., Aalto, S., et al. 2015, *A&A*, **580**, A52

Falstad, N., González-Alfonso, E., Aalto, S., & Fischer, J. 2017, *A&A*, **597**, A105

Fixsen, D. J., Bennett, C. L., & Mather, J. C. 1999, *ApJ*, **526**, 207

Fu, H., Cooray, A., Feruglio, C., et al. 2013, Nature, **498**, 338

Fudamoto, Y., Ivison, R. J., Oteo, I., et al. 2017, *MNRAS*, **472**, 2028

Fumagalli, M., Prochaska, J. X., Kasen, D., et al. 2011, *MNRAS*, **418**, 1796

Gao, Y., & Solomon, P. M. 2004, *ApJ*, **606**, 271

Genzel, R., Baker, A. J., Tacconi, L. J., et al. 2003, *ApJ*, **584**, 633

Genzel, R., Tacconi, L. J., Gracia-Carpio, J., et al. 2010, *MNRAS*, **407**, 2091

Godard, B., Falgarone, E., & Pineau des Forêts, G. 2014, *A&A*, **570**, A27

Goicoechea, J. R., Rodríguez-Fernández, N. J., & Cernicharo, J. 2004, *ApJ*, **600**, 214

González-Alfonso, E., Fischer, J., Graciá-Carpio, J., et al. 2012, *A&A*, **541**, A4

González-Alfonso, E., Fischer, J., Aalto, S., & Falstad, N. 2014, *A&A*, **567**, A91

González-Alfonso, E., Fischer, J., Spoon, H. W. W., et al. 2017, *ApJ*, **836**, 11

González-Alfonso, E., Fischer, J., Bruderer, S., et al. 2018, *ApJ*, **857**, 66

Gralla, M. B., Marriage, T. A., Addison, G., et al. 2020, *ApJ*, **893**, 104

Harris, A. I., Baker, A. J., Frayer, D. T., et al. 2012, *ApJ*, **752**, 152

Hensley, B. S., & Draine, B. T. 2021, *ApJ*, **906**, 73

Herrera-Camus, R., Sturm, E., Graciá-Carpio, J., et al. 2020, *A&A*, **633**, L4

Hodge, J. A., & da Cunha, E. 2020, *R. Soc. Open Sci.*, **7**, 200556

Ibar, E., Ivison, R. J., Best, P. N., et al. 2010, *MNRAS*, **401**, L53

Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, *A&A*, **457**, 841

Ilbert, O., Capak, P., Salvato, M., et al. 2009, *ApJ*, **690**, 1236

Indriolo, N., Neufeld, D. A., Gerin, M., et al. 2015, *ApJ*, **800**, 40

Indriolo, N., Bergin, E. A., Falgarone, E., et al. 2018, *ApJ*, **865**, 127

Ivison, R. J., Smail, I., Le Borgne, J. F., et al. 1998, *MNRAS*, **298**, 583

Ivison, R. J., Smail, I., Papadopoulos, P. P., et al. 2010, *MNRAS*, **404**, 198

Ivison, R. J., Swinbank, A. M., Smail, I., et al. 2013, *ApJ*, **772**, 137

Ivison, R. J., Page, M. J., Cirasuolo, M., et al. 2019, *MNRAS*, **489**, 427

James, A., Dunne, L., Eales, S., & Edmunds, M. G. 2002, *MNRAS*, **335**, 753

Jin, S., Daddi, E., Magdis, G. E., et al. 2019, *ApJ*, **887**, 144

Jones, A. P., Fanciullo, L., Köhler, M., et al. 2013, *A&A*, **558**, A62

Jones, A. P., Köhler, M., Ysard, N., Bocchio, M., & Verstraete, L. 2017, *A&A*, **602**, A46

Jones, G. C., Maiolino, R., Caselli, P., & Carniani, S. 2019, *A&A*, **632**, L7

Keeton, C. R. 2001, *ArXiv e-prints* [arXiv:astro-ph/0102340]

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

Kennicutt, R. C., Jr. 1998b, *ApJ*, **498**, 541

Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, *MNRAS*, **363**, 2

Lada, C. J. 2015, in *Galaxies in 3D across the Universe*, eds. B. L. Ziegler, F. Combes, H. Dannerbauer, & M. Verdugo, 309, 31

Lehner, N., Howk, J. C., Tripp, T. M., et al. 2013, *ApJ*, **770**, 138

Lequeux, J. 1983, *A&A*, **125**, 394

Leroy, A. K., Bolatto, A., Gordon, K., et al. 2011, *ApJ*, **737**, 12

Li, A., & Draine, B. T. 2001, *ApJ*, **554**, 778

Li, J., Wang, R., Riechers, D., et al. 2020, *ApJ*, **889**, 162

Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, *ApJ*, **772**, 119

Liu, D., Schinnerer, E., Groves, B., et al. 2019, *ApJ*, **887**, 235

Lu, N., Zhao, Y., Díaz-Santos, T., et al. 2017, *ApJS*, **230**, 1

Lutz, D., Poglitsch, A., Altieri, B., et al. 2011, *A&A*, **532**, A90

Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, *ApJ*, **760**, 6

Magnelli, B., Saintonge, A., Lutz, D., et al. 2012, *A&A*, **548**, A22

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

Marsden, D., Gralla, M., Marriage, T. A., et al. 2014, *MNRAS*, **439**, 1556

Meyer, D. M., & Roth, K. C. 1991, *ApJ*, **376**, L49

Müller, H. S. P. 2010, *A&A*, **514**, L6

Müller, H. S. P., Schlöder, F., Stutzki, J., & Winnewisser, G. 2005, *J. Mol. Struct.*, **742**, 215

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

Negrello, M., Amber, S., Amvrosiadis, A., et al. 2017, *MNRAS*, **465**, 3558

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

Oke, J. B., & Gunn, J. E. 1983, *ApJ*, **266**, 713

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

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

Oteo, I., Ivison, R. J., Dunne, L., et al. 2016, *ApJ*, **827**, 34

Oteo, I., Zhang, Z. Y., Yang, C., et al. 2017, *ApJ*, **850**, 170

Peeples, M. S., & Shankar, F. 2011, *MNRAS*, **417**, 2962

Peng, Y.-J., & Maiolino, R. 2014, *MNRAS*, **443**, 3643

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, *AJ*, **124**, 266

Perley, R. A., & Butler, B. J. 2017, *ApJS*, **230**, 7

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, *A&A*, **518**, L1

Planck Collaboration XXI. 2011, *A&A*, **536**, A21

Planck Collaboration VI. 2020, *A&A*, **641**, A6

Planck Collaboration Int. XXVII. 2015, *A&A*, **582**, A30

Planck Collaboration Int. XXVIII. 2015, *A&A*, **582**, A31

Planck Collaboration Int. XXIX. 2016, *A&A*, **586**, A132

Polletta, M., Tajer, M., Maraschi, L., et al. 2007, *ApJ*, **663**, 81

Priestley, F. D., & Whitworth, A. P. 2020, *MNRAS*, **494**, L49

Rangwala, N., Maloney, P. R., Glenn, J., et al. 2011, *ApJ*, **743**, 94

Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, *A&A*, **563**, A31

Reuter, C., Vieira, J. D., Spilker, J. S., et al. 2020, *ApJ*, **902**, 78

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

Riechers, D. A., Weißb., A., Walter, F., & Wagg, J. 2010, *ApJ*, **725**, 1032

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

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

Riechers, D. A., Leung, T. K. D., Ivison, R. J., et al. 2017, *ApJ*, **850**, 1

Rizzo, F., Vegetti, S., Fraternali, F., & Di Teodoro, E. 2018, *MNRAS*, **481**, 5606

Rubin, K. H. R., Prochaska, J. X., Koo, D. C., & Phillips, A. C. 2012, *ApJ*, **747**, L26

Sánchez Almeida, J., Elmegreen, B. G., Muñoz-Tuñón, C., & Elmegreen, D. M. 2014, *A&ARv*, **22**, 71

Schmidt, M. 1959, *ApJ*, **129**, 243

Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, *A&A*, **575**, A74

Sharon, C. E., Baker, A. J., Harris, A. I., et al. 2015, *ApJ*, **798**, 133

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

Solomon, P. M., & Vanden Bout, P. A. 2005, *ARA&A*, **43**, 677

Somerville, R. S., & Davé, R. 2015, *ARA&A*, **53**, 51
Spilker, J. S., Marrone, D. P., Aguirre, J. E., et al. 2014, ApJ, 785, 149
Spilker, J. S., Aravena, M., Béthermin, M., et al. 2018, Science, 361, 1016
Staniszewski, Z., Ade, P. A. R., Aird, K. A., et al. 2009, ApJ, 701, 32
Stanley, F., Knudsen, K. K., Aalto, S., et al. 2021, A&A, in press
https://doi.org/10.1051/0004-6361/202037888
Stone, M., Veilleux, S., Meléndez, M., et al. 2016, ApJ, 826, 111
Swinbank, A. M., Papadopoulos, P. P., Cox, P., et al. 2011, ApJ, 742, 11
Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179
Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157
Tagore, A. S., & Jackson, N. 2016, MNRAS, 457, 3066
Tagore, A. S., & Keeton, C. R. 2014, MNRAS, 445, 694
Thomson, A. P., Ivison, R. J., Simpson, J. M., et al. 2014, MNRAS, 442, 577
Tody, D. 1986, in Instrumentation in Astronomy VI, ed. D. L. Crawford, SPIE
Conf. Ser., 627, 733
Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017, ARA&A, 55, 389
Valentino, F., Daddi, E., Puglisi, A., et al. 2020, A&A, 641, A155
Valiante, E., Smith, M. W. L., Eales, S., et al. 2016, MNRAS, 462, 3146
van Dishoeck, E. F., Jansen, D. J., Schilke, P., & Phillips, T. G. 1993, ApJ, 416,
L83
Vieira, J. D., Crawford, T. M., Switzer, E. R., et al. 2010, ApJ, 719, 763
Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296
Weiß, A., Henkel, C., Downes, D., & Walter, F. 2003, A&A, 409, L41
Weiß, A., Walter, F., & Scoville, N. Z. 2005, A&A, 438, 533
Weiß, A., Downes, D., Neri, R., et al. 2007, A&A, 467, 955
Weiß, A., De Breuck, C., Marrone, D. P., et al. 2013, ApJ, 767, 88
Whitaker, K. E., Franx, M., Leja, J., et al. 2014, ApJ, 795, 104
Wiseman, P., Perley, D. A., Schady, P., et al. 2017, A&A, 607, A107
Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
Yabe, K., Ohta, K., Akiyama, M., et al. 2015, ApJ, 798, 45
Yang, C., Gao, Y., Omont, A., et al. 2013, ApJ, 771, L24
Yang, C., Omont, A., Beelen, A., et al. 2016, A&A, 595, A80
Yang, C., Omont, A., Beelen, A., et al. 2017, A&A, 608, A144
Yang, C., González-Alfonso, E., Omont, A., et al. 2020, A&A, 634, L3
Zhang, Z.-Y., Gao, Y., Henkel, C., et al. 2014, ApJ, 784, L31

<footnote n="1">¹ Institut de Radioastronomie Millimétrique (IRAM), 300 Rue de la Piscine, 38400 Saint-Martin-d'Hères, France
e-mail: berta@iram.fr</footnote>

<footnote n="2">² Department of Physics and Astronomy, Rutgers, The State University of New Jersey, 136 Frelinghuysen Road, Piscataway, NJ 08854-8019, USA</footnote>

<footnote n="3">³ Sorbonne Université, UPMC Université Paris 6 & CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98b Boulevard Arago, 75014 Paris, France</footnote>

<footnote n="4">⁴ Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, The School of Natural Sciences, The University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="5">⁵ School of Physics and Astronomy, Cardiff University, The Parade, Cardiff CF24 3AA, UK</footnote>

<footnote n="6">⁶ Instituto Astrofísica de Canarias (IAC), 38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="7">⁷ Universidad de La Laguna, Dpto. Astrofísica, 38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="8">⁸ University of Cape Town, Department of Astronomy, Private Bag X3, Rondebosch 7701, Cape Town, South Africa</footnote>

<footnote n="9">⁹ Department of Physics and Astronomy, University of the Western Cape, Private Bag X17, Bellville 7535, Cape Town, South Africa</footnote>

<footnote n="10">¹⁰ Istituto Nazionale di Astrofisica, Istituto di Radioastronomia, Via Gobetti 101, 40129 Bologna, Italy</footnote>

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

<footnote n="12">¹² Department of Astronomy, Cornell University, Space Sciences Building, Ithaca, NY 14853, USA</footnote>

<footnote n="13">¹³ Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="14">¹⁴ Leiden University, Leiden Observatory, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="15">¹⁵ Division of Particle and Astrophysical Science, Graduate School of Science, Nagoya University, Aichi 464-8602, Japan</footnote>

<footnote n="16">¹⁶ National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="17">¹⁷ European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching, Germany</footnote>

<footnote n="18">¹⁸ Aix-Marseille Université, CNRS & CNES, Laboratoire d'Astrophysique de Marseille, 38, Rue Frédéric Joliot-Curie, 13388 Marseille, France</footnote>

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

<footnote n="20">²⁰ School of Physics and Astronomy, University of Nottingham, University Park, Nottingham NG7 2RD, UK</footnote>

<footnote n="21">²¹ Department of Astronomy, University of Maryland, College Park, MD 20742, USA</footnote>

<footnote n="22">²² Instituto Nacional de Astrofísica, Óptica y Electrónica, Astrophysics Department, Apdo 51 y 216, Tonantzintla, Puebla 72000, Mexico</footnote>

<footnote n="23">²³ Instituto de Astrofísica e Ciâncias do Espaco, Tapada da Ajuda, Edifício Leste, 1349-018 Lisboa, Portugal</footnote>

<footnote n="24">²⁴ Department of Physical Sciences, The Open University, Milton Keynes MK7 6AA, UK</footnote>

<footnote n="25">²⁵ National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA</footnote>

<footnote n="26">²⁶ Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany</footnote>
# Appendix A: Details of lens model

The statistical uncertainties on the lens model parameters, for both the foreground lens object and the HerBS-89a lensed source were computed via two MCMC runs, focusing on the lens and the source, respectively (see Sect. 5). Figures A.1 and A.2 show the posterior distribution of the lens and source parameters, as obtained with MCMC.

In order to perform the kinematic analysis described in Sect. 7, we apply the lens model to high resolution emission line

data preserving the velocity information by fitting Sérsic profiles to channel images in the data cubes. To obtain a sufficient S/Ns, we bin the 40 km s⁻¹ channels of ¹²CO(9–8) and H₂O by a factor of 3 and 6, respectively. Despite this rebinning, the S/Ns for some of the channels is still very low.

Figures A.3 and A.4 present the results of this channel by channel modeling, for the ¹²CO(9–8) and H₂O(2₀₂–1₁₁) emission lines, respectively.

<figure><img src="image_21.png" /><figcaption>**Fig. A.1.** Posterior distribution of the parameters of the foreground lens object, as obtained by the lens modeling and MCMC runs (see Sect. 5 for details). The constraints are given on the following parameters: the position (ΔRA and ΔDec in arcsec), the Einstein radius (b), the ellipticity (e), and the position angle (θe) with the orientation measured east of north. The contours show the 0.5, 1, 1.5, and 2σ of the integrated likelihood. The mean and the ±1σ interval around the mean are displayed as histograms along the diagonal of the plot.</figcaption></figure>
<figure><img src="image_22.png" /></figure>

**Fig. A.2.** Posterior distribution of the parameters of the HerBS-89a lensed source, as obtained by the lens modeling and MCMC runs (see Sect. 5 for details). The constraints are given on the following parameters: the position (ΔRA and ΔDec in arcsec), the Einstein radius (b), the ellipticity (e), the scale radius (R_s), the position angle (θ_e) with the orientation measured east of north and the Spheric index (n). See Fig. A.1 for further details.
<figure><img src="image_23.png" /><figcaption>**Fig. A.3.** Lens modeling of the $^{12}$CO(9–8) spectral channels. The velocity ranges covered by each channel are marked in the right-hand panels. From left to right: observed image (with the synthesized beam shown in the lower left corner) and critical curve (red line); reconstructed model image convolved with the PSF; residuals, reconstructed model image at full angular resolution; and, separated and with a different angular scale, reconstructed source-plane images and caustic curve. The cross shows the central coordinate RA 13:16:11.52 and Dec +28:12:17.7 (J2000).</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Fig. A.4.** Lens modeling of the H₂O(2₀₂−1₁₁) spectral channels. The velocity ranges covered by each channel are marked in the right-hand panels. See Fig. A.3 for further details.</figcaption></figure>
## Appendix B: Review of $\kappa_v$ values

Although the use of the MBB model and Eq. (5) to derive dust masses is common and well-established in the literature, the actual value of $\kappa_v$ is far from being well defined. Moreover, it is not granted that a universal value shared by all galaxies does exist; evidence rather points to the opposite. We here review the current status on $\kappa_v$ taking into account recent studies.

The common preference found in the literature is to adopt the $\kappa_v$ values tabulated by Li & Draine (2001) or its updated version by Draine (2003). At 850 µm, the reference value is $\kappa_{850} = 0.038 \text{ m}^2 \text{ kg}^{-1}$ with a dependence on frequency $\sim \nu^{2.07}$ at $\lambda = 100$ and 600 µm (Berta et al. 2016). To reconcile with a small conflict between the dust to hydrogen mass ratio predicted by their model and the elemental abundances in the solar neighborhood ISM, Draine et al. (2014) re-normalized their dust properties such that $\kappa_{850} = 0.047$ and $\kappa_{500} = 0.116 \text{ m}^2 \text{ kg}^{-1}$. Hensley & Draine (2021) discuss more recent finding, that could lead to further possible corrections toward larger $\kappa_v$. A larger value $\kappa_{850} = 0.072 \text{ m}^2 \text{ kg}^{-1}$ come also from the THEMIS dust model (Jones et al. 2013, 2017), which is in agreement with the measured all sky extinction from SDSS quasars (Planck Collaboration Int. XXIX 2016).

The tabulated $\kappa_v$ values mentioned above are theoretical (see, e.g., Weingartner & Draine 2001) and an alternative path is to adopt an empirical value of $\kappa_v$, measured at a reference frequency and then scale it to any other frequency using the power law with index $\beta$ as usual. Studying a sample of 22 local galaxies with known metallicities and available sub-mm data, James et al. (2002) derived $\kappa_{850} = 0.07 \pm 0.02 \text{ m}^2 \text{ kg}^{-1}$.

Using Herschel observations of local galaxies M99 and M100, Eales et al. (2010b) derive a value of $\kappa_v \sim 0.06 \text{ m}^2 \text{ kg}^{-1}$ at 350 µm, corresponding to $\kappa_{850} = 0.01 \text{ m}^2 \text{ kg}^{-1}$ adopting a dust emissivity index $\beta = 2$. From early solar-neighborhood Planck results, a value $\kappa_{850} = 0.066 \text{ m}^2 \text{ kg}^{-1}$ is derived<fnref n="9" />, roughly twice the commonly adopted value by Draine (2003). Recently, Clark et al. (2016) assembled a sample of 22 galaxies from the Herschel Reference Survey (Boselli et al. 2010) and derived a value of $\kappa_{500} = 0.051_{-0.026}^{+0.070} \text{ m}^2 \text{ kg}^{-1}$ at 500 µm, corresponding to $\kappa_{850} = 0.018_{-0.01}^{+0.023}$. Significant variations of $\kappa_v$ at sub-galactic scales are shown by Clark et al. (2019).

Such very low values of $\kappa_v$ obtained for external galaxies are nevertheless subject to several caveats related to the derivation method, and thus need to be taken with caution. To mention a few: the sight-lines intercept a range of media, and possibly multiple temperature dust components; the assumptions made about dust-to-metals ratios need solid metallicity measurements and might depend on the adopted calibration; a variable $X_{CO}$ brings further uncertainties. For further discussions, see for example Priestley & Whitworth (2020).

Finally, Clark et al. (2016, 2019) present a thorough compilation of the values of $\kappa_v$ found in the literature to date, and remarkably collect values spanning a range of more that three orders of magnitude (at their reference 500 µm wavelength; see the respective Fig. 1 in the above mentioned papers).

With these results in mind, it is clear that a proper comparison of dust masses to literature values must take into account the different assumptions underlying the SED modeling. In order to follow a popular choice in this field, we here adopted the Draine (2003) $\kappa_v$, applying the correction prescribed by Draine et al. (2014) to dust masses (as derived with both the MBB and DL07 modeling).

<footnote n="9">⁹ Planck Collaboration XXI (2011) measure $\tau_{250}/N_H = 0.92 \times 10^{-25} \text{ cm}^2$ for the solar neighborhood. Using the relation $\tau_v = \kappa_0 (\nu/\nu_0)^\beta R_{DG}\mu_H N_H$ (Bianchi et al. 2019; Planck Collaboration Int. XXVIII 2015), this leads to $\kappa_{850} = 0.066 \text{ m}^2 \text{ kg}^{-1}$. Here $\beta = 1.8$ and a dust to hydrogen mass ratio $R_{DG} = 0.0091$ (Draine et al. 2014) were adopted. The quantity $\mu_H$ is the mass of the hydrogen atom. It is worth to note the strong dependence of $\kappa_v$ on $\beta$ and $R_{DG}$: different choices lead to significantly different values of $\kappa_{850}$. Typical alternatives to the values mentioned above are, for example, $\beta = 2$ (often found in the literature) and $R_{DG} = 0.0074$ (Jones et al. 2017). See Bianchi et al. (2019) for an in depth description.</footnote>