© S. Berta et al. 2021

1. Introduction

In the last two decades, surveys in the far-infrared and submillimeter wavebands have opened up a new window for our understanding of the formation and evolution of galaxies, revealing a population of massive, dust-enshrouded galaxies forming stars at enormous rates in the early Universe (see, e.g., Blain et al. 2002; Carilli & Walter 2013; Casey et al. 2014; Hodge & da Cunha 2020). In particular, the extragalactic imaging surveys done with the Herschel Space Observatory (Pilbratt et al. 2010), such as Herschel-ATLAS (Eales et al. 2010a), HerMES (Oliver et al. 2012), and PEP (Lutz et al. 2011), have increased the number of dust-obscured star-forming galaxies from hundreds to several hundred thousand. Together with other large-area surveys, such as the all-sky Planck-HFI (Planck Collaboration Int. XXVII 2015), the South Pole Telescope (SPT; Carlstrom et al. 2011) cosmological survey (Staniszewski et al. 2009; Vieira et al. 2010), and the Atacama Cosmology Telescope (ACT; Marsden et al. 2014; Gralla et al. 2020), today we have vast samples of luminous dusty star-forming galaxies (DSFGs) that are amongst the brightest galaxies in the Universe, including numerous examples of strongly lensed systems (e.g., Negrello et al. 2010, 2017; Cox et al. 2011; Bussmann et al. 2013; Spilker et al. 2014; Cañameras et al. 2015; Reuter et al. 2020) and rare cases of galaxies with intrinsic infrared luminosities, L_FIR ⪆ 10¹³ L_⊙, and star formation rates (SFRs) in excess of 1000 M_⊙ yr⁻¹, known as Hyper-Luminous Infrared Galaxies (HyLIRGs, see, e.g., Ivison et al. 1998, 2013, 2019; Fu et al. 2013; Oteo et al. 2016; Riechers et al. 2013, 2017).

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

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

One of these galaxies, H-ATLAS J131611.5+281219 (hereafter HerBS-89a), at z = 2.9497, displays a very strong 2 mm continuum (with a flux density S_159 GHz = 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⁻¹. Both the 2 mm continuum and ¹²CO(5 − 4) line emission were resolved by the imaging, with an extension of in the east-west direction. The corresponding infrared luminosity of HerBS-89a (between 8 and 1000 μm in the rest-frame; uncorrected for amplification) is estimated to be (2.9 ± 0.2) × 10¹³ L_⊙, which suggested that HerBS-89a could be a HyLIRG.

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

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

The NOEMA data are complemented by Hubble Space Telescope (HST) imaging that traces the foreground massive lensing galaxy in the near-infrared (using the F110W filter around 1.1 μm), and by optical/near-infrared data obtained with the Gran Telescopio Canarias (GTC) that provide constraints on the redshift of the lensing galaxy. In addition, we also present data obtained with the Karl G. Jansky Very Large Array (VLA) on the ¹²CO(1 − 0) emission line, which allow us to derive the mass of molecular gas in HerBS-89a and give an anchoring point for an analysis of its CO spectral line energy distribution.

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

Throughout this paper we adopt a spatially flat ΛCDM cosmology with H₀ = 67.4 km s⁻¹ Mpc⁻¹ and Ω_M = 0.315 (Planck Collaboration VI 2020) and assume a Chabrier (2003) initial mass function (IMF). At the redshift z = 2.9497 of HerBS-89a, one arc-second corresponds to 7.9 kpc and the luminosity distance to the source is D_L = 2.5 × 10⁴ Mpc.

2. Observations and data reduction

2.1. NOEMA

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

The first project was a Discretionary Directorial Time (DDT) project, labeled E18AE (P.I.: R. Neri), observed on February 5, 2019 with ten antennas using the extended A-configuration, yielding an angular resolution of , for a total observing time of 4.3 h on-source time. This project, which was specifically tailored to measure the ¹²CO(9 − 8) and H₂O(2₀₂ − 1₁₁) emission lines, also enabled the detection of two strong absorption lines of the molecular ion OH⁺ corresponding to the redshifted ground state transitions of OH⁺(1₁ − 0₁) and OH⁺(1₂ − 0₁) (ν_rest = 1033.118 GHz and 971.803 GHz, respectively)¹.

The second project, labeled W19DE (P.I.: S. Berta), was completed on March 30th, 2020 with ten antennas using the intermediate C-configuration, yielding an angular resolution of , for a total on-source observing time of 4.2 h. This project was a follow-up of the DDT to measure the third ground state transition of OH⁺(1₀ − 0₁) (ν_rest = 909.158 GHz)¹ and the ground state transition of the molecular ion CH⁺(1 − 0) (ν_rest = 835.137 GHz)².

Observing conditions were excellent for both projects with an atmospheric phase stability of 20° rms and 0.8 mm of precipitable water vapor. The correlator was operated in the low resolution mode to provide spectral channels with nominal resolution of 2 MHz. The NOEMA antennas were equipped with 2SB receivers that cover a spectral window of 7.744 GHz in each sideband and polarization. For the first series of observations, we covered the frequency range from 244.9 to 252.6 GHz and 260.4 to 268.1 GHz; for the second series, the frequency ranges were from 209.4 to 217.1 GHz and 224.9 to 232.6 GHz.

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

2.2. VLA

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

2.3. HST

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

2.4. GTC

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

The second series of observations was performed on June 7, 2020 using the visitor instrument HiPERCAM for a total observing time of 1.25 h under average seeing conditions (). HiPERCAM is a quintuple-beam CCD imager enabling the five Sloan filter ugriz to be observed simultaneously over a field of view of with a pixel size of (Dhillon et al. 2018, and references therein). The HiPERCAM team designed a dedicated data reduction tool³ that applies standard procedures for the reduction and calibration of the raw images.

3. NOEMA and VLA results

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

3.1. NOEMA results

3.1.1. Continuum emission

Figure 1 presents the two continuum maps of HerBS-89a obtained by merging the lower and upper side-bands of each set of observations, resulting in a high-angular resolution image centered at ∼254.6 GHz and a lower angular resolution image at ∼220.8 GHz. We reach a sensitivity of σ = 29 and 33 μJy beam⁻¹ in the two bands, respectively. Merging the two side-bands of the high-resolution continuum images improves the S/Ns and resulted, by applying a uniform weight, in a final image with beam size of , to be compared to the achieved by the A-configuration that was used for these observations with natural weighting.

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

Figure 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.

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

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

To the east of HerBS-89a, the weak unresolved source (labeled HerBS-89b), which was already detected in 2 mm continuum by Neri et al. (2020), is detected in the 1.2 mm continuum emission in both images (Fig. 1), about 6″ east of the phase center of our observations. The detection reported here confirms the authenticity of this source. However, there is no corresponding source in the SDSS catalog at that position. The faintness of HerBS-89b precludes the extraction of spectroscopic information from the available NOEMA data.

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

3.1.2. Molecular emission and absorption lines

The large bandwidth of the NOEMA receivers allows us to search in HerBS-89a for redshifted emission and absorption lines over a wide range in frequency. The complete spectra of HerBS-89a, integrated over the areas subtended by the lines (see below) and normalized by the continuum, are displayed in Fig. 3. Each of the panels covers the frequency range of one sideband (LSB or USB) of one of the two observational projects discussed above, namely between 245 and 268 GHz (two upper panels) and between 209 and 233 GHz (two lower panels). Two strong molecular emission lines are present, ¹²CO(9 − 8) and H₂O(2₀₂ − 1₁₁), which both display wide profiles similar to those seen in the ¹²CO(3 − 2) and ¹²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₁ − 0₁) line, which is adjacent in frequency to the ¹²CO(9 − 8) emission line, and the (1₂ − 0₁) and (1₀ − 0₁) lines, which are detected for the first time in a high-z galaxy. Blue-shifted from the strong OH⁺(1₂ − 0₁) absorption line, another weaker absorption line is detected that is due to HCN(11 − 10) (ν_rest = 974.487 GHz) and/or NH(1₂ − 0₁) (ν_rest = 974.471 GHz) (see Sect. 8.1 for a detailed discussion). Next to the OH⁺(1₀ − 0₁) absorption line, another line is seen in emission, which corresponds to o − NH₂(2₀₂ − 1₁₁)(5/2 − 3/2) (ν_rest = 907.433 GHz). A second NH₂ line is detected with lower signal-to-noise ratio (S/N) at higher frequency, the NH₂(2₂₀ − 1₁₁)(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⁻¹) emission component is likely present, of which the extreme red and blue wings are detected at low S/N at either side of the absorption line.

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 rms of 0.6 and 0.4 mJy per channel, respectively. The rest-frame and observed frequencies are given on the upper and lower horizontal axes, respectively. All the detected molecular emission and absorption lines are identified (with the arrows positioned at the line frequencies). In addition, the redshifted positions of molecular lines (at z = 2.9497) that fall within the observed frequency range, but were not detected, are indicated (in italics): , H2O+(312 − 303), SH+(2 − 1), and o − NH2(20, 2 − 11, 1)(3/2 − 1/2).

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₀₂ − 1₁₁)(3/2 − 1/2) (ν_rest = 902.210 GHz), H₂O⁺(3₁₂ − 3₀₃) (ν_rest = 982.955 GHz), and (ν_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 ¹²CO(9 − 8) and H₂O(2₀₂ − 1₁₁) 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 ¹²CO(9 − 8) and H₂O(2₀₂ − 1₁₁) velocity averaged images, under the assumption that they are probing roughly the same reservoir of high density and temperature molecular gas in HerBS-89a. The increased S/N of the combined image (see upper right panel in Fig. 4) reveals the molecular emission to the west, joining the southern and northern peaks, and resembling the image of the dust continuum emission.

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. Upper right panel: combined 12CO(9 − 8) and H2O(202 − 111) image that enhances the low-level line emission, in particular along the western arc of the partial Einstein ring. Right panel in the second row: combined map of OH+ ((11 − 01) and (12 − 01)) absorption lines. The synthesized beam is shown in the lower left corner of each of the images.

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

The faint absorption lines of our lower angular resolution data are detected mainly toward the southern component of HerBS-89a; the CH⁺ absorption seems to be more extended than the OH⁺ absorption. The image of the HCN(11 − 10) and/or NH(1₂ − 0₁) 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₀₂ − 1₁₁)(5/2 − 3/2) emission line peaks westward of the continuum peak, as expected for red-shifted gas (see below).

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

Fig. 5. Spectra of the molecular emission and absorption lines detected in HerBS-89a reported in this study, together with the 12CO(5 − 4) and 12CO(3 − 2) spectra from Neri et al. (2020). Each molecular line is identified in the lower corner of each panel (labeled in black) and, in the case of a line overlap, as 12CO(9 − 8)/OH+(11 − 01), OH+(12 − 01)/HCN(11 − 10) and OH+(10 − 01)/o − NH2(202 − 111)(5/2 − 3/2), the second molecular line is labeled in red. The spectra are displayed with the continuum subtracted and, in each panel, the molecular line (labeled in black) is plotted relative to the zero velocity corresponding to its rest frequency. Fits to the emission and absorption molecular lines are shown as solid lines, and, when multi-component and/or multiple lines have been fitted, they are displayed individually as dashed lines (see text for details).

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⁻¹) and central velocities shifted by ±289 km s⁻¹ from the nominal line frequency at the redshift of HerBS-89a (z = 2.9497) (Table 2). We adopt here the same approach for the ¹²CO(9 − 8), H₂O(2₀₂ − 1₁₁) and NH₂(2₀₂ − 1₁₁)(5/2 − 3/2) emission lines, fixing their component velocities to the values found for ¹²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₂₀ − 1₁₁)(5/2 − 3/2) emission line, the detection has too low S/N to enable a reasonable fit and we therefore here report only upper limits to its line flux.

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

The OH⁺(1₂ − 0₁) 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⁻¹. Next to this OH⁺ absorption line, the broad weak absorption line, which corresponds to HCN(11 − 10) and/or NH(1₂ − 0₁), is centered at zero velocity and has a width of ∼800 km s⁻¹. 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₁) and NH₂(2₀₂ − 1₁₁)(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⁻¹). In the observed spectrum (Fig. 5) there is an apparent shift for both species to the red by ∼250 km s⁻¹ with respect to their nominal frequencies. For the OH⁺(1₀ − 0₁) absorption, this shift is larger than for the two other OH⁺ transitions, and in the case of NH₂, the profile seems to be inconsistent with those seen for other molecular emission lines tracing dense gas (CO and H₂O). Since NH₂ (amidogen) is a photodissociation product of NH₃ (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₂ emission line has been entirely suppressed by the OH⁺ absorption line. We have therefore fitted simultaneously the overlapping OH⁺ absorption and NH₂ emission lines using three Gaussian profiles, namely: (i) two components for NH₂, where the velocities are derived from the ¹²CO(5 − 4) and ¹²CO(3 − 2) emission lines and the ratio of the blue and red peaks similar to that of the double Gaussian profile of the ¹²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₂ − 0₁) 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.

3.1.3. 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

(1)

where S_ν is the flux density of the spectrum and that of the continuum only. The optical depth is evaluated at the central frequency ν₀ of the line. Table 3 reports the values of τ(ν₀) obtained from the integrated spectra shown in Fig. 3, that is, averaged over the angles covered by the lines on the sky. In the case of the two strongest absorption lines (the OH⁺(1₁ − 0₁) and (1₂ − 0₁) transitions), maps of the opacity are displayed in Fig. 6.

Fig. 6. Maps of the optical depth of the OH+(11 − 01) and OH+(12 − 01) absorption lines, evaluated as Fν/Fcont = −e−τν at the central frequency ν0 of each line.

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, 2018):

(2)

The quantity ∫τ_νdV and the speed of light c can be expressed in units of [cm s⁻¹], 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⁻²]. 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₂ − 0₁), listing the velocity integration ranges, the integrated optical depths, and the column densities. Uncertainties are computed from the dispersions of the spectra, via standard error propagation.

3.1.4. Velocity fields

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

Fig. 7. Velocity field maps (shown in colors) in each of the molecular emission and absorption lines detected in HerBS-89a superimposed on the continuum at the individual frequencies of each molecular line (in contours of 10σ spacing). 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. Upper right panel: combined 12CO(9 − 8) and H2O(202 − 111) image that reveals the kinematics along the western arm of the partial Einstein ring. Right panel in the second row: combined velocity field of the two OH+ absorption lines ((11 − 01) and (12 − 01)). The synthesized beam is shown in the lower left corner of each of the images.

3.2. VLA results: The ¹²CO(1–0) emission line

Figure 5 shows the VLA ¹²CO(1 − 0) emission line spectrum, resampled to a velocity resolution of ∼120 km s⁻¹ to highlight the key properties of the line profile. Because of the limited S/Ns, we fit the observed line with a single Gaussian profile. The position of the ¹²CO(1 − 0) emission line is consistent with the nominal redshift of HerBS-89a. Its FWHM is ΔV ∼ 1400 km s⁻¹, 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⁻¹ (Table 2).

The emission line map, which is shown in Fig. 8, was obtained by integrating the VLA data cube between −800 and +500 km s⁻¹ relative to the frequency of ¹²CO(1 − 0) at z = 2.9497. The ¹²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 ¹²CO(9 − 8) and (5 − 4) lines). It is noteworthy that, contrary to what is seen in the ¹²CO(9 − 8) and dust continuum emission, the northern component in ¹²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).

Fig. 8. 12CO(1 − 0) emission in HerBS-89a. The 1σ sensitivity is 32.1 μJy beam−1 and the beam size (shown in the lower left corner) is (PA −62.3°, E of N). The emission is resolved in the east-west direction. Contours are in 1σ steps starting from the 3σ level.

Due to the limited sensitivity of the VLA data, the underlying continuum of HerBS-89a was not detected. We measure an rms of ∼15.4 and ∼20.1 μJy beam⁻¹ in the continuum maps at 29 GHz (band width 0.74 GHz) and 38 GHz (band width 0.93 GHz) in the observed frame, with beam sizes of and , respectively.

4. The foreground lensing galaxy

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

4.1. Imaging results (HST)

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

Fig. 9. View of the foreground lensing galaxy in the HST WFC3 F110W filter (1.15 μm) image, with the 254.6 GHz dust continuum measured with NOEMA superposed as red contours.

4.2. Photometric redshift (GTC)

The foreground lensing galaxy is also detected in our GTC OSIRIS r-band imaging with a magnitude r  =  24.5 ± 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, 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 ( 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 a photometric redshift.

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

The best photometric redshift of the foreground lensing galaxy thus obtained is z_phot = 0.9, with a 99% confidence interval of 0.4 < z_phot, 3σ < 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 χ² distribution. The photometric redshift is best constrained by the r − i color, tracing the D4000 break⁴. The lack of rest-frame near-infrared coverage and the upper limit in the u band restrict the accuracy of the photo-z measurement. Deeper data and a wider wavelength coverage would be necessary to improve its precision.

Fig. 10. Best derivation of the photometric redshift of the foreground lensing galaxy, using the Le PHARE code (Ilbert et al. 2006; Arnouts et al. 1999). The filled triangle shows the HST F110W flux density and the filled circles represent the GTC measurements (see text). The WISE data (open circles) have not been included in the fit. The curve shows the best fit model, a Sb spiral template (Polletta et al. 2007; Ilbert et al. 2009). The inset displays the distribution of χ2 as a function of redshift; the dotted horizontal lines correspond to the 1, 2, 3σ levels, defined as . The best photometric redshift for the foreground galaxy is where the errors are 3σ.

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). However, none produced a better fit to the data.

4.3. Spectroscopic redshift (NOEMA)

An independent indication of the redshift of the lensing foreground can possibly be derived from the 2 mm spectral scan that was obtained to derive the spectroscopic resdhift of HerBS-89a in Neri et al. (2020). In the frequency range 139.3−147.0 GHz, where the ¹²CO(5 − 4) emission line of HerBS-89a was detected, a weaker narrow (∼200 km s⁻¹) 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 ¹²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.

Fig. 11. Spectrum between 139.29 and 147.04 GHz toward HerBS-89a, from Neri et al. (2020), showing the 12CO(5 − 4) emission line of HerBS-89a (at 145.902 GHz) and the weaker, narrow (∼200 km s−1) emission line at 143.74 GHz, possibly belonging to the foreground lensing galaxy (see Sect. 4.3 for further details).

5. Lens modeling

In order to recover the source-plane morphology and the intrinsic properties of HerBS-89a, we derive a lens model from the high-angular resolution NOEMA continuum data as shown in Fig. 2 (with initial lens parameters informed by HST observations) and proceed to apply it to the molecular line data. We perform the lens modeling and source reconstruction using the lensmodel package (Keeton 2001) with the pixsrc extension (Tagore & Keeton 2014; Tagore & Jackson 2016). For the sake of computational efficiency, all modeling has been done in the sky plane, although we recognize that fitting directly to visibility data can prevent specific choices in the interferometric imaging process (e.g., related to uv weighting and deconvolution) from having undue impact on a model. Although the quality of our NOEMA observations is excellent, the limited number of resolution elements across the source prevents us from leveraging the full power of the nonparametric approach of Tagore & Keeton (2014), instead limiting us to the use of parametric (Sérsic profile) components in reconstructing the sources.

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

Fig. 12. Lens modeling results for the 1.2 mm dust continuum emission and spectral lines in HerBS-89a (identified in the right most panels). In each row, from left to right: observed image (with the synthesized beam shown in the lower left corner) and critical curve (red line); reconstructed model image convolved with the PSF; residuals; reconstructed model image at full angular resolution; and (separated and with a different angular scale) reconstructed source-plane image and caustic. The residuals are within ±2.5σ, showing that the model agrees well with the data. The cross shows the coordinates of the phase center of NOEMA observations (RA 13:16:11.52 and Dec +28:12:17.7 in J2000). The OH+ absorption lines have been inverted to allow for a positive normalization of the Sérsic profile reconstruction.

The foreground deflector is assumed to have a singular isothermal elliptical mass distribution (SIEMD) described by five parameters: Einstein radius (b), RA, Dec, ellipticity (e)⁵, 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₀, ellipticity e, position angle θ_e, scale radius R_s, and Sérsic index n:

(3)

Using the estimated photometric redshift of the foreground galaxy and its uncertainty (see Sect. 4.2), we estimate that the deflector mass enclosed within the Einstein radius is M_⊙ (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

(4)

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₂O(2₀₂ − 1₁₁) zeroth-moment map, which shows faint (∼5σ, compared to the ∼10σ peak) emission to the west and northwest that motivated some exploration of a two-component model. We adopt here a single-component fit for this line, due to numerical instability in the two-component fit, but there are some systematic features in the residuals associated with the peak of the emission. The results are listed in Table 6, along with the corresponding magnification factors μ, for the high resolution data and for ¹²CO(1 − 0). Because of the lower resolutions and S/Ns, we do not model the spectral lines and continuum that were measured at lower frequencies (from 209 to 232 GHz).

The magnification factor μ 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 ¹²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 ¹²CO(1 − 0) emission would be helpful to better constrain the amplification factor as well as to study in greater detail the morphology (extent and kinematics) of the cold gas reservoir of HerBS-89a.

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

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

6.1. Continuum spectral energy distribution

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

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

Fig. 13. Spectral energy distribution. Top panel: observed SED of HerBS-89a (filled symbols) and 89b (open symbols). The HerBS-89a data are fit with three different dust models: DL07 (blue lines), optically-thin MBB (green lines); MBB in general form (red lines). The dotted lines represent the models as they are, and the dashed lines include the correction for the effects of CMB (see text for details). The light-blue continuous lines represent the synchrotron emission with a spectral index α = −0.8, normalized on the basis of the radio-FIR correlation subject to the redshift evolution found by Magnelli et al. (2015, upper line) and Delhaize et al. (2017, lower line). The data include SPIRE flux densities (black dots; from Bakx et al. 2018); SCUBA-2 photometry (blue dot; from Bakx et al. 2020); the NOEMA derived continuum flux densities at 1 mm continuum (green dots) from this work and 2 and 3 mm (red dots) from Neri et al. (2020); and 3σ upper limits at 29 and 38 GHz (10 and 7.8 mm, respectively) are derived from the VLA data (arrows). The data for HerBS-89b are taken from this paper for 1 mm and Neri et al. (2020) for 2 mm flux densities. Bottom panel: HerBS-89a NOEMA spectra (Neri et al. 2020, and this work) compared to the general form MBB model with CMB correction.

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

In the case of the DL07 modeling, interstellar dust is described as a mixture of carbonaceous and amorphous silicate grains, whose size distributions are chosen to reproduce different observed extinction laws. We defer to Draine & Li (2001, 2007), Li & Draine (2001), and Berta et al. (2016) for a detailed description of the model and its implementation. Despite the richness of the long-wavelength data, the SED is still lacking measurements in the mid-infrared, so that the fine details of the DL07 model, for example, the possible contribution of the 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:

(5)

where B_ν(T_dust) is the Planck function, T_dust is the dust temperature, and κ_ν = κ₀(ν/ν₀)^β 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 τ_ν = τ₀(ν/ν₀)^β, instead of the ν^β term. Here ν₀ is the frequency at which τ = τ₀ = 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, ν₀ is also a free parameter.

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

Berta et al. (2016) and Bianchi (2013) present thorough discussions of the proper use of κ_ν 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² kg⁻¹ 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).

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

The observed ¹²CO(1 − 0) emission line intensity can be translated into a luminosity (in K km s⁻¹ pc²), from which an estimate of the molecular gas mass can be made, using the following relation (Solomon & Vanden Bout 2005):

(6)

where S_LineΔV = I_Line is the integrated line flux in units of Jy km s⁻¹; ν_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:

(7)

We define the molecular and total gas masses as

(8)

(9)

where M_He is the mass of helium, M_H2 is the mass of molecular hydrogen and M_HI the mass of atomic hydrogen. Hereafter, we assume that M_mol ≫ M_HI and thus that M_gas ≈ M_mol.

The molecular hydrogen gas mass is computed from the ¹²CO(1 − 0) luminosity applying a conversion factor α_CO in units of M_⊙ (K km s⁻¹ pc²)⁻¹:

(10)

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 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. (2021) used the far-infrared continuum of a small but statistically complete sample of 12 low redshift galaxies together with CO and CI measurements to determine the CO − H₂ conversion factor α_CO. In another paper, Dunne et al. (in prep.) combined this sample with others measurements found in the literature for a total of more than 70 far-infrared and sub-mm selected galaxies with CO and CI detection and 230 with CO only. In these studies, a value of κ₈₅₀ = 0.065 m² kg⁻¹ is adopted (in agreement with the results by Planck Collaboration XXI 2011, see Appendix B for details) and it is found that a value α_CO ≃ 3 is consistent with the behavior of the three gas tracers in a diverse population of galaxies, ranging from high-z powerful star-forming galaxies to local disks.

For the computation of the molecular gas mass of HerBS-89a, we here adopt α_CO = 3.0 M_⊙ (K km s⁻¹ pc²)⁻¹, 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

(11)

(12)

The molecular gas mass derived from the ¹²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:

(13)

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):

(14)

where the rate of star formation is computed as SFR = 1.09 × 10⁻¹⁰ L_{8 − 1000 μm} with the infrared luminosity expressed in units of L_⊙ (Kennicutt 1998a, after modification for a Chabrier 2003 IMF).

Assuming that the observed infrared luminosity is entirely associated with the ongoing star formation activity of HerBS-89a, we derive an intrinsic star formation rate SFR = 614 ± 59 M_⊙ yr⁻¹ and a depletion timescale τ_depl = (3.4 ±  1.0) × 10⁸ 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 ¹²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⁻¹ pc²)⁻¹, 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⁻¹ 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⁷ 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.

6.3. ¹²CO spectral line energy distribution

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

Fig. 14. 12CO Spectral Line Energy Distribution (SLED) for HerBS-89a (red dots), compared to other well-studied local and high-z galaxies (or samples of galaxies). The CO SLEDs are normalized to the J = 1 − 0 transition, after correcting for the magnification, and the line fluxes of HerBS-89a have been corrected for amplification (see Sect. 5 and Table 6) before normalizing. The three samples of galaxies are from Bothwell et al. (2013) including 32 DSFGs at 1.2 < z < 4.1, Spilker et al. (2014), including 22 SPT-selected DFSGs at 2.0 < z < 5.7, and Valentino et al. (2020), including 78 “main sequence” star-forming galaxies. References to individual sources are: Cosmic Eyelash (Swinbank et al. 2011; Danielson et al. 2011); Cloverleaf (Barvainis et al. 1997; Weiß et al. 2003; Bradford et al. 2009; Riechers et al. 2011); GN20 (Carilli et al. 2010; Cortzen et al. 2020, and in prep.); M 82 (Weiß et al. 2005); and the Milky Way disk (Fixsen et al. 1999).

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_H2), the column density of CO per unit velocity gradient (N_CO/dV), and the solid angle (Ω_app) of the source. The CO SLED only depends on T_kin, n_H2 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 , the kinetic temperature K and the CO column density , where the errors are 3σ.

Fig. 15. Large Velocity Gradient (LVG) model fit to the 12CO SLED of HerBS-89a (after correcting for lens magnification). Top panel: data are shown as black dots and the solid curve shows the best fit from the single component (MCMC) model corresponding to the maximum posterior possibility (see Sect. 6.3 for details). Bottom panel: two-components fit, including an additional cold and extended component (darker, dot-dashed line).

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 cm⁻³, K, and cm⁻² km s⁻¹, where the errors are 3σ. The warm component is in broad agreement with the previous single component fitting results.

6.4. The para- H₂O(2₀₂–1₁₁) emission line

Water is one of the most abundant molecules after H₂ and CO in the gaseous interstellar medium (ISM). It serves as a unique diagnostic for probing the physical conditions of the ISM in both local (e.g., González-Alfonso et al. 2014) and high-redshift galaxies (e.g., Omont et al. 2013; Riechers et al. 2013; Yang et al. 2016, 2020, and references therein). In particular, since the E_upper > 200 K levels of H₂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₂O(2₀₂ − 1₁₁) emission line (ν_rest = 987.9268 GHz) in HerBS-89a is an illustration thereof: the line width is broad (ΔV ∼ 1100 km s⁻¹) and comparable to those of the mid/high-J CO lines, consistent with the finding that the J = 2 H₂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_H2O = 1.59 ± 0.37 Jy km s⁻¹, and comparable to the line flux of the ¹²CO(9 − 8) emission line (Table 2). The CO and H₂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_{H2O(202 − 111)} = (6.4 ± 1.5) × 10⁷ L_⊙, yielding a ratio L_{H2O(202 − 111)}/L_IR = (1.12 ± 0.36) × 10⁻⁵, in terms of the 8 to 1000 μm L_IR (Table 7). HerBS-89a lies on the correlation between L_IR and L_H2O, that is derived by Yang et al. (2016), as shown in Fig. 16. This result is in line with the fact that far-infrared pumping is the dominant mechanism for the excitation of the submillimeter H₂O emission that originates in very dense, heavily obscured star-formation dominated regions (see, e.g., González-Alfonso et al. 2014; Yang et al. 2016). The steeper-than-linear growth is likely the result of an increase of the optical depth in the dust continuum (at ∼100 μm) with increasing infrared luminosity. Finally, the ratio of L_{H2O(202 − 111)}/L_IR also indicates that an AGN component, if present in HerBS-89a, is not contributing much to the infrared luminosity, in contrast to a system like Mrk 231 or APM 08279+5455 where the AGN is dominant (and references therein Yang et al. 2016).

Fig. 16. Correlation between LIR and LH2O(202 − 111) for local (black symbols) and high-redshift (blue symbols) galaxies. HerBS-89a is labeled with a red star and is located along the tight relationship between the infrared and water luminosity (shown as a light blue line). Figure is adapted from Yang et al. (2016) and references therein.

7. Molecular gas kinematics

In an ideal situation, analysis of spectral line data cubes like those presented here for HerBS-89a can enable detailed analyses of source-plane kinematics. The lens model derived in Sect. 4 can in principle be applied to the emission in a particular line across many independent velocity channels, allowing them to be reconstructed in the source plane, used to compute moment maps, and fits to dynamical models (e.g., Sharon et al. 2019). In practice, although the quality of our NOEMA observations is excellent, the S/N and the number of resolution elements (i.e., synthesized beams) across the observed size of HerBS-89a are not high enough to fully achieve this goal. Several factors are at play here. First, as noted in Sect. 5 above, the combination of S/N and effective resolution limits us to the use of parametric (Sérsic profile) sources in reconstructing individual channels. Second, reconstruction of a given source-plane channel map tends to (a) favor sources that are very (sometimes implausibly) compact as a result of the use of a constrained emission profile, and (b) couple with the noise fluctuations in that channel in such a way that introduces uncertainties on the centroids of the reconstructed sources. Finally, comparison of a delensed source-plane reconstruction with an observed channel map necessarily involves the convolution of the former to the (coarser) resolution of the latter, meaning that many possible source-plane reconstructions are formally compatible with the observed data.

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

Fig. 17. Source plane reconstruction split by spectral channels. Left panel: positions (colored thick dots) and shapes (ellipses) of the reconstructed profiles of the individual 12CO(9 − 8) channels, over-plotted on the map of the dust continuum emission in the source plane. The ellipses are plotted at 1/20th of the peak flux of each reconstructed channel map. Central panel: zeroth moment map of the 12CO(9 − 8) emission line derived from the reconstruction of the seven individual channels shown in the left panel. The contours show the dust continuum emission, starting at 1/20 of the peak flux in steps of 1/20th (see also Fig. A.3). Right panel: zeroth moment map obtained splitting the 12CO(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.

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

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

Figure 17 provides a further illustration of the degree to which small noise fluctuations in the observed channel maps can perturb the centroids of the source-plane reconstructions. Initially inspired by the apparent bent morphology (and nominal associated velocity gradients; see below) in the central panel of Fig. 17, we constructed new, roughly half-line moment maps from the bluest four and reddest three of our seven 120 km s⁻¹ channels, and reconstructed these using the lens model optimized in Sect. 5. The sum of these two reconstructions (i.e., another zeroth moment map) is presented in the right panel of Fig. 17. The apparent double morphology might at first right suggest a pre-coalescence merging system; however, the ‘doubleness’ of the morphology stems at least in part from the general compactness of our Sérsic profile reconstructions. Moreover, notable aspects of the internal source structure (e.g., the distinct southeastern offset of the reconstructed −312 to −192 km s⁻¹ channel map) in the center panel do not translate to the reconstructions used to derive the right panel. This level of inconsistency raises the question of whether detailed information on the intrinsic or the internal source structure (as distinct from the overall extent of emission) can be robustly recovered from this dataset, and suggests that the apparent small-scale structure in the resulting zeroth-moment maps has to be interpreted with caution.

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

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

(15)

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

The individual dynamical masses of these components, assuming rotational support as above, would be M_dyn sin² i ≈ 2.6 × 10¹⁰ M_⊙ within R = 0.85 kpc for the eastern component and M_dyn sin² i ≈ 2.2 × 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.

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

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

8. Molecular lines other than CO and water

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

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

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

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

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

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

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

Between J = 5−4 and 11−10, HCN flips from emission to absorption as shown in the case of Arp 220 where Rangwala et al. (2011) detect the J = 12−11 to 17−16 transitions in absorption⁶. These transitions of HCN have critical densities around 10⁹ − 10¹⁰ cm⁻³ and probe very dense (n_H2 ≈ 3 × 10⁶ cm⁻³) 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.24 ± 0.12 yielding, in the case of HCN, a column density of (0.38 ± 0.1) × 10¹⁴ cm⁻² (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₂) to the formation and destruction of NH₃. It was first detected in the diffuse interstellar medium (Meyer & Roth 1991) and later in the dense molecular gas of SgrB2 (Cernicharo et al. 2000; Goicoechea et al. 2004). The Herschel SPIRE-FTS spectra of Mrk 231 and Arp 220 (González-Alfonso et al. 2018; Rangwala et al. 2011) show the presence of NH with various transitions detected. NH has a spectrum similar to that of OH⁺, since these molecules are isoelectronic. Many of the transitions of NH overlap with OH⁺ transitions, in particular NH(1₂ − 1₀) with OH⁺(1₀ − 0₁), with similar column densities derived for OH⁺ and NH in Mrk 231 (González-Alfonso et al. 2018). In the case of HerBS-89a, attributing all of the 246.72 GHz absorption line to NH implies a column density of (6.7 ± 1.9) × 10¹⁴ cm⁻², which is also comparable to the column density of OH⁺ (Table 3). Together with NH₂, the likely detection of NH in HerBS-89a opens up the possibility of exploring the chemistry of N-bearing molecules in a high-z galaxy by measuring other transitions of NH and searching for NH₃, and exploring whether (as is the case in SgrB2) the derived column densities are compatible with grain-surface chemistry and sputtering by shocks (see, e.g., Goicoechea et al. 2004).

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

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

In HerBS-89a, all three OH⁺ ground-state 1_J − 0₁ lines (J = 0, 1, 2) are seen in absorption (Fig. 5 middle panel). The (1₁ − 0₁) and (1₂ − 0₁) transitions are the strongest, while the (1₀ − 0₁) 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⁻¹ and the inferred OH⁺ column density is N_OH⁺ ∼ 10¹⁵ cm⁻² (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⁻¹ relative to the systemic velocity, including the (1₀ − 0₁) transition after the contamination of the NH₂ emission line is accounted for (see Sect. 3.1.2).

CH⁺. The methylidyne cation, CH⁺, has been shown to be a sensitive tracer of feedback mechanisms in high-z DSFGs (Falgarone et al. 2017), who report the detection of CH⁺(1 − 0) in six lensed galaxies at redshifts z ≈ 2.5, including the Cosmic Eyelash and SDP.17⁷. The CH⁺(1 − 0) line is found both in emission and absorption, providing critical information on the mechanisms of mechanical energy dissipation and/or strong UV irradiation in these high-z galaxies. Due to its high critical density for excitation (∼10⁷ cm⁻³ for the J = 1 − 0 transition), the CH⁺ line is seen in emission in dense (> 10⁵ cm⁻³) gas, probing shocked regions powered by galactic winds, and in absorption, in large (> 10 kpc) reservoirs of turbulent, cold (∼100 K), low-density (∼100 cm⁻³) 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⁻¹, comparable to that of the OH⁺ lines, and is also red-shifted to ≈ + 100 km s⁻¹ (Fig. 5). The inferred CH⁺ column density is estimated to be N_CH⁺ ∼ 10¹⁴ cm⁻² (Sect. 8), similar to what has been derived in the case of the Cosmic Eyelash (Falgarone et al. 2017).

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

9. Inflowing gas

In HerBS-89a, the three OH⁺ ground-state 1_J − 0₁ 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⁻¹. 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 ¹²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⁸ shows a red-shifted line, namely the Cosmic Eyelash (Falgarone et al. 2017); for this source, the CH⁺(1 − 0) absorption line is shifted by ΔV ∼ +150 km s⁻¹ with respect to the systemic velocity as defined by the CO emission lines. A similar velocity shift in the Cosmic Eyelash is also observed in the OH⁺(1₁ − 0₁) and H₂O⁺(1₁₁ − 0₀₀) absorption lines (Indriolo et al. 2018). Recently, Butler et al. (in prep.) report OH⁺(1₁ − 0₁) absorption lines toward ten gravitationally lensed Herschel-selected high-z galaxies; for all of them, the line is invariably seen in absorption and blue-shifted with respect to the systemic velocity, with the exception of one source where the absorption line peaks at the systemic velocity. Together, the currently available observations of OH⁺ (and CH⁺) have established that these molecular ions are robust tracers of molecular gas outflows in high-z galaxies, confirming the results of observations of OH⁺ or OH in local ultra-luminous galaxies (e.g., Lu et al. 2017; González-Alfonso et al. 2017, and references therein). In addition, recent studies of absorption lines in high-z lensed DSFGs in the OH 119 μm doublet, in SPT-2319−55 at z = 5.29 (Spilker et al. 2018) and the quasar ULASJ131911+095051 at z = 6.13 (Herrera-Camus et al. 2020), and in the H₂O(3₃₀ − 3₂₁) and (4₂₃ − 4₁₄) transitions in SPT-0346−52 at z = 5.656 (Jones et al. 2019), have revealed blue-shifted lines tracing powerful outflows (with rates of ∼200 − 500 M_⊙ yr⁻¹).

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₁ − 0₁) (or other tracers such as OH, H₂O or H₂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⁻¹, suggesting a low excitation inflow in this source (Falstad et al. 2017); Circinus, which shows an unambiguous OH 119 μm inverted P-Cygni profile (Stone et al. 2016); and IRAS 10091+4704, the starburst-dominated system at z = 0.24, were a clear inverted P-Cygni OH profile indicates the presence of a fast inflow of molecular gas with a rate of ∼100 M_⊙ yr⁻¹ (Herrera-Camus et al. 2020).

Theoretically, gas flowing into less massive (M ≤ 10¹² M_⊙) galaxies should be dynamically and thermally cold, while more massive halos receive most of their baryons as hotter (T > 10⁵ 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:

(16)

Fig. 19. Schematic of the simplified geometry of the inflow.

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₁ − 0₁) and (1₂ − 0₁) is (9.6 − 11.2 × 10¹⁴) cm⁻² (Table 3). The mass of the OH⁺ molecular ion is m_OH⁺ = 2.824 × 10⁻²³ 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². The observed velocity of the flow is 100 km s⁻¹. The distance that the inflowing gas travels from the observed position to the center of the galaxy is one of the most uncertain quantities in Eq. (16). Typical distances found in the literature are on the order of one to a few kpc. Using the positions of lines and continuum in the source plane (Table 6), we can only compute the distances projected on the sky as between ∼0.1 and ∼0.3 kpc, which are of the same order than the systematic positional uncertainty computed via Eq. (4) (∼0.08 kpc). Assuming a geometry similar to that of OH⁺, because no source plane reconstruction has been attempted for the low angular resolution NOEMA data, a similar reasoning can be applied to the CH⁺ line. Based on these values and estimates and using Eq. (16), we find

(17)

and

(18)

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⁻⁸ (computed as the Galactic average value from Table 5 of Indriolo et al. 2015) and f_CH⁺ = 7.6 × 10⁻⁹ (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⁻¹ from OH⁺, and ∼6 − 24 M_⊙ yr⁻¹ 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:

(19)

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, for example due to AGN feedback, possible mergers or other external causes, are here summarized by the term Y.

Different scenarios can be at play in the case of HerBS-89a. In the least likely scenario, we might be witnessing a flow channeling gas from the intergalactic medium onto the galaxy, described by the term Φ 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², 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⁻¹ at z = 2.2. Model predictions reproduce these values for dark matter halos masses of M_DH ∼ 2 × 10¹² 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.

10. Summary and concluding remarks

We report new millimeter NOEMA and VLA observations of the lensed starburst galaxy HerBS-89a at z = 2.95 that are complemented with near-infrared and optical data (obtained with the HST and GTC) that reveal the foreground lensing galaxy. The ability to process large instantaneous bandwidths with NOEMA enables an in-depth exploration of a sizeable part of the 1 mm spectrum of HerBS-89a, covering in total 121 GHz in the rest-frame, leading to the detection of transitions of molecular ions and molecular species that have not been previously observed in a high-z galaxy. In addition to the strong emission lines of ¹²CO(9 − 8) and H₂O(2₀₂ − 1₁₁), we detect the three OH⁺ ground-state 1_J − 0₁ 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₂ in emission, and HCN(11 − 10) and/or NH(1₂ − 0₁) lines, seen in absorption. modeling the lensing of HerBS-89a, we reconstructed the dust continuum emission and the molecular emission lines in the source plane deriving the morphology and the kinematics of this high-z galaxy.

The main conclusions of this paper are as follows:

• The high-angular resolution () images of HerBS-89a in dust continuum and the emission lines of ¹²CO(9 − 8) and H₂O(2₀₂ − 1₁₁) 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 diameter structure reminiscent of an Einstein ring. HerBS-89a is gravitationally amplified by a foreground massive galaxy at , with an estimated deflector mass of . Using lens modeling, we have reconstructed the source-plane morphologies of HerBS-89a in dust continuum and molecular line emission. The magnification factor is estimated to be μ ∼ 5 for the dust continuum and μ ∼ 4 − 5 for the emission lines.

• HerBS-89a is a high-z dusty star forming galaxy. From the ¹²CO(1 − 0) measurements, we derive a mass for the molecular gas reservoir of M_mol = (2.1 ± 0.4) × 10¹¹ M_⊙. modeling the available CO lines (J = 1 − 0, 3−2, 5−4 and 9−8) with an LVG approach indicates that the molecular gas in HerBS-89a has an average density of n_H2 ∼ 500 cm⁻³ and a kinetic temperature T_kin ∼ 200 K. Using the fine sampling of the dust continuum flux densities observed from 250 μm to 3.4 mm (or 63 to 860 μm in the rest-frame), we derive an infrared luminosity of L_IR = (4.6 ± 0.4) × 10¹² L_⊙, and a dust mass of M_dust = (2.6 ± 0.2) × 10⁹ M_⊙, yielding a dust-to-gas ratio δ_GDR ≈ 80. The star formation rate is estimated to be SFR = 614 ± 59 M_⊙ yr⁻¹ and the depletion timescale τ_depl = (3.4 ± 1.0) × 10⁸ years.

• Lens modeling suggests that the warmest and densest material in HerBS-89a (i.e., as traced by dust, ¹²CO(9–8), and H₂O(2₀₂ − 1₁₁) emission) is distributed across a region a few kpc across. We cannot distinguish definitively between a single-disk and a two-component description of this distribution on the basis of existing observations, but dynamical and molecular gas mass estimates are consistent in either scenario, and the two-component scenario in particular provides a natural explanation for the very large observed CO velocity widths.

• A series of molecular emission and absorption lines are detected in HerBS-89a that trace very dense gas (up to n_H2 ∼ 10⁹ − 10¹⁰ cm⁻³). NH₂ is seen in emission and detected in two transitions, while an absorption feature at 246.72 GHz corresponds to HCN(11–10), NH(1₂ − 1₀) 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⁻¹ with respect to the systemic velocity of HerBS-89a. The molecular ions OH⁺ and CH⁺ both trace low-density (n_H2 ∼ 100 cm⁻³) gas. The fact that the velocity red-shifted absorption lines cover substantially the distribution of the CO and water emission lines in the image (and source) planes indicates that this low-density gas must be kinematically unrelated to the dense gas. We here argue that we are tracing a rare case of inflow of low-density gas, indicating that HerBS-89a is accreting matter from its surroundings. Based on the available data, we have estimated that the mass inflow rate is on the order of ≈10 − 20 M_⊙ yr⁻¹. Various scenarios are possible to account for the observations including inflowing gas from the intergalactic medium, a galactic fountain or a stream of matter between the two merging components.

The results here described highlight the importance of processing large instantaneous bandwidth in the study of high-z galaxies, as such observations allow us to efficiently measure multiple lines and explore in an unbiased way the entire available frequency range for the presence of molecular species. In the case of HerBS-89a, we were able to probe this system from the very highest to the lowest molecular gas densities with high-quality spectral data in the 830−1060 GHz rest-frame frequency range, which includes key molecular ionic transitions tracing feedback mechanisms. The findings here reported show that it is possible to unravel the complex physics at play within this starburst (likely to be a merger) system in the early Universe from the kinematics of its inner most regions to the low-density inflowing molecular gas. They provide the foundation for future follow-up observations to further study the nature of HerBS-89a and better constrain its morphology and kinematics, which will need higher-angular resolution and higher S/N data, as well as for future studies of larger samples of high-z galaxies, with properties akin to those of HerBS-89a, selected from on-going large redshift surveys.

Acknowledgments

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

References

Appendix A: Details of lens model

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

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

Fig. A.2. Posterior distribution of the parameters of the HerBS-89a lensed source, as obtained by the lens modeling and MCMC runs (see Sect. 5 for details). The constraints are given on the following parameters: the position (ΔRA and ΔDec in arcsec), the Einstein radius (b), the ellipticity (e), the scale radius (Rs), the position angle (θe) with the orientation measured east of north and the Sérsic index (n). See Fig. A.1 for further details.

In order to perform the kinematic analysis described in Sect. 7, we apply the lens model to high resolution emission line data preserving the velocity information by fitting Sérsic profiles to channel images in the data cubes. To obtain a sufficient S/Ns, we bin the 40 km s⁻¹ channels of ¹²CO(9 − 8) and H₂O by a factor of 3 and 6, respectively. Despite this rebinning, the S/Ns for some of the channels is still very low.

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

Fig. A.3. Lens modeling of the 12CO(9 − 8) spectral channels. The velocity ranges covered by each channel are marked in the right-hand panels. From left to right: observed image (with the synthesized beam shown in the lower left corner) and critical curve (red line); reconstructed model image convolved with the PSF; residuals, reconstructed model image at full angular resolution; and, separated and with a different angular scale, reconstructed source-plane images and caustic curve. The cross shows the central coordinate RA 13:16:11.52 and Dec +28:12:17.7 (J2000).

Fig. A.4. Lens modeling of the H2O(202 − 111) spectral channels. The velocity ranges covered by each channel are marked in the right-hand panels. See Fig. A.3 for further details.

Appendix B: Review of κ_ν values

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

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

The tabulated κ_ν values mentioned above are theoretical (see, e.g., Weingartner & Draine 2001) and an alternative path is to adopt an empirical value of κ_ν, measured at a reference frequency and then scale it to any other frequency using the power law with index β as usual. Studying a sample of 22 local galaxies with known metallicities and available sub-mm data, James et al. (2002) derived κ₈₅₀ = 0.07 ± 0.02 m² kg⁻¹.

Using Herschel observations of local galaxies M 99 and M 100, Eales et al. (2010b) derive a value of κ_ν ∼ 0.06 m² kg⁻¹ at 350 μm, corresponding to κ₈₅₀ = 0.01 m² kg⁻¹ adopting a dust emissivity index β = 2. From early solar-neighborhood Planck results, a value κ₈₅₀ = 0.066 m² kg⁻¹ is derived⁹, roughly twice the commonly adopted value by Draine (2003). Recently, Clark et al. (2016) assembled a sample of 22 galaxies from the Herschel Reference Survey (Boselli et al. 2010) and derived a value of m² kg⁻¹ at 500 μm, corresponding to . Significant variations of κ_ν at sub-galactic scales are shown by Clark et al. (2019).

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

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

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

All Tables

All Figures

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

In the text

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

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 rms of 0.6 and 0.4 mJy per channel, respectively. The rest-frame and observed frequencies are given on the upper and lower horizontal axes, respectively. All the detected molecular emission and absorption lines are identified (with the arrows positioned at the line frequencies). In addition, the redshifted positions of molecular lines (at z = 2.9497) that fall within the observed frequency range, but were not detected, are indicated (in italics): , H2O+(312 − 303), SH+(2 − 1), and o − NH2(20, 2 − 11, 1)(3/2 − 1/2).

In the text

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. Upper right panel: combined 12CO(9 − 8) and H2O(202 − 111) image that enhances the low-level line emission, in particular along the western arc of the partial Einstein ring. Right panel in the second row: combined map of OH+ ((11 − 01) and (12 − 01)) absorption lines. The synthesized beam is shown in the lower left corner of each of the images.

In the text

Fig. 5. Spectra of the molecular emission and absorption lines detected in HerBS-89a reported in this study, together with the 12CO(5 − 4) and 12CO(3 − 2) spectra from Neri et al. (2020). Each molecular line is identified in the lower corner of each panel (labeled in black) and, in the case of a line overlap, as 12CO(9 − 8)/OH+(11 − 01), OH+(12 − 01)/HCN(11 − 10) and OH+(10 − 01)/o − NH2(202 − 111)(5/2 − 3/2), the second molecular line is labeled in red. The spectra are displayed with the continuum subtracted and, in each panel, the molecular line (labeled in black) is plotted relative to the zero velocity corresponding to its rest frequency. Fits to the emission and absorption molecular lines are shown as solid lines, and, when multi-component and/or multiple lines have been fitted, they are displayed individually as dashed lines (see text for details).

In the text

	Fig. 6. Maps of the optical depth of the OH+(11 − 01) and OH+(12 − 01) absorption lines, evaluated as Fν/Fcont = −e−τν at the central frequency ν0 of each line.
In the text

Fig. 7. Velocity field maps (shown in colors) in each of the molecular emission and absorption lines detected in HerBS-89a superimposed on the continuum at the individual frequencies of each molecular line (in contours of 10σ spacing). 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. Upper right panel: combined 12CO(9 − 8) and H2O(202 − 111) image that reveals the kinematics along the western arm of the partial Einstein ring. Right panel in the second row: combined velocity field of the two OH+ absorption lines ((11 − 01) and (12 − 01)). The synthesized beam is shown in the lower left corner of each of the images.

In the text

	Fig. 8. 12CO(1 − 0) emission in HerBS-89a. The 1σ sensitivity is 32.1 μJy beam−1 and the beam size (shown in the lower left corner) is (PA −62.3°, E of N). The emission is resolved in the east-west direction. Contours are in 1σ steps starting from the 3σ level.
In the text

	Fig. 9. View of the foreground lensing galaxy in the HST WFC3 F110W filter (1.15 μm) image, with the 254.6 GHz dust continuum measured with NOEMA superposed as red contours.
In the text

Fig. 10. Best derivation of the photometric redshift of the foreground lensing galaxy, using the Le PHARE code (Ilbert et al. 2006; Arnouts et al. 1999). The filled triangle shows the HST F110W flux density and the filled circles represent the GTC measurements (see text). The WISE data (open circles) have not been included in the fit. The curve shows the best fit model, a Sb spiral template (Polletta et al. 2007; Ilbert et al. 2009). The inset displays the distribution of χ2 as a function of redshift; the dotted horizontal lines correspond to the 1, 2, 3σ levels, defined as . The best photometric redshift for the foreground galaxy is where the errors are 3σ.

In the text

	Fig. 11. Spectrum between 139.29 and 147.04 GHz toward HerBS-89a, from Neri et al. (2020), showing the 12CO(5 − 4) emission line of HerBS-89a (at 145.902 GHz) and the weaker, narrow (∼200 km s−1) emission line at 143.74 GHz, possibly belonging to the foreground lensing galaxy (see Sect. 4.3 for further details).
In the text

Fig. 12. Lens modeling results for the 1.2 mm dust continuum emission and spectral lines in HerBS-89a (identified in the right most panels). In each row, from left to right: observed image (with the synthesized beam shown in the lower left corner) and critical curve (red line); reconstructed model image convolved with the PSF; residuals; reconstructed model image at full angular resolution; and (separated and with a different angular scale) reconstructed source-plane image and caustic. The residuals are within ±2.5σ, showing that the model agrees well with the data. The cross shows the coordinates of the phase center of NOEMA observations (RA 13:16:11.52 and Dec +28:12:17.7 in J2000). The OH+ absorption lines have been inverted to allow for a positive normalization of the Sérsic profile reconstruction.

In the text

Fig. 13. Spectral energy distribution. Top panel: observed SED of HerBS-89a (filled symbols) and 89b (open symbols). The HerBS-89a data are fit with three different dust models: DL07 (blue lines), optically-thin MBB (green lines); MBB in general form (red lines). The dotted lines represent the models as they are, and the dashed lines include the correction for the effects of CMB (see text for details). The light-blue continuous lines represent the synchrotron emission with a spectral index α = −0.8, normalized on the basis of the radio-FIR correlation subject to the redshift evolution found by Magnelli et al. (2015, upper line) and Delhaize et al. (2017, lower line). The data include SPIRE flux densities (black dots; from Bakx et al. 2018); SCUBA-2 photometry (blue dot; from Bakx et al. 2020); the NOEMA derived continuum flux densities at 1 mm continuum (green dots) from this work and 2 and 3 mm (red dots) from Neri et al. (2020); and 3σ upper limits at 29 and 38 GHz (10 and 7.8 mm, respectively) are derived from the VLA data (arrows). The data for HerBS-89b are taken from this paper for 1 mm and Neri et al. (2020) for 2 mm flux densities. Bottom panel: HerBS-89a NOEMA spectra (Neri et al. 2020, and this work) compared to the general form MBB model with CMB correction.

In the text

Fig. 14. 12CO Spectral Line Energy Distribution (SLED) for HerBS-89a (red dots), compared to other well-studied local and high-z galaxies (or samples of galaxies). The CO SLEDs are normalized to the J = 1 − 0 transition, after correcting for the magnification, and the line fluxes of HerBS-89a have been corrected for amplification (see Sect. 5 and Table 6) before normalizing. The three samples of galaxies are from Bothwell et al. (2013) including 32 DSFGs at 1.2 < z < 4.1, Spilker et al. (2014), including 22 SPT-selected DFSGs at 2.0 < z < 5.7, and Valentino et al. (2020), including 78 “main sequence” star-forming galaxies. References to individual sources are: Cosmic Eyelash (Swinbank et al. 2011; Danielson et al. 2011); Cloverleaf (Barvainis et al. 1997; Weiß et al. 2003; Bradford et al. 2009; Riechers et al. 2011); GN20 (Carilli et al. 2010; Cortzen et al. 2020, and in prep.); M 82 (Weiß et al. 2005); and the Milky Way disk (Fixsen et al. 1999).

In the text

Fig. 15. Large Velocity Gradient (LVG) model fit to the 12CO SLED of HerBS-89a (after correcting for lens magnification). Top panel: data are shown as black dots and the solid curve shows the best fit from the single component (MCMC) model corresponding to the maximum posterior possibility (see Sect. 6.3 for details). Bottom panel: two-components fit, including an additional cold and extended component (darker, dot-dashed line).

In the text

	Fig. 16. Correlation between LIR and LH2O(202 − 111) for local (black symbols) and high-redshift (blue symbols) galaxies. HerBS-89a is labeled with a red star and is located along the tight relationship between the infrared and water luminosity (shown as a light blue line). Figure is adapted from Yang et al. (2016) and references therein.
In the text

Fig. 17. Source plane reconstruction split by spectral channels. Left panel: positions (colored thick dots) and shapes (ellipses) of the reconstructed profiles of the individual 12CO(9 − 8) channels, over-plotted on the map of the dust continuum emission in the source plane. The ellipses are plotted at 1/20th of the peak flux of each reconstructed channel map. Central panel: zeroth moment map of the 12CO(9 − 8) emission line derived from the reconstruction of the seven individual channels shown in the left panel. The contours show the dust continuum emission, starting at 1/20 of the peak flux in steps of 1/20th (see also Fig. A.3). Right panel: zeroth moment map obtained splitting the 12CO(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.

In the text

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

In the text

	Fig. 19. Schematic of the simplified geometry of the inflow.
In the text

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

In the text

Fig. A.2. Posterior distribution of the parameters of the HerBS-89a lensed source, as obtained by the lens modeling and MCMC runs (see Sect. 5 for details). The constraints are given on the following parameters: the position (ΔRA and ΔDec in arcsec), the Einstein radius (b), the ellipticity (e), the scale radius (Rs), the position angle (θe) with the orientation measured east of north and the Sérsic index (n). See Fig. A.1 for further details.

In the text

Fig. A.3. Lens modeling of the 12CO(9 − 8) spectral channels. The velocity ranges covered by each channel are marked in the right-hand panels. From left to right: observed image (with the synthesized beam shown in the lower left corner) and critical curve (red line); reconstructed model image convolved with the PSF; residuals, reconstructed model image at full angular resolution; and, separated and with a different angular scale, reconstructed source-plane images and caustic curve. The cross shows the central coordinate RA 13:16:11.52 and Dec +28:12:17.7 (J2000).

In the text

	Fig. A.4. Lens modeling of the H2O(202 − 111) spectral channels. The velocity ranges covered by each channel are marked in the right-hand panels. See Fig. A.3 for further details.
In the text