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

(Abridged) Exploiting the sensitivity and broad band width of NOEMA, we have studied the molecular gas and dust in the galaxy HerBS-89a, at z=2.95. High angular resolution images reveal a partial 1.0"diameter Einstein ring in the dust continuum emission and the molecular emission lines of 12CO(9-8) and H2O(2_02-1_11). We report the detection of the three fundamental transitions of the molecular ion OH+, seen in absorption; the molecular ion CH+(1-0) seen in absorption (and tentatively in emission); two transitions of amidogen (NH2), seen in emission; and HCN(11-10) and/or NH(1_2-0_1) seen in absorption. The NOEMA data are complemented with VLA data tracing the 12CO(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 HST imaging that reveals the foreground lensing galaxy in the near-infrared. Together with data from the GTC, we derive a photometric redshift of z(phot)~0.9 for the foreground lensing galaxy. Modelling the lensing of HerBS-89a, we reconstruct the dust continuum and molecular emission lines (magnified by a factor ~4-5) in the source plane. The 12CO(9-8) and H2O 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 accounts for the observed broad line widths. HerBS-89a is a powerful star forming galaxy with a dust-to-gas ratio delta(GDR)~80, a SFR = 614 +/- 59 Msun/yr and a depletion timescale tau(depl) = (3.4 +/- 1.0) 1e8 years. The OH+ and CH+ absorption lines, all have their main velocity component red-shifted by \Delta(V)~100 km/s relative to the global CO reservoir. We argue that these absorption lines trace a rare example of gas inflow towards the center of the galaxy.


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 starforming galaxies from hundreds to several hundred thousand. Together with other large-area surveys, like the allsky Planck-HFI (Planck Collaboration et al. 2015a), 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), we have today 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. 2010Negrello et al. , 2017Cox et al. 2011;Bussmann et al. 2013;Spilker et al. 2014;Cañameras et al. 2015;Reuter et al. 2020) and rare cases of galaxies with intrinsic infrared luminosities, L FIR 10 13 L ⊙ , and star formation rates (SFRs) in excess of 1000 M ⊙ yr −1 , known as Hyper-Luminous Infrared Galaxies (HyLIRGs, see, e.g., Ivison et al. 1998Ivison et al. , 2013Ivison et al. , 2019Fu et al. 2013;Oteo et al. 2016;Riechers et al. 2013Riechers et al. , 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 sub/millimeter spectroscopy has proved to be the most reliable method for measuring accurate redshifts, an approach that 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;Fudamoto et al. 2017;Neri et al. 2020;Reuter et al. 2020, and references therein). A&A proofs: manuscript no. berta_HerBS89a_v2p1 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µ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 Programme, 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 Programme sample display CO emission line widths in excess of 800 km s −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 159GHz = 4.56 ± 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 −1 . Both the 2-mm continuum and 12 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 13 L ⊙ , which suggested that HerBS-89a could be a HyLIRG.
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 12 CO(9 − 8), H 2 O (2 02 − 1 11 ), NH 2 and HCN(11−10) (overlapping with NH(1 2 −0 1 )), to the low H 2 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 12 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; Section 3 presents the main results, including the morphology of the source, the properties of the molecular emission and absorption lines and the underlying continuum; Section 4 describes the characteristics of the foreground lensing galaxy; Section 5 presents the lensing model and the morphology of HerBS-89a in the source plane; and Section 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 2 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, Section 10 summarizes the main conclusions of this paper and outlines future prospects.
Throughout this paper we adopt a spatially flat ΛCDM cosmology with H 0 = 67.4 km s −1 Mpc −1 and Ω M = 0.315 (Planck Collaboration et al. 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 × 10 4 Mpc.

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, labelled 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 hours on-source time. This project, which was specifically tailored to measure the 12 CO(9 − 8) and H 2 O(2 02 − 1 11 ) emission lines, also enabled the detection of two strong absorption lines of the molecular ion OH + corresponding to the redshifted ground state transitions of OH + (1 1 −0 1 ) and OH + (1 2 −0 1 ) (ν rest = 1033.118 GHz and 971.803 GHz, respectively) 1 .
The second project, labelled 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 hours. This project was a follow-up of the DDT to measure the third ground state transition of OH + (1 0 − 0 1 ) (ν rest = 909.158 GHz) 1 and the ground state transition of the molecular ion CH + (1 − 0) (ν rest = 835.137 GHz) 2 .
Observing conditions were excellent for both projects with an atmospheric phase stability of 20 o 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 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 LkHα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σ error on the absolute position of HerBS-89a is estimated to be < 0 ′′ . 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, i.e., 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 hours (with 62 minutes 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 −1 ) 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 used 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 (P.A. −62.3 o ) 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%.

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 characterising 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 seconds. Data were reduced using the Python AstroDrizzle package, with parameters optimised to improve the final image quality. The pixel scale of the Infrared-Camera is 0 ′′ . 128, but the image was resampled to  Neri et al. (2020). 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σ.
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.

GTC
HerBS-89a was observed with the GTC 10.4 m telescope in two observing runs (program ID: PI: H. Dannerbauer). First, optical imaging was obtained using the instrument OSIRIS on April 6th, 2019 . The observations were conducted in service mode under clear skies (although with non-photometric conditions), integrating for 10 minutes 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 7th, 2020 using the visitor instrument HiPERCAM for a total observing time of 1.25 hours under average seeing conditions (0 ′′ . 9). HiPERCAM is a quintuple-beam CCD imager enabling the five Sloan filter ugriz to be observed simultaneously over a field of view of 2 ′ .8×1 ′ .4 with a pixel size of 0 ′′ . 081 (Dhillon et al. 2018, and references therein). The HiPERCAM team designed a A&A proofs: manuscript no. berta_HerBS89a_v2p1  dedicated data reduction tool 3 that applies standard procedures for the reduction and calibration of the raw images.

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 highangular resolution observations for both the continuum and the 3 http://deneb.astro.warwick.ac.uk/phsaap/hipercam/docs/html/ 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.3). 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 σ = 29 and 33 µJy/beam in the two bands, respectively. Merging the two side-bands of the high-resolution continuum images improves the S/N ratio and resulted, by applying a uniform weight, in a final image with beam size of 0 ′′ . 43 × 0 ′′ . 22, to be compared to the 0 ′′ . 55 × 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.

Continuum emission
For the lower frequency data obtained using the mediumcompact 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 × 0 ′′ . 66.
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 catalogue at that position. The faintness of 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 −1 for the 245-268 GHz data and 60 km s −1 for the 209-233 GHz data, reaching r.m.s. 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), i.e., o-H 3 O + (0 − 0 − 1 + 0 ), H 2 O + (3 12 − 3 03 ), SH + (2 − 1), and o-NH 2 (2 0,2 − 1 1,1 )(3/2 − 1/2).
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. Note that the latter should be added in quadrature to the 10% absolute flux calibration uncertainty (see Sect. 2.1). 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σ 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. The upper right panel presents the 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. The right panel in the second row displays the 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.

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, Fig. 5. The spectra of the molecular emission and absorption lines detected in HerBS-89a reported in this study, together with the 12 CO(5 − 4) and 12 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, i.e., 12 CO(9 − 8)/OH + (1 1 − 0 1 ), OH + (1 2 − 0 1 )/HCN(11 − 10) and OH + (1 0 − 0 1 )/o − NH 2 (2 02 − 1 11 )(5/2 − 3/2), the second molecular line is labelled 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).
which are detected for the first time in a high-z galaxy. Blueshifted from the strong OH + (1 2 − 0 1 ) absorption line, another weaker absorption line is detected that is due to HCN(11 − 10) (ν rest = 974.487 GHz) and/or NH(1 2 − 0 1 ) (ν 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) (ν rest = 907.433 GHz). A second NH 2 line is detected with lower signal-to-noise ratio at higher frequency, the NH 2 (2 20 − 1 11 )(5/2 − 3/2) transition at ν 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 signal-to-noise 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) (ν rest = 893.152 GHz), o-NH 2 (2 02 − 1 11 )(3/2 − 1/2) (ν rest = 902.210 GHz), H 2 O + (3 12 − 3 03 ) (ν rest = 982.955 GHz), and o-H 3 O + (0 − 0 − 1 + 0 ) (ν 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 signal-to-noise ratio 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 towards 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, A&A proofs: manuscript no. berta_HerBS89a_v2p1 including the 12 CO(3−2) and 12 CO(5−4) emission lines reported in Neri et al. (2020).
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 ratio 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 ∆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 frequencies. For the OH + (1 0 − 0 1 ) absorption, this shift is larger than for the two other OH + transitions, and in the case of NH 2 , the profile seems to be inconsistent with those seen for other molecular emission lines tracing dense gas (CO and H 2 O). Since NH 2 (amidogen) is a photodissociation product of NH 3 (ammonia), it should probe the dense gas traced in the high-J CO or water lines, and therefore display a similar double Gaussian profile. However, we only see the red-shifted component of the 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) emission 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 width at half maximum (FWHM), the integrated intensity, and the intrinsic frequency. All the uncertainties are propagated into the derived quantities.

Optical depths and column densities of lines seen in absorption
The optical depth of an absorption line (τ ν ) can be evaluated from the observed spectrum normalized by the continuum as follows where S ν is the flux density of the spectrum and S cont ν that of the continuum only. The optical depth is evaluated at the central frequency ν 0 of the line. Table 3 reports the values of τ(ν 0 ) obtained from the integrated spectra shown in Fig. 3, i.e., averaged over the angles covered by the lines on the sky. In the case of the two strongest absorption lines (the OH + (1 1 − 0 1 ) and (1 2 − 0 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 τ ν 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;: . The velocity maps were obtained above 3σ 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. The upper right panel presents the 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. The right panel in the second row displays the 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.
The quantity τ ν dV and the speed of light c can be expressed in units of [cm s −1 ], and the transition frequency ν in [Hz] (see Table 2 for the corresponding frequencies of each line), so that N l is directly computed in units of [cm −2 ]. 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 ∆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, CH + (1 − 0) and the absorption line observed at 246.72 GHz that is due to HCN(11−10) and/or NH(1 2 −0 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.

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 12 CO(9−8) and H 2 O(2 02 −1 11 ) 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 12 CO(9−8) and H 2 O(2 02 −1 11 ) 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 disc-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 (2 02 − 1 11 )(5/2 − 3/2) emission line is clearly seen 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.3. 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/N ratio, 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 ∆V ∼ 1400 km s −1 , comparable to those of the three high-J CO emission lines, and the observed total line intensity is I ν = 0.64 ± 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 S. Berta et al.: Close-up View of a z = 2.95 Star-Forming Galaxy +500 km s −1 relative to the frequency of 12 CO(1 − 0) at z = 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σ 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).

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. Figure 9 displays the HST wide-J filter F110W image (at 1.15 µm) towards 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 North-West 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 data, the foreground galaxy has an estimated flux density of ∼ 4.36 ± 0.88 µJy at 1.15 µm.

Photometric redshift (GTC)
The foreground lensing galaxy is also detected in our GTC OSIRIS r-band imaging with a magnitude r = 24.5 ± 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σ 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 µ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 µ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 an photometric redshift.
A&A proofs: manuscript no. berta_HerBS89a_v2p1 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 WIS E 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 χ 2 as a function of redshift; the dotted horizontal lines correspond to the 1, 2, 3 σ levels, defined as χ 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σ.
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σ < 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 χ 2 distribution. The photometric redshift is best constrained by the r−i color, tracing the D4000 break 4 . The lack of rest-frame nearinfrared 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 ∼ 0.7, namely [Oii] (λ3727) and Hα (λ6563), but none produced a better fit to the data.

Spectroscopic redshift (NOEMA)
An independent indication of the redshift of the lensing foreground can possibly be derived from the 2-mm spectral scan that 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σ level. was obtained to derive the spectroscopic resdhift of HerBS-89a in Neri et al. (2020). In the frequency range 139.3 to 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σ 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.

Lens Modelling
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 observations) and proceed to apply it to the molecular line data. We perform the lens modelling and source reconstruction using Notes. The position angle's orientation is measured east of north.
the lensmodel package (Keeton 2001) with the pixsrc extension (Tagore & Keeton 2014;Tagore & Jackson 2016). For the sake of computational efficiency, all modelling 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 non-parametric approach of 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 non-parametric 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), R.A., Dec., ellipticity (e), 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 , elliptic-5 Ellipticity is defined as e = 1 − q, where q is the axis ratio. ity e, position angle θ e , scale radius R s , and Sérsic index n: 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 lens = 9.8 +3.2 −4.9 × 10 10 M ⊙ (where the upper and lower bounds are at 3σ). 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 HPBW min HPBW maj 2 in terms of the synthesized beam dimensions HPBW min and HPBW max , the S/N of the observed detection σ I /I line , and the magnification µ (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 2 O(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 note that 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 µ, for the high resolution data and for 12 CO(1 − 0). Because of the lower resolutions and S/N ratios, we do not model the spectral lines and continuum that were measured at lower frequencies (from 209 to 232 GHz).
The magnification factor µ 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 ± 0.03(±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 CO(1 − 0) line, we infer a lensing magnification µ = 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 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.
A&A proofs: manuscript no. berta_HerBS89a_v2p1 Notes. (a) 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, θ e for the position angle (measured east of north), R s for the scale radius, n for the Sérsic index, and µ for the amplification factor. Surface brightness units of Jy/beam are inherited from the surface brightness units of the data. Strong covariance between R S and n (see, e.g., Fig. A.1) means that the 12 CO(9 − 8) and 12 CO(1 − 0) data are formally compatible with (larger) R S and (smaller) n like those recovered for the other reconstructions.

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 Sections 6.1 and 6.2; the CO line excitation is investigated in Sect. 6.3, and the properties of the water emission line are examined in Sect. 6.4.

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(Bakx et al. , 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 modelling, 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 , , 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, e.g., the possible contribution of the polycyclic 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 dust is described as: where B ν (T dust ) is the Planck function, T dust is the dust temperature, and κ ν = κ 0 (ν/ν 0 ) β is the mass absorption coefficient of dust at rest frequency ν.
The general form of the MBB differs from the optically thin approximation by the factor (1 − e −τ ν ), with τ ν = τ 0 (ν/ν 0 ) β , instead of the ν β term. Here ν 0 is the frequency at which τ = τ 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 , β and the model normalization. In the general form, ν 0 is also a free parameter.
The MBB fits are limited to the data with rest-frame wavelength λ rest > 50 µm, thus avoiding biases towards 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 κ ν and β. Appendix B summarizes the current κ ν panorama. Following the choice often found in the literature, we adopt the Draine (2003) κ ν and apply the correction prescribed by Draine et al. (2014): κ ν = 0.047 m 2 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 β. 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σ 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 ν ∝ ν −0.8 is plotted, 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).

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 of the molecular gas mass can be made, using the following relation (Solomon & Vanden Bout 2005): where S Line ∆V = I Line is the integrated line flux in units of Jy km s −1 ; ν rest is the line rest frequency in GHz; and D L is the luminosity distance in Mpc. The line luminosity, expressed in units of L ⊙ , can be derived as: We define the molecular and total gas masses as where M He is the mass of helium, M H 2 is the mass of molecular hydrogen and M HI the mass of atomic hydrogen. Hereafter, we will assume that M mol ≫ M HI and thus that M gas ≈ M mol . The molecular hydrogen gas mass is computed from the 12 CO(1 − 0) luminosity applying a conversion factor α CO in units of M ⊙ (K km s −1 pc 2 ) −1 : In the Milky Way and nearby star-forming galaxies with near-solar metallicity, the empirical conversion factor α CO is α MW = 4.4 ± 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  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σ 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. Lower panel: HerBS-89a NOEMA spectra (Neri et al. 2020, and this work) compared to the general form MBB model with CMB correction. starbursts, there is good evidence from dynamical arguments that α 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 α CO =0.8 (Carilli & Walter 2013). Applying a 1.36× correction factor to account for the helium contribution, the latter value becomes α CO = 1.09. Dunne et al. (2020) 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 α CO . In another paper, Dunne et al. (in prep.) combined this sample will 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 κ 850 = 0.065 m 2 kg −1 is adopted (in agreement with the results by Planck Collaboration et al. 2011, see Appendix B for details.) and it is found that a value α CO ≃ 3 is consistent with the behaviour 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 α CO = 3.0 M ⊙ (K km s −1 pc 2 ) −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 ′ CO(1-0) = (5.15 ± 1.05) × 10 10 K km s −1 pc 2 , 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 (δ GDR ) in HerBS-89a: where the assumption is made that the two estimates are representative for the whole galaxy, and that M mol ≫ M HI . The data lead to a value of δ GDR between 80 and 174, depending on whether the DL07 or the MBB (general form) dust mass estimate is used, with an uncertainty of ±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 δ GDR ∼ 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 α CO The total molecular gas mass (M mol ) and SFR relate to each other via the so-called depletion timescale (in years): where the rate of star formation is computed as SFR = 1.09 × 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 ± 59 M ⊙ yr −1 and a depletion timescale τ depl = (3.4 ± 1.0) × 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 α CO = 3.0 × 1.36, a value that is close to the Milky Way value α MW = 4.4 M ⊙ (K km s −1 pc 2 ) −1 , commonly used for MS galaxies, rather than α SB = 1.09, often associated to starbursts or outliers in the main sequence. This choice leads to δ GDR and τ 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 ∼ 600 M ⊙ 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 α SB , the gas-to-dust ratio and the depletion time scale of HerBS-89a would both decrease by a factor of 3.74, yielding δ GDR ≃ 22 − 47 and τ depl = 9.1 ± 2.7 × 10 7 yr. Realistically, the actual value of α 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.

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 velocityintegrated 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 M82 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 µ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" starforming 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 ∼5-6 and a falloff towards 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 gradient (N CO /dV), and the solid angle (Ω app ) of the source. The 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 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 Cortzen et al. in prep.); M82 (Weiß et al. 2005); and the Milky Way disk (Fixsen et al. 1999).
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 /cm −3 ) = 2.79 +0.72 −0.54 , the kinetic temperature T kin = 174 +157 −92 K and the CO column density log 10 (N CO/dV /(cm −2 km s −1 )) = 17.59 +0.55 −0.59 , where the errors are 3σ.
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 +1.06 −0.57 cm −3 , T kin = 20 +20 −7 K, and log 10 (N CO/dV ) = 15.63 +1.14 −0.79 cm −2 km s −1 , where the errors are 3σ. The warm component is in broad agreement with the previous single component fitting results.

The para-H 2 O(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 highredshift galaxies (e.g., Omont et al. 2013;Riechers et al. 2013;Yang et al. 2016, 2020, andreferences therein). In particular, since the E upper > 200 K levels of H 2 O 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, submillimeter water lines are the second strongest molecular emitter after the CO lines (e.g., Yang et al. 2013). The detection of the para-H 2 O( 202111) emission line (ν rest = 987.9268 GHz) in HerBS-89a is an illustration thereof: the line width is broad (∆V ∼ 1100 km s −1 ) and comparable to those of the mid/high-J CO lines, consistent with the finding that the J = 2 H 2 O 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 2 O = 1.59 ± 0.37 Jy km s −1 , and comparable to the line flux of the 12 CO(9 − 8) emission line ( Table 2). The CO and H 2 O 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 2 O(2 02 −1 11 ) = (6.4 ± 1.5) × 10 7 L ⊙ , yielding a ratio L H 2 O(2 02 −1 11 ) /L IR = (1.12 ± 0.36) × 10 −5 , in terms of the 8 to 1000 µm L IR (Table 7). HerBS-89a lies on the correlation between L IR and L H 2 O , i.e. L H 2 O(2 02 −1 11 ) ∝ L 1.06±0.19 IR 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 sub-millimeter H 2 O emission that originates in very dense, heavily obscured starformation 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 ∼ 100 µm) with increasing infrared luminosity. Fi-

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 Section 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 play Fig. 16. Correlation between L IR and L H 2 O(2 02 −1 11 ) 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). 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 in the Appendix) some of the results of our lens modelling efforts. Our main focus is the 12 CO(9-8) line, which has higher S/N than the H 2 O(2 02 − 1 11 ) line and is therefore more promising as a kinematic probe. In order to improve the S/N, we rebin the 40 km s −1 data to a resolution of 120 km s −1 , yielding a total of seven independent velocity channels that are suitable for modelling, 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 Figure 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 −1 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 12 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 recon-A&A proofs: manuscript no. berta_HerBS89a_v2p1 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. structions 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érsicprofile 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 preparation) 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 Sec. 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 recon-structions 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 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 where i is the inclination angle of the rotating disk with respect to the line of sight, and V 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 ≃ 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 dyn sin 2 i ≈ 5.1 × 10 10 M ⊙ 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 two regions of strongest integrated emission (compare also to Figs. 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/20 th . 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. 17 and A.3). These two regions' velocity scales have been artificially offset by ±289 km s −1 , 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 −1 shift.
The individual dynamical masses of these components, assuming rotational support as above, would be M dyn sin 2 i ≈ 2.6 × 10 10 M ⊙ within R = 0.85 kpc for the eastern component and M dyn sin 2 i ≈ 2.2 × 10 10 M ⊙ 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.
A&A proofs: manuscript no. berta_HerBS89a_v2p1 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 HI can be neglected (see, e.g., Tacconi et al. 2018), the dynamical mass could yield an estimate of the stellar mass M * = f baryon M dyn − M mol − M HI and of the molecular gas fraction f mol = M mol /( f baryon M dyn − M 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 2 O(2 02 −1 11 ) emission line, adopting the same approach as above, using a coarser velocity re-binning to boost the S/N ratio (see Appendix A). As already suggested by the slightly different observed and reconstructed source-plane zeroth moment maps (see Fig. 12), the individual 12 CO(9 − 8) and H 2 O(2 02 − 1 11 ) 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 2 O(2 02 − 1 11 ) channel maps shows a velocity gradient similar to the one shown in the top left panel of Fig. A.3 for the 12 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 12 CO(9 − 8) emission line, we have limited the kinematic modelling 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 12 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.

Molecular lines other than CO and water
In the following sub-sections, we present and discuss the various molecular absorption and emission lines other than CO and H 2 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 2 (Sect. 8.1), followed by the two molecular ions, OH + and CH + , which probe the low density molecular gas with low H 2 fractions (Sect. 8.2).
As the data are of low signal-to-noise, a detailed study of NH 2 , first reported here in a high-z galaxy, remains out of reach. However, from the velocity-integrated emission map of the highangular resolution data of the NH 2 (2 20 −2 11 )(5/2−3/2) emission line (Fig.4), it is clear that the NH 2 emission closely resembles the distributions of the 12 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 2 − 0 1 ). With the current data, there is no way to distinguish between these molecular species, so we will briefly describe each molecule separately. The absorption feature is relatively weak and detected at a low signal-to-noise ratio. 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 H 2 3×10 4 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;Canameras et al. 2020, 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 6 . These transitions of HCN have critical densities around 10 9 − 10 10 cm −3 and probe very dense (n H 2 ≈ 3 × 10 6 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 τ(ν 0 ) ∼ 0.24 ± 0.12 yielding, in the case of HCN, a column density of (0.38 ± 0.1) × 10 14 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 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 2 ) to the formation and destruction of NH 3 . 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 2 −1 0 ) with OH + (1 0 −0 1 ), 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 14 cm −2 , which is also comparable to the column density of OH + ( Table 3). Together with NH 2 , 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 3 , 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) andMrk 231 (González-Alfonso et al. 2018), and in highz galaxies with NOEMA in HLFS 3 (Riechers et al. 2013) and ALMA (Indriolo et al. 2018, who reported OH + (1 1 − 0 1 ) 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 2 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 2 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 1 − 0 1 ) 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. (2020), 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 1 lines (J = 0, 1, 2) are seen in absorption (Fig. 5 middle panel). The (1 1 −0 1 ) and (1 2 −0 1 ) transitions are the strongest, while the (1 0 −0 1 ) 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 ∆V ≈ 500 km s −1 and the inferred OH + col-umn density is N OH + ∼ 10 15 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 ≈ +100 km s −1 relative to the systemic velocity, including the (1 0 − 0 1 ) transition after the contamination of the NH 2 emission line is accounted for (see Sec. 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 ≈ 2.5, including the Cosmic Eyelash and SDP.17 7 . 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 (∼ 10 7 cm −3 for the J = 1 − 0 transition), the CH + line is seen in emission in dense (> 10 5 cm −3 ) gas, probing shocked regions powered by galactic winds, and in absorption, in large (> 10 kpc) reservoirs of turbulent, cold (∼ 100 K), lowdensity (∼ 100 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 ≈ 400 km s −1 , comparable to that of the OH + lines, and is also red-shifted to ≈ +100 km s −1 (Fig. 5). The inferred CH + column density is estimated to be N CH + ∼ 10 14 cm −2 (Table 8), similar to what has been derived in the case of the Cosmic Eyelash (Falgarone et al. 2017).
The low signal-to-noise, tentative detection of CH + emission could indicate that in HerBS-89a there is also dense gas (> 10 5 cm −3 ) with a significant velocity dispersion (the line has a full width at zero intensity of 1000 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.

Inflowing gas
In HerBS-89a, the three OH + ground-state 1 J − 0 1 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 ≈ +100 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 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 8 shows A&A proofs: manuscript no. berta_HerBS89a_v2p1 a red-shifted line, namely the Cosmic Eyelash (Falgarone et al. 2017); for this source, the CH + (1 − 0) absorption line is shifted by ∆V ∼ +150 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 1 − 0 1 ) and H 2 O + (1 11 − 0 00 ) absorption lines (Indriolo et al. 2018). Recently, Butler et al. (in prep.) report OH + (1 1 − 0 1 ) absorption lines towards ten gravitationally lensed Herschel-selected highz 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 DS-FGs in the OH 119 µm doublet, in SPT-2319−55 at z = 5.29 (Spilker et al. 2018)  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 Eyelash (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 2 O or H 2 O + ) 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 ∼ 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 ∼ 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-α 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 starburstspect 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.). Theoretically, gas flowing into less massive (M ≤ 10 12 M ⊙ ) 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 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 ∼ 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 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 inflow , the mass inflow rate can be computed as follows: Here X denotes either the OH + or CH + molecular ion; N l is the column density of the inflowing gas; m X is the mass of the molecular ion; A obs is the de-magnified projected area covered by the inflowing gas in the sky; φ is the angle between the line of sight and the inflow direction; and V obs is the measured red-shifted velocity of the line. With these definitions, the term V obs / cos φ represents the de-projected inflow velocity V inflow , and the cross-sectional area of the cylinder perpendicular to the flow is given by A inflow = A obs / sin φ. The N l (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 × 10 14 ) cm −2 ( Table 3). The mass of the OH + molecular ion is m OH + = 2.824 × 10 −23 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 kpc; Table 6) yielding A obs = 0.13 − 0.17 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 of 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 and modulo the angle φ 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 × 10 −8 (computed as the Galactic average value from Table 5 of Indriolo et al. 2015) and f CH + = 7.6 × 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Ṁ inflow sin φ cos φ ∼ 3 − 14 M ⊙ yr −1 from OH + , and ∼ 6 − 24 M ⊙ 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: where Φ is the inflow rate; SFR = ǫ M 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; Ψ = λSFR is the outflow rate, linked to the SFR by the mass loading factor λ; finally, η rec is the fraction of gas ejected by outflows falling back onto the galaxy. Additional factors, e.g., 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, i.e., the term Φ 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 3× their SFRs. Turning into physical units,Ṁ inflow is a few solar masses per year at z = 0 growing to ∼ 100 M ⊙ yr −1 at z = 2.2. Model predictions reproduce these values for dark matter halos masses of M DH ∼ 2×10 12 M ⊙ (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 τ depl ) come into play. How they relate toṀ 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 ηǫ M gas in Eq. 19). In this case, an outflow should also be detected, as described by −λǫ M 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.

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 restframe, 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 + groundstate 1 J − 0 1 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 A&A proofs: manuscript no. berta_HerBS89a_v2p1 HCN(11−10) and/or NH(1 2 −0 1 ) lines, seen in absorption. Modelling 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 ×10 10 M ⊙ . Using lens modelling, we have reconstructed the sourceplane morphologies of HerBS-89a in dust continuum and molecular line emission. The magnification factor is estimated to be µ ∼ 5 for the dust continuum and µ ∼ 4 − 5 for the emission lines.
-Lens modelling suggests that the warmest and densest material in HerBS-89a (i.e., as traced by dust, 12 CO(9-8), and H 2 O(2 02 − 1 11 ) 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 ∼ 10 9 − 10 10 cm −3 ). NH 2 is seen in emission and detected in two transitions, while an absorption feature at 246.72 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 ≈ +100 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 ∼ 100 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 of the order of ≈ 10 − 20 M ⊙ 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 highquality spectral data in the 830 to 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 signal-to-noise ratio 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.   . Lens modelling 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 R.A. 13:16:11.52 and Dec. +28:12:17.7 (J2000).