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. The reduced datacubes are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http: //cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/646/A122 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. 2010Negrello et al. , 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 13 L , and star formation rates (SFRs) in excess of 1000 M yr -1 , known as Hyper-Luminous Infrared Galaxies (HyLIRGs, see, e.g., Ivison et al. 1998Ivison et al. , 2013Ivison et al. , 2019;;Fu et al. 2013;Oteo et al. 2016;Riechers et al. 2013Riechers et al. , 2017)).
Detailed follow-up studies of these galaxies to investigate their nature and physical properties require robust estimates of their distances. Due to the dust obscuration in these objects, searching for CO emission lines via 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 -foot_0 . Based on their estimated locations relative to the L CO(1-0) versus ∆V(CO) relationship of Harris et al. (2012), several of the sources are inferred to be gravitationally amplified, while a number of them appear to belong to the rare class of hyper-luminous infrared galaxies (Neri et al. 2020).
One of these galaxies, H-ATLAS J131611.5+281219 (hereafter HerBS-89a), at z = 2.9497, displays a very strong 2 mm continuum (with a flux density S 159 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 -1 . Both the 2 mm continuum and 12 CO(5-4) line emission were resolved by the 1. 2 imaging, with an extension of 0. 9 ± 0. 1 in the eastwest direction. The corresponding infrared luminosity of HerBS-89a (between 8 and 1000 µm in the rest-frame; uncorrected for amplification) is estimated to be (2.9 ± 0.2) × 10 13 L , which suggested that HerBS-89a could be a HyLIRG.
In 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 12 CO(9-8) and H 2 O(2 02 -1 11 ); the three lines of the ground state of the molecular ion OH + , (1 1 -0 1 ), (1 2 -0 1 ) and (1 0 -0 1 ), 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 (2 02 -1 11 )(5/2-3/2) seen in emission; and and/or NH(1 2 -0 1 ) seen in absorption, here also reported for the first time in a high-z galaxy.
In HerBS-89a, the 0. 3 images of the molecular emission lines and the dust continuum reveal a partial 1. 0 diameter Einstein ring. All of the observed transitions of the molecular ions OH + and CH + are seen in absorption against the strong dust continuum. Together, the molecular emission and absorption lines allow us to probe simultaneously in HerBS-89a regions that have very different properties, ranging from the (very) dense molecular gas traced by 12 CO(9-8), H 2 O (2 02 -1 11 ), NH 2 and HCN(11-10) (overlapping with NH(1 2 -0 1 )), to the low H 2 fraction, diffuse gas traced in OH + and the reservoirs of turbulent, cold and low-density molecular gas traced by CH + , thereby providing a unique view at sub-kpc spatial resolution (in the source plane) of the physical properties and feedback activity in this high-z system.
The NOEMA data are complemented by Hubble Space Telescope (HST) imaging that traces the foreground massive lensing galaxy in the near-infrared (using the F110W filter around 1.1 µm), and by optical/near-infrared data obtained with the Gran Telescopio Canarias (GTC) that provide constraints on the redshift of the lensing galaxy. In addition, we also present data obtained with the Karl G. Jansky Very Large Array (VLA) on the 12 CO(1-0) emission line, which allow us to derive the mass of molecular gas in HerBS-89a and give an anchoring point for an analysis of its CO spectral line energy distribution.
The structure of the paper is as follows. Section 2 describes the NOEMA, VLA, HST, and GTC observations, and the reduction of the respective data sets; 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 2 O are outlined in Sect. 8; and a discussion of the gas inflow suggested by the red-shifted molecular absorption lines of OH + and CH + is presented in Sect. 9. Finally, 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 0 = 67.4 km s -1 Mpc -1 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 4 Mpc.
We used NOEMA to target high-frequency molecular lines in HerBS-89a, redshifted into the 1 mm band. The observations were carried out under two separate projects.
The first project was a Discretionary Directorial Time (DDT) project, labeled E18AE (P.I.: R. Neri), observed on February 5, 2019 with ten antennas using the extended A-configuration, yielding an angular resolution of 0. 3, for a total observing time of 4.3 h on-source time. This project, which was specifically tailored to measure the 12 CO(9-8) and H 2 O(2 02 -1 11 ) emission lines, also enabled the detection of two strong absorption lines of the molecular ion OH + corresponding to the redshifted ground state transitions of OH + (1 1 -0 1 ) and OH + (1 2 -0 1 ) (ν rest = 1033.118 GHz and 971.803 GHz, respectively) 1 .
The second project, labeled W19DE (P.I.: S. Berta), was completed on March 30th, 2020 with ten antennas using the intermediate C-configuration, yielding an angular resolution of 1. 0 × 0. 6, for a total on-source observing time of 4.2 h. This project was a follow-up of the DDT to measure the third ground state transition of OH + (1 0 -0 1 ) (ν rest = 909.158 GHz) 1 and the ground state transition of the molecular ion CH + (1-0) (ν rest = 835.137 GHz) 2 .
Observing conditions were excellent for both projects with an atmospheric phase stability of 20 • 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 <0. 2.
The National Radio Astronomy Observatory's (NRAO) VLA was used to observe the 12 CO(1-0) emission line in HerBS-89a when the array was in the C configuration. The observations are part of a larger project (program I.D.:VLA/20A-083 -P.I.: D. Riechers) whose goal was to measure the 12 CO(1-0) emission lines in the sample of Herschel-selected galaxies studied in Neri et al. (2020); the complete results of that project will be presented in a forthcoming paper (Stanley et al., in prep.). The data were acquired on March 22, 2020 under stable atmospheric conditions. The 0.9 cm Ka band receivers were tuned (1 GHz bandwidth per IF) to the expected frequency of the CO emission line, that is 29.1848 GHz based on the redshift determined by Neri et al. (2020), and to 38.499 GHz for the second baseband. In total, we observed for 1.6 h (with 62 min on-source) recording 1024×2 MHz dual-polarization channels across a total bandwidth of 2 GHz, which was chosen to maximize line sensitivity while retaining as much bandwidth as possible. The 2 GHz bandwidth setup was used to maximize the potential for stacking of faint lines (for the entire VLA project), while at same time retaining sufficient spectral resolution (2 MHz, i.e., 16-21 km s -1 ) to finely sample the broad CO emission line of HerBS-89a. The 8bit samplers were selected to maximize sensitivity. The source 3C286 was observed to determine accurate complex gain and bandpass correction solutions and was also observed to set the absolute scale flux density based on the Perley & Butler (2017) models; the pointing accuracy was checked regularly. The data were calibrated, averaged in polarization, mapped using natural baseline weighting and analyzed in the CASA (Common Astronomy Software Applications) package. The resulting line map has a spatial resolution of 1. 24 × 0. 79 (PA -62.3 • ) and a rms noise level of 32.1 µJy beam -1 over a band width of 0.13 GHz. The absolute flux scale is estimated to be accurate to within 10%.
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 .
HerBS-89b quality. The pixel scale of the Infrared-Camera is 0. 128, but the image was resampled to a finer pixel scale of 0. 064 by exploiting the adopted sub-pixel dither pattern. The astrometry of the resulting image was calibrated by matching the positions of 11 stars from SDSS DR12 and 2MASS and is accurate to within 0. 1.
HerBS-89a was observed with the GTC 10.4 m telescope in two observing runs (program ID: GTC03-19ADDT and GTC09-19ADDT; PI: H. Dannerbauer). First, optical imaging was obtained using the instrument OSIRIS on April 6, 2019.
The observations were conducted in service mode under clear skies (although with non-photometric conditions), integrating for 10 min in the Sloan r-band filter with a seeing of 0. 8. The field of view is 7 .5 × 6 .0 and the pixel size 0. 254. Standard procedures for data reduction and calibration of the raw images were performed with the IRAF data reduction package (Tody 1986).
The second series of observations was performed on June 7, 2020 using the visitor instrument HiPERCAM for a total observing time of 1.25 h under average seeing conditions (0. 9). HiPERCAM is a quintuple-beam CCD imager enabling the five A122, page 3 of 33 Sloan filter ugriz to be observed simultaneously over a field of view of 2 .8 × 1 .4 with a pixel size of 0. 081 (Dhillon et al. 2018, and references therein). The HiPERCAM team designed a dedicated data reduction toolfoot_2 that applies standard procedures for the reduction and calibration of the raw images.
In this section, we describe the new data obtained on HerBS-89a, outline the properties derived from the NOEMA 1 mm highangular resolution observations for both the continuum and the molecular emission and absorption lines (Sect. 3.1), and present the results on the 12 CO(1-0) emission line measured with the VLA (Sect. 3.2).
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 -1 in the two bands, respectively. Merging the two side-bands of the high-resolution continuum images improves the S/Ns and resulted, by applying a uniform weight, in a final image with beam size of 0. 43 × 0. 22, to be compared to the 0. 55 × 0. 33 achieved by the Aconfiguration that was used for these observations with natural weighting.
Figure 2 presents a zoom-in on this high-resolution continuum image of HerBS-89a. It displays an Einstein ring-like morphology, with a double source linked by weak arc structures, indicating that HerBS-89a is lensed. The high quality of the image reveals the details of the morphology of the dust continuum, such as the differences between the northern and southern continuum peaks as well as the weaker emission extending between them.
For the lower frequency data obtained using the mediumcompact C-configuration of NOEMA, the source is only marginally resolved. Combining the data from the upper and lower side-bands resulted in a beam of 1. 1 × 0. 66.
To the east of HerBS-89a, the weak unresolved source (labeled HerBS-89b), which was already detected in 2 mm continuum by Neri et al. (2020), is detected in the 1.2 mm continuum emission in both images (Fig. 1), about 6 east of the phase center of our observations. The detection reported here confirms the authenticity of this source. However, there is no corresponding source in the SDSS 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).
The large bandwidth of the NOEMA receivers allows us to search in HerBS-89a for redshifted emission and absorption lines over a wide range in frequency. The complete spectra of HerBS-89a, integrated over the areas subtended by the lines (see below) and normalized by the continuum, are displayed in Fig. 3. Each of the panels covers the frequency range of one sideband (LSB or USB) of one of the two observational projects discussed above, namely between 245 and 268 GHz (two upper panels) and between 209 and 233 GHz (two lower panels). Two strong molecular emission lines are present, 12 CO(9-8) and H 2 O(2 02 -1 11 ), which both display wide profiles similar to those seen in the 12 CO(3-2) and 12 CO(5-4) emission lines (Neri et al. 2020). In addition, the three lines of the ground state of the molecular ion OH + are all detected in absorption: the (1 1 -0 1 ) line, which is adjacent in frequency to the 12 CO(9-8) emission line, and the (1 2 -0 1 ) and (1 0 -0 1 ) lines, which are detected for the first time in a high-z galaxy. Blueshifted from the strong OH + (1 2 -0 1 ) absorption line, another weaker absorption line is detected that is due to HCN(11-10) (ν rest = 974.487 GHz) and/or NH(1 2 -0 1 ) (ν rest = 974.471 GHz) (see Sect. 8.1 for a detailed discussion). Next to the OH + (1 0 -0 1 ) absorption line, another line is seen in emission, which corresponds to o -NH 2 (2 02 -1 11 )(5/2-3/2) (ν rest = 907.433 GHz). A second NH 2 line is detected with lower signal-to-noise ratio (S/N) at higher frequency, the NH 2 (2 20 -1 11 )(5/2-3/2) transition at ν rest = 993.3226 GHz. Finally, the molecular ion CH + is detected in the ground transition, CH + (1-0); its profile is dominated by an absorption line that has a width similar to those of the OH + absorption lines and at the same red-shifted velocity; in addition, a weak and broad (∼850 km s -1 ) emission component is likely present, of which the extreme red and blue wings are detected at low S/N at either side of the absorption line.
Also identified in Fig. 3 are the positions of red-shifted molecular lines that fall within the observed frequency range but remained undetected in HerBS-89a, including: SH + (2-1) (ν rest = 893.152 GHz), o -NH 2 (2 02 -1 11 )(3/2-1/2) (ν rest = 902.210 GHz), H 2 O + (3 12 -3 03 ) (ν rest = 982.955 GHz), and o -H 3 O + (0 - 0 -1 + 0 ) (ν rest = 984.708 GHz), most of which are seen in the spectrum of the local merger Arp 220 (Rangwala et al. 2011).
Figure 4 shows the angular extents of the eight molecular emission and absorption lines detected with NOEMA in HerBS-89a, each compared to the continuum extracted in the respective sideband to which it belongs. The 12 CO(9 -8) andH 2 O(2 02 -1 11 ) emission lines are detected in both lensed components of HerBS-89a, but are slightly shifted in position with respect to each other, with the water line being well centered on the continuum emission peaks. In order to recover more information on the distribution of the molecular emission lines, we combined the 12 CO(9-8) and H 2 O(2 02 -1 11 ) velocity averaged images, under the assumption that they are probing roughly the same reservoir of high density and temperature molecular gas in HerBS-89a. The increased S/N of the combined image (see upper right panel in Fig. 4) reveals the molecular emission to the west, joining the southern and northern peaks, and resembling the image of the dust continuum emission.
The two strong OH + absorption lines detected in our higher angular resolution data, namely the (1 1 -0 1 ) and (1 2 -0 1 ) transitions, also coincide with the southern and northern peaks of the continuum emission. The combined OH + (1 1 -0 1 ) and (1 2 -0 1 ) map, integrated between -400 and +400 km s -1 from the line center in order to avoid contamination by neighboring lines (right panel in the second row of Fig. 4), further highlights this spatial coincidence.
The faint absorption lines of our lower angular resolution data are detected mainly toward the southern component of HerBS-89a; the CH + absorption seems to be more extended than the OH + absorption. The image of the HCN(11-10) and/or NH(1 2 -0 1 ) absorption line shows two peaks exactly centered on the dust continuum emission peaks with no indication of emission extending further out. Finally, the NH 2 (2 02 -1 11 )(5/2-3/2) emission line peaks westward of the continuum peak, as expected for red-shifted gas (see below).
In Fig. 5, we show the integrated spectra of all the individual molecular emission and absorption lines detected in HerBS-89a, including the 12 CO(3-2) and 12 CO(5-4) emission lines reported in Neri et al. (2020). In the following, we discuss the properties of each of these molecular lines.
For the two low-J CO emission lines, Neri et al. (2020) fitted each of the double-peaked line profiles with two Gaussian components, with identical widths (FWHM = 631 km s -1 ) and central velocities shifted by ±289 km s -1 from the nominal line frequency at the redshift of HerBS-89a (z = 2.9497) (Table 2). We adopt here the same approach for the 12 CO(9-8), H 2 O(2 02 -1 11 ) and NH 2 (2 02 -1 11 )(5/2-3/2) emission lines, fixing their component velocities to the values found for 12 CO(5-4) and (3-2). We also force the two Gaussian components (of each line) to have the same width (see Table 2). For the NH 2 (2 20 -1 11 )(5/2-3/2) emission line, the detection has too low S/N to enable a reasonable fit and we therefore here report only upper limits to its line flux.
As mentioned above, the CH + (1-0) line seems to display weak emission features on both sides of the absorption line. The fit to the profile shows the overlapping deep absorption feature and the less-well defined, broad component seen in emission.
The OH + (1 2 -0 1 ) absorption line displays strong, red-shifted absorption, whose profile is well reproduced by a narrow and deep absorption line, red-shifted by ∼100 km s -1 . Next to this OH + absorption line, the broad weak absorption line, which corresponds to HCN(11-10) and/or NH(1 2 -0 1 ), is centered at zero velocity and has a width of ∼800 km s -1 . As shown in Fig. 4, this broad absorption line is very concentrated spatially, peaking on the northern and southern peaks of the dust continuum emission.
The OH + (1 0 -0 1 ) and NH 2 (2 02 -1 11 )(5/2-3/2) lines are close to each other in frequency (i.e., 230.18 and 229.74 GHz, corresponding to a velocity difference of ∆V ∼ 150 km s -1 ). In the observed spectrum (Fig. 5) there is an apparent shift for both species to the red by ∼250 km s -1 with respect to their nominal 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): o -
frequencies. For the OH + (1 0 -0 1 ) absorption, this shift is larger than for the two other OH + transitions, and in the case of NH 2 , the profile seems to be inconsistent with those seen for other molecular emission lines tracing dense gas (CO and H 2 O). Since NH 2 (amidogen) is a photodissociation product of NH 3 (ammonia), it should probe the dense gas traced in the high-J CO or water lines, and therefore display a similar double Gaussian profile. However, we only see the red-shifted component of the double-peaked line profile, most likely because the blue-shifted part of the NH 2 emission line has been entirely suppressed by the OH + absorption line. We have therefore fitted simultaneously the overlapping OH + absorption and NH 2 emission lines using three Gaussian profiles, namely: (i) two components for NH 2 , where the velocities are derived from the 12 CO(5-4) and 12 CO(3-2) emission lines and the ratio of the blue and red peaks similar to that of the double Gaussian profile of the 12 CO(5-4) emis-sion line; (ii) one Gaussian for the OH + absorption line whose width and the center velocity are derived from the narrow and deep OH + (1 2 -0 1 ) absorption line. The fit shown in Fig. 5 reproduces very well the observed spectrum including the velocity shift described above. Table 2 summarizes the results of the spectral line fitting. For each molecular line, the main properties of the best Gaussian fits are listed, including the central velocity and frequency, the full Notes. The data for the 12 CO(3-2) and 12 CO(5-4) emission lines are from Neri et al. (2020). The fitting of each line is shown in Fig. 5. When one parameter is quoted without uncertainty, its value was fixed during the fitting procedure. The observed line luminosities µL line have been derived using the standard relation given by Solomon & Vanden Bout (2005) -see Sect. 6.2 and Eq. ( 6). ( * ) Combination of the two Gaussian components. ( * * ) Only the visible red component. ( * * * ) Net absorption only. ( †) The two NH 2 transitions are (2 02 -1 11 ) (5/2-3/2) and (2 20 -2 11 ) (5/2-3/2), respectively. In the case of the latter transition, we provide a 3σ upper limit. ( ‡) In the case of the absorption line observed at 246.72 GHz, we here provide the derived properties for HCN and NH separately, noting that both molecules likely do contribute to the observed absorption (see Sect. 8.1).
Fig. 5. Spectra of the molecular emission and absorption lines detected in HerBS-89a reported in this study, together with the 12 CO(5-4) and 12 CO(3-2) spectra from Neri et al. (2020). Each molecular line is identified in the lower corner of each panel (labeled in black) and, in the case of a line overlap, as 12 CO(9-8)/OH + (1 1 -0 1 ), OH + (1 2 -0 1 )/HCN(11-10) and OH + (1 0 -0 1 )/o -NH 2 (2 02 -1 11 )(5/2-3/2), the second molecular line is 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). width at half maximum (FWHM), the integrated intensity, and the intrinsic frequency. All the uncertainties are propagated into the derived quantities.
The optical depth of an absorption line (τ ν ) can be evaluated from the observed spectrum normalized by the continuum as follows
where S ν is the flux density of the spectrum and S cont ν that of the continuum only. The optical depth is evaluated at the central frequency ν 0 of the line. Table 3 reports the values of τ(ν 0 ) obtained from the integrated spectra shown in Fig. 3, 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 1 -0 1 ) and (1 2 -0 1 ) transitions), maps of the opacity are displayed in Fig. 6.
Using Eq. ( 1) and the continuum-normalized spectra, we compute "optical depth spectra" that describe the variation of τ ν along the absorption lines, as a function of frequency or velocity. The column density N l of the absorption line is derived by integrating over the velocity range of the line along the line of sight (e.g., Comito et al. 2003; see also Indriolo et al. 2015Indriolo et al. , 2018)):
The quantity τ ν dV and the speed of light c can be expressed in units of [cm s -1 ], and the transition frequency ν in [Hz] (see Table 2 for the corresponding frequencies of each line), so that N l is directly computed in units of [cm -2 ]. The quantities g l and g u are the lower and upper state statistical weights, and A ul is the spontaneous emission coefficient. We compute the column density for the strongest of the hyperfine transitions for the specific rotational ∆J, and we assume that nearly all molecules are in the ground rotational state. The values of the statistical weights and the emission coefficients are taken from the Cologne Database for Molecular Spectroscopy (CDMS; Müller et al. 2005). Table 3 summarizes the results for the absorption lines observed in HerBS-89a, i.e, the three OH + transitions, Table 3. Optical depths and column densities for the molecular absorption lines.
Line
0.36 ± 0.13 -145 to +500 91.4 ± 22.8 2.6 ± 0.7 HCN(11-10) 0.24 ± 0.12 -570 to +430 66.1 ± 18.3 0.38 ± 0.1 NH(1 2 -0 1 ) 0.24 ± 0.12 -570 to +430 66.1 ± 18.3 6.7 ± 1.9
Notes. The frequency ν 0 corresponds to the deepest dip of each lines. In the case of the absorption line observed at 246.72 GHz, we provide two possible column densities obtained assuming that the line is a transition of either HCN or NH separately, although we note that both lines likely contribute at the same time (see Sect. 8.1).
CH + (1-0) and the absorption line observed at 246.72 GHz that is due to HCN(11-10) and/or NH(1 2 -0 1 ), listing the velocity integration ranges, the integrated optical depths, and the column densities. Uncertainties are computed from the dispersions of the spectra, via standard error propagation.
Detailed views on the dynamics of HerBS-89a are presented in Fig. 7, which displays the velocity fields of all the molecular emission and absorption lines detected in this source. The emission lines of 12 CO(9-8) and H 2 O(2 02 -1 11 ) show similar, regular east-west velocity gradients, from blue to red along the southern arc, and reversed in the northern peak, as expected from gravitational lensing (see Sect. 5). The combined image of the 12 CO(9-8) and H 2 O(2 02 -1 11 ) emission lines reveals further details, in particular the velocity distribution along the western arc. The velocity field of these dense gas tracers is consistent with the presence of kinematically distinct components that could indicate a merger system, or with a single rotating 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 (2 02 -1 11 )(5/2-3/2) emission line is clearly seen west of the continuum peak, as expected from the velocity field revealed in the higher angular resolution data for the 12 CO(9-8) and H 2 O(2 02 -1 11 ) emission lines.
The 12 CO(1-0) emission line
Figure 5 shows the VLA 12 CO(1 -0) emission line spectrum, resampled to a velocity resolution of ∼120 km s -1 to highlight the key properties of the line profile. Because of the limited S/Ns, we fit the observed line with a single Gaussian profile. The position of the 12 CO(1-0) emission line is consistent with the nominal redshift of HerBS-89a. Its FWHM is ∆V ∼ 1400 km s -1 , comparable to those of the three high-J CO emission lines, and the observed total line intensity is I ν = 0.64 ± 0.13 Jy km s -1 (Table 2). The emission line map, which is shown in Fig. 8, was obtained by integrating the VLA data cube between -800 and +500 km s -1 relative to the frequency of 12 CO(1-0) at z = 2.9497. The 12 CO(1-0) line emission is slightly resolved and displays an east-west extension, similar to what is observed in the other CO lines for which the emission has been resolved (i.e., the 12 CO(9-8) and (5-4) lines). It is noteworthy that, contrary to what is seen in the 12 CO(9-8) and dust continuum emission, the northern component in 12 CO(1-0) appears slightly stronger than the southern component; however, the difference is at a <1σ level. If confirmed, this inconsistency may be the result of differential magnification and indicate that the molecular gas reservoir traced in CO (1-0) is more extended than that traced by the high-J CO and water transitions (see Sect. 5 for a detailed discussion).
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 -1 in the continuum maps at 29 GHz (band width 0.74 GHz) and 38 GHz (band width 0.93 GHz) in the observed frame, with beam sizes of 1. 23×0. 79 and 0. 96 × 0. 59, respectively.
In this section, we derive the characteristics of the foreground lensing galaxy from the near-infrared HST image and the GTC photometry, supplemented by NOEMA data.
Figure 9 displays the HST wide-J filter F110W image (at 1.15 µm) toward HerBS-89a compared to the 1.2 mm dust continuum emission. The HST image reveals a galaxy along the line of sight to HerBS-89a, which lies close to its northern component but is offset by ∼0. 2 southwest from the northern dust continuum peak. The lensed galaxy itself remains undetected in the HST image. The near-infrared emission of the foreground lensing galaxy is resolved in the HST image. It is dominated by a bright bulge-like structure. It also displays a weak emission that extends in the east-west direction and that could resemble a disk. A bright spot to the northwest of the spheroidal component seems to wrap around the northern peak of the dust continuum, and could be a small satellite galaxy. From the HST data, the foreground galaxy has an estimated flux density of ∼4.36 ± 0.88 µJy at 1.15 µm.
The foreground lensing galaxy is also detected in our GTC OSIRIS r-band imaging with a magnitude r = 24.5±0.3 AB mag. The deeper imaging with HiPERCAM detects the foreground galaxy in the griz bands, while in the u band the object remains undetected. The target lies very close to a bright (r = 19.5) point-like object, 3. 3 to the east, and the photometric measurement needs to take properly into account the wings of the bright neighbor.
The HiPERCAM maps were calibrated on SDSS DR12 stars present in the field. The code GALFIT (Peng et al. 2002) was used to subtract the profile of the bright point-like nearby object using a model PSF. Aperture photometry was performed on the residual image to measure the emission of the target lensing galaxy. Uncertainties were derived measuring the noise level in an annulus around the object, thus taking into account the presence of the nearby source. Finally, systematic zero-point uncertainties were propagated in the estimate of the magnitude uncertainties. Table 4 lists the derived magnitudes and the 3σ u band upper limit in the AB system (Oke & Gunn 1983).
The foreground galaxy is also detected by WISE (Wright et al. 2010) at 3.4 and 4.6 µm. However, the WISE photometry likely includes contributions from the lensing galaxy, the lensed HerBS-89a galaxy, and the nearby bright point-like source. Because of the broad WISE beam (6. 1 at 3.4 µm), it is impossible to disentangle the different components. Therefore only the GTC and HST data points have been retained in this analysis.
At r = 24.5 it is far from easy to obtain an optical spectroscopic redshift, but thanks to the sensitivity of HiPERCAM and the capability to observe simultaneously the five ugriz bands, together with the HST F110W imaging results, it is possible to derive a photometric redshift.
A122, page 11 of 33 (Polletta et al. 2007;Ilbert et al. 2009). The inset displays the distribution of χ 2 as a function of redshift; the dotted horizontal lines correspond to the 1, 2, 3σ levels, defined as χ 2 min + 1, 2.71, 6.63. The best photometric redshift for the foreground galaxy is z phot = 0.9 +0.3 -0.5 where the errors are 3σ.
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), andMAGPHYS-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. 2003Berta et al. , 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 χ 2 distribution. The photometric redshift is best constrained by the r-i color, tracing the D4000 breakfoot_3 . The lack of rest-frame near-infrared coverage and the upper limit in the u band restrict the accuracy of the photo-z measurement. Deeper data and a wider wavelength coverage would be necessary to improve its precision.
In addition, various attempts were made using templates that include emission lines (e.g., Babbedge et al. 2004), with the aim to verify if the r and F110W bands could be contaminated by lines at z ∼ 0.7, namely [Oii] (λ3727) and Hα (λ6563). However, none produced a better fit to the data.
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 12 CO(5-4) emission line of HerBS-89a was detected, a weaker narrow (∼200 km s -1 ) emission line is seen at 143.74 GHz, close to the center position of HerBS-89a (see Fig. 11). The narrowness of this emission line excludes the possibility that it could be related to HerBS-89a. If it traces instead the foreground lensing galaxy, the line could correspond to a 12 CO emission line, either detected in the J = 2-1 transition, in which case z lens = 0.604, or in the J = 3-2 transition for which z lens = 1.406. The former is consistent within the 2σ uncertainties with the photometric redshift derived from the GTC and HST data. A search for the counterpart of this emission line (at 1 mm) could allow us to verify the above identification, and eventually derive a precise spectroscopic redshift of the foreground galaxy gravitationally amplifying HerBS-89a.
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.
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) 5 , and ellipticity position angle (θ e ; measured east of north). Table 5 lists the best fit values of these parameters and their uncertainties as determined from a Markov chain Monte Carlo (MCMC) exploration of parameter space. The position of the lens is given relative to the reference position RA = 13:16:09.79 and Dec = 28:12:18.17. Because of the very high S/N of the continuum data, we use the lens model derived from the continuum data when reconstructing the spectral line data.
The left panels of Fig. 12 include the critical curve (i.e., the locus of infinite magnification in the image plane; red lines). In the case of a SIEMD with no external shear, this curve is also an iso-density contour of the deflector's mass distribution whose position and geometry can be compared to the HST image of the lensing galaxy (Fig. 9). The corresponding source-plane caustic is plotted in the right panels of Fig. 12.
The continuum and (almost) all of the spectral line zeroth moment maps are fit well by single Sérsic components in the source plane (see Table 6 and Fig. 12 for fitting results, and below for more discussion of the line data). The profile is parametrized in terms of its position, normalization I 0 , ellipticity e, position angle θ e , scale radius R s , and Sérsic index n:
Using the estimated photometric redshift of the foreground galaxy and its uncertainty (see Sect. 4.2), we estimate that the deflector mass enclosed within the Einstein radius is M lens = 9.8 +3.2 -4.9 × 10 10 M (where the upper and lower bounds are at 3σ). Statistical uncertainties in the lens and continuum source model were evaluated with two MCMC cycles. The first fixed the source model and determined the uncertainties of the lens parameters; the second fixed the lens model to its optimal parameters 5 Ellipticity is defined as e = 1q, where q is the axis ratio. and let the source parameters vary. This approach was chosen for computing speed and for ease of comparison between continuum and line results (because the latter assume the lens model derived from the continuum data, i.e., the lens model is fixed in this case). However, it should be noted that the MCMC approach underestimates the lens+source errors. Appendix A presents the conditional posterior distributions of the lens and source model parameters.
To these statistical uncertainties, we add fiducial systematic uncertainties of 10% associated with the choices of lensing potential, source brightness profile, interpolation errors, etc. The only exceptions are the source model positional errors, for which we adopt systematic uncertainties given by HPBW min HPBW maj 2
in terms of the synthesized beam dimensions HPBW min and HPBW max , the S/N of the observed detection σ I /I line , and the magnification µ (see, e.g., Downes et al. 1999). For the lens model, we estimate systematic positional errors as the astrometric errors from the continuum imaging. Because we expect systematic errors to dominate and because of computational constraints, we do not perform MCMC on the line spectral source parameters.
As mentioned above, nearly all of the spectral line data are fit well by single Sérsic components; the exception is the H 2 O(2 02 -1 11 ) zeroth-moment map, which shows faint (∼5σ, compared to the ∼10σ peak) emission to the west and northwest that motivated some exploration of a two-component model. We adopt here a single-component fit for this line, due to numerical instability in the two-component fit, but there are some systematic features in the residuals associated with the peak of the emission. The results are listed in Table 6, along with the corresponding magnification factors µ, for the high resolution data and for 12 CO(1-0). Because of the lower resolutions and S/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 12 CO(1-0) line, we infer a lensing magnification µ = 5.0, with a large systematic uncertainty mainly due to possible different choices in how the zeroth-moment map is derived (also reflecting the limited S/N of the current data set). Deeper and higher angular resolution observations of the 12 CO(1-0) emission would be helpful to better constrain the amplification factor as well as to study in greater detail the morphology (extent and kinematics) of the cold gas reservoir of HerBS-89a.
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 Notes. (a) The dust continuum emission is measured at 254.5 GHz, and the uncertainties are ±1σ. For each parameter, the systematic uncertainties are tabulated separately in parentheses on a second line (see Sect. 5). In the heading of the table, e stands for the ellipticity, θ e for the position angle (measured east of north), R s for the scale radius, n for the Sérsic index, and µ for the amplification factor. Surface brightness units of Jy beam -1 are inherited from the surface brightness units of the data. Strong covariance between R s and n (see, e.g., Fig. A.1) means that the 12 CO(9-8) and 12 CO(1-0) data are formally compatible with (larger) R s and (smaller) n like those recovered for the other reconstructions.
Sect. 6.3, and the properties of the water emission line are examined in Sect. 6.4.
The new NOEMA continuum measurements (see Table 1), combined with the data from Neri et al. (2020) and Herschel and SCUBA-2 photometry (Bakx et al. 2018(Bakx et al. , 2020)), provide an exquisite sampling of the far-infrared (FIR) and sub-mm spectral energy distribution (SED) of HerBS-89a. A total of 16 broad band measurements cover the range from 250 µm to 3 mm in the observed frame, 12 of which are spread along the Rayleigh-Jeans tail.
Figure 13 shows the continuum data available for HerBS-89a as well as the continuum measurements extracted for the nearby source HerBS-89b. The latter source shows an increase of flux at wavelengths above ∼2 mm with respect to the dust emission, which could be due to non-thermal processes and rules out the possibility that HerBS-89b is a third lensed image of HerBS-89a.
Following and extending the work by Neri et al. (2020) and Berta et al. (2016), we model the SED of HerBS-89a using three different approaches: (i) the Draine & Li (2007, hereafter DL07) dust model; (ii) a single-temperature modified black body (MBB) in the optically-thin approximation; and (iii) an MBB model in its general form.
In the case of the DL07 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), andBerta 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:
where B ν (T dust ) is the Planck function, T dust is the dust temperature, and κ ν = κ 0 (ν/ν 0 ) β is the mass absorption coefficient of dust at rest frequency ν.
The general form of the MBB differs from the optically thin approximation by the factor (1e -τ ν ), with τ ν = τ 0 (ν/ν 0 ) β , instead of the ν β term. Here ν 0 is the frequency at which τ = τ 0 = 1.0 (see Berta et al. 2016). Therefore, the general form of the MBB tends to produce best fit solutions with higher dust temperature than the optically thin case. In the optically thin approximation, three free parameters are at play: T , β and the model normalization. In the general form, ν 0 is also a free parameter.
The MBB fits are limited to the data with rest-frame wavelength λ rest > 50 µm, thus avoiding biases 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 2 kg -1 at 850 µm. The same correction is also applied to the dust masses derived with the DL07 approach.
The best fit models are shown in Fig. 13, and the results are listed in Table 7. We include all the values obtained with our fits using the continuum flux densities not corrected for the lens magnification, as well as our preferred source-plane values for the infrared luminosity and the dust mass based on the DL07 fit corrected for gravitational magnification (see Sect. 5). The CMB effect produces a steepening of the millimetric slope by ∼+0.15, comparable to the uncertainty on β. At the redshift z = 2.95 and for the dust temperature of HerBS-89a, the effect is milder than what found at z = 3.6-5.8 by Jin et al. (2019).
Finally, the VLA 3σ upper limits for HerBS-89a that are shown in Fig. 13 have been computed over an area equivalent to the NOEMA 263.3 GHz continuum (see Table 1). For comparison, a power-law spectrum S ν ∝ ν -0.8 is plotted, representing a synchrotron emission (e.g., Ibar et al. 2010;Thomson et al. 2014), normalized on the basis of the radio-farinfrared correlation (Magnelli et al. 2015;Delhaize et al. 2017).
The observed 12 CO(1-0) emission line intensity can be translated into a luminosity (in K km s -1 pc 2 ), from which an estimate of the molecular gas mass can be made, using the following relation (Solomon & Vanden Bout 2005):
where S Line ∆V = I Line is the integrated line flux in units of Jy km s -1 ; ν rest is the line rest frequency in GHz; and D L is the luminosity distance in Mpc. The line luminosity, expressed in units of L , can be derived as:
We define the molecular and total gas masses as
where M He is the mass of helium, M H 2 is the mass of molecular hydrogen and M HI the mass of atomic hydrogen. Hereafter, we assume that M mol M HI and thus that M gas ≈ M mol . The molecular hydrogen gas mass is computed from the 12 CO(1-0) luminosity applying a conversion factor α CO in units of M (K km s -1 pc 2 ) -1 :
In the Milky Way and nearby star-forming galaxies with near-solar metallicity, the empirical conversion factor α CO is α MW = 4.4 ± 0.9, already including a contribution due to helium (see, e.g., Magnelli et al. 2012;Bolatto et al. 2013;Carilli & Walter 2013;Tacconi et al. 2020). For extreme local 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 2 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 κ 850 = 0.065 m 2 kg -1 is adopted (in agreement with the results by Planck Collaboration XXI 2011, see Appendix B for details) and it is found that a value α 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 -1 pc 2 ) -1 , and we multiply it by the usual factor 1.36 to account for the helium contribution. Possible consequences of using a different value are discussed at the end of this section.
The resulting line luminosity and molecular gas mass, corrected for lensing magnification, are thus
The molecular gas mass derived from the 12 CO(1-0) emission line and the SED-based estimate of the dust mass (see Sect. 6.1) provide the basis for deriving the average gas-to-dust ratio (δ GDR ) in HerBS-89a:
where the assumption is made that the two estimates are representative for the whole galaxy, and that M mol M HI . The data lead to a value of δ GDR between 80 and 174, depending on whether the DL07 or the MBB (general form) dust mass estimate is used, with an uncertainty of ±21 based on error propagation. Interestingly, this value is consistent within the errors with those seen in typical star-forming galaxies of solar metallicity, near the so-called "Main Sequence of star formation" (MS), for which δ GDR ∼ 100 (e.g., Leroy et al. 2011;Magdis et al. 2012;Rémy-Ruyer et al. 2014). This result is nevertheless subject to the adopted value for α CO .
The total molecular gas mass (M mol ) and SFR relate to each other via the so-called depletion timescale (in years):
where the rate of star formation is computed as SFR = 1.09 × 10 -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 -1 and a depletion timescale τ depl = (3.4 ± 1.0) × 10 8 years. For its current formation rate, HerBS-89a would exhaust its molecular gas reservoir in only 340 million years, not considering any mass return to the interstellar medium or further gas inflow (see Sect. 9).
The value for M mol was derived from the 12 CO(1-0) luminosity using α CO = 3.0 × 1.36, a value that is close to the Milky Way value α MW = 4.4 M (K km s -1 pc 2 ) -1 , used for MS galaxies, rather than α SB = 1.09, often associated to starbursts or outliers in the main sequence. This choice leads to δ GDR and τ depl values that would position HerBS-89a close to or on the MS star-forming galaxies on the integrated Kennicutt-Schmidt relation of star formation (e.g., Lada 2015;Liu et al. 2019), which suggests that HerBS-89a is dominated by a "secular" star-forming mode.
However, despite the evolution of the MS normalization (increasing as a function of redshift, e.g., Elbaz et al. 2011;Whitaker et al. 2014;Schreiber et al. 2015), and because of the bending of the MS at the high stellar mass end, even at z = 3, a SFR ∼ 600 M yr -1 would hardly be associated with a galaxy on the main sequence. It is therefore more likely that HerBS-89a is located in between the MS and the starbursts region. Using α SB , the gas-to-dust ratio and the depletion time scale of HerBS-89a would both decrease by a factor of 3.74, yielding δ GDR 22-47 and τ depl = 9.1 ± 2.7 × 10 7 yr. Realistically, the actual value of α CO , on which these results depend, is likely between the two above mentioned extremes.
Further insights into the position of HerBS-89a with respect to the MS would require the knowledge of its stellar mass M * . Detecting HerBS-89a and distinguishing it from its deflector seem out of reach for currently existing optical/NIR facilities, but the James Webb Space Telescope (JWST) could provide the needed sensitivity.
Together with the measurements reported in Neri et al. (2020), we have detected in total four 12 CO emission lines in HerBS-89a, namely the J = 1-0, 3-2, 5-4 and 9-8 transitions. The 12 CO Spectral Line Energy Distribution (SLED) of HerBS-89a, normalized to the J = 1-0 transition after correcting for the lens magnification (see Table 6), is shown in Fig. 14. The velocityintegrated fluxes of the 12 CO lines increase monotonically with rotational quantum number up to J = 5-4, with a turnover around that transition, before decreasing at higher J values, with the J = 9-8 emission line about a factor of 8 times weaker than the J = 5-4 line. The peak of the CO SLED for HerBS-89a should therefore occur around the J = 5-4 and J = 6-5 transitions. This result makes HerBS-89a comparable to other starburst galaxies as shown in Fig. 14, where various SLEDs of several individual DSFGs and one quasar (the Cloverleaf) are displayed together with the center of 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.
To investigate the CO line excitation further and constrain the physical conditions of the molecular gas in HerBS-89a, we model the CO line fluxes using a Large Velocity Gradient (LVG) statistical equilibrium method (following the approach described in Yang et al. 2017, and references therein). The free parameters are the kinetic temperature of the molecular gas (T kin ), the volume density (n H 2 ), the column density of CO per unit velocity Fig. 14. 12 CO Spectral Line Energy Distribution (SLED) for HerBS-89a (red dots), compared to other well-studied local and high-z galaxies (or samples of galaxies). The CO SLEDs are normalized to the J = 1-0 transition, after correcting for the magnification, and the line fluxes of HerBS-89a have been corrected for amplification (see Sect. 5 and Table 6) before normalizing. The three samples of galaxies are from Bothwell et al. (2013) including 32 DSFGs at 1.2 < z < 4.1, Spilker et al. (2014), including 22 SPT-selected DFSGs at 2.0 < z < 5.7, and Valentino et al. (2020), including 78 "main sequence" star-forming galaxies. References to individual sources are: Cosmic Eyelash (Swinbank et al. 2011;Danielson et al. 2011); Cloverleaf (Barvainis et al. 1997;Weiß et al. 2003;Bradford et al. 2009;Riechers et al. 2011); GN20 (Carilli et al. 2010;Cortzen et al. 2020, and in prep.); M 82 (Weiß et al. 2005); and the Milky Way disk (Fixsen et al. 1999). gradient (N CO /dV), and the solid angle (Ω app ) of the source. The CO SLED only depends on T kin , n H 2 and N CO /dV. A Bayesian approach is used to fit the line fluxes generated from the model and the code emcee is adopted to perform a MCMC calculation. The resulting single component fit to the CO SLED of HerBS-89a is shown in Fig. 15. The SLED peaks around J up = 6-5, and, from the single excitation component fitting, the molecular gas density is found to be log 10 (n H 2 /cm -3 ) = 2.79 +0.72 -0.54 , the kinetic temperature T kin = 174 +157 -92 K and the CO column density log 10 (N CO/dV /(cm -2 km s -1 )) = 17.59 +0.55 -0.59 , where the errors are 3σ.
A second fit including an additional, colder and more extended component is shown in the lower panel of Fig. 15. Because of the paucity of available data points, the cold component is poorly constrained, with log 10 (n H 2 ) = 2.30 +1.06 -0.57 cm -3 , T kin = 20 +20 -7 K, and log 10 (N CO/dV ) = 15.63 +1.14 -0.79 cm -2 km s -1 , where the errors are 3σ. The warm component is in broad agreement with the previous single component fitting results.
Water is one of the most abundant molecules after H 2 and CO in the gaseous interstellar medium (ISM). It serves as a unique diagnostic for probing the physical conditions of the ISM in both local (e.g., González-Alfonso et al. 2014) and highredshift galaxies (e.g., Omont et al. 2013;Riechers et al. 2013;Yang et al. 2016, 2020, andreferences therein). In particular, since the E upper > 200 K levels of H 2 O are primarily excited through absorption of far-infrared photons from dust emission in warm dense regions, it is a useful diagnostic of the far-infrared radiation field independent of gas conditions. In addition, submillimeter water lines are the second strongest molecular emitter after the CO lines (e.g., Yang et al. 2013). The detection of the para-H 2 O(2 02 -1 11 ) emission line (ν rest = 987.9268 GHz) in HerBS-89a is an illustration thereof: the line width is broad (∆V ∼ 1100 km s -1 ) and comparable to those of the mid/high-J CO lines, consistent with the finding that the J = 2 H 2 O lines are spatially co-located with the CO lines (e.g., Omont et al. 2013;Yang et al. 2016). The observed line flux is bright, with I H 2 O = 1.59 ± 0.37 Jy km s -1 , and comparable to the line flux of the 12 CO(9-8) emission line (Table 2). The CO and H 2 O lines have morphologies and velocity fields that are broadly comparable (see Figs. 4 and7) indicating that they both trace the same warm high-density gas in HerBS-89a. Using Eqs. ( 6) and ( 7), we derive an intrinsic line luminosity (corrected for lensing magnification) of L H 2 O(2 02 -1 11 ) = (6.4 ± 1.5) × 10 7 L , yielding a ratio L H 2 O(2 02 -1 11 ) /L IR = (1.12 ± 0.36) × 10 -5 , in terms of the 8 to 1000 µm L IR (Table 7). HerBS-89a lies on the correlation between 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 2 O emission that originates in very dense, heavily obscured starformation dominated regions (see, e.g., González-Alfonso et al. 2014;Yang et al. 2016). The steeper-than-linear growth is likely the result of an increase of the optical depth in the dust continuum (at ∼100 µm) with increasing infrared luminosity. Finally, the ratio of L H 2 O(2 02 -1 11 ) /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).
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 12 CO(9-8) line, which has higher S/N than the H 2 O(2 02 -1 11 ) line and is therefore more promising as a kinematic probe. In order to improve the S/N, we rebin the 40 km s -1 data to a resolution of 120 km s -1 , yielding a total of seven independent velocity channels that are suitable for 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 A.3 indicates that these reconstructions are formally acceptable, although the five channels blue-ward of 168 km s -1 have scale radii that are notably smaller than those of the two most highly redshifted channels. The centroids of all seven channels fall within the envelope defined by the reconstructed dust continuum emission, and in turn by the reconstruction of the 12 CO(9-8) zeroth moment map (second row of Fig. 12), giving confidence that the overall spatial extent of the warmest and densest material in HerBS-89a is understood.
The center panel of Fig. 17 presents a source-plane zeroth moment map constructed by summing the source-plane 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érsicprofile sources in most of the velocity channels, the moment map derived from them suggests that we have reliable information about the internal source structure on much finer scales than is actually the case. It is also worth noting that when the moment map derived from the reconstructed channel maps is lensed forward into the image plane, it provides a poorer match to the observed 12 CO(9-8) moment map than does the single Sérsic-component reconstruction in the second row of Fig. 12. This apparent paradox occurs because small noise fluctuations in the observed channel maps have unavoidably affected the source centroids, shapes, and orientations of the corresponding source-plane reconstructions. Similar inconsistencies have been seen by other authors who have separately delensed both channel maps and integrated moment maps (e.g., Dong et al. 2019). Use of prior information on the expected relationships between emitting regions in adjacent velocity channels (e.g., Rizzo et al. 2018;Young et al., in prep.) would be one way to mitigate this tendency, but we have not used such an approach here.
Figure 17 provides a further illustration of the degree to which small noise fluctuations in the observed channel maps can perturb the centroids of the source-plane reconstructions. Initially inspired by the apparent bent morphology (and nominal associated velocity gradients; see below) in the central panel of Fig. 17, we constructed new, roughly half-line moment maps from the bluest four and reddest three of our seven 120 km s -1 channels, and reconstructed these using the lens model optimized in Sect. 5. The sum of these two reconstructions (i.e., another zeroth moment map) is presented in the right panel of Fig. 17. The apparent double morphology might at first right suggest a pre-coalescence merging system; however, the 'doubleness' of the morphology stems at least in part from the general compactness of our Sérsic profile reconstructions. Moreover, notable aspects of the internal source structure (e.g., the distinct southeastern offset of the reconstructed -312 to -192 km s -1 channel map) in the center panel do not translate to the reconstructions used to derive the right panel. This level of inconsistency raises the question of whether detailed information on the intrinsic or the internal source structure (as distinct from the overall extent of emission) can be robustly recovered from this dataset, and suggests that the apparent small-scale structure in the resulting zeroth-moment maps has to be interpreted with caution.
With the above discussion as preamble, we turn now to Fig. 18, which shows source-plane velocity fields derived from the seven-channel (120 km s -1 resolution) reconstructions of Fig. A.3 in two different ways. The left panel presents an intensity-weighted first moment map calculated in the usual way. We find that the dominant velocity gradient is in an east-west direction, roughly consistent with the east-west elongation of the reconstructed dust continuum emission. The comparison of the center and right panels in Fig. 17 discussed above suggests that the perturbation from an "exactly" aligned east-west velocity gradient due to the most highly blue-shifted emission may not be robust.
Assuming a single rotationally supported disk, the radial profile of the dynamical mass of the system is given by
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 (1e), since the velocity gradient does not fully correspond to the continuum Sérsic profile in the source plane, we opt to keep the sin(i) term explicit in our computations. The constant f is a dimensionless scale factor that depends on the structure of the galaxy and has value f 1.0 for a thin disk embedded in a massive spheroid (Lequeux 1983).
Along the major axis of the continuum profile (red line in the top-left panel of Fig. 18), we extract the velocity curve V(R) shown in the middle-left panel. The corresponding bottom-left panel presents the radial dynamical mass profile, which leads to an estimate of M dyn sin 2 i ≈ 5.1 × 10 10 M within R = 1.25 kpc.
In contrast, the top-right panel of Fig. 18 presents a first moment map in which masking has been used to highlight the Fig. 18. Kinematic and dynamical analysis of the 12 CO(9-8) emission line in the source plane of HerBS-89a. Top panels: left: first moment map, compared to the Sérsic profile of the dust continuum emission. The contours start from 1/20 of the peak intensity with steps of 1/20th. The red line marks the continuum major axis. The reference position is defined where the rest-velocity (V = 0 km s -1 ) and is denoted with a star. From this position we calculate the radial distance. The cross at the lower right corner marks the phase center of the NOEMA observations. Right: based on Figs. A.3 and 17, two possible components are identified; the first moment maps of the two components are here shown, after shifting the velocity scale by ±289 km s -1 for clarity. The axes of the two components are shown as green and orange lines. The dashed line marks the contour of the velocity field shown in the top-left panel. Middle panels: velocity curves extracted along the axes marked in the top panels. Colors match those of the respective axes, with magenta points extracted from the top left moment map and the green and orange points extracted from the top right moment map. In the two right-most panels, the data are shown both before (small symbols) and after (large symbols) shifting by ±289 km s -1 . The amplitude of the error bars is shown in the left-most panel. Bottom panels: dynamical mass profile as a function of radius, as derived from the velocity curves. two regions of strongest integrated emission (compare also to Figs. 17 andA.3). These two regions velocity scales have been artificially offset by ±289 km s -1 , corresponding to the separation between the two peaks in the integrated CO spectrum (see Fig. 5) in order to better present them as "candidate" rotating disks. Velocity curves extracted along the axes of these two possible structures (green and orange lines) are shown in the middle row of Fig. 18, both before (small symbols) and after (large symbols) applying the ±289 km s -1 shift.
The individual dynamical masses of these components, assuming rotational support as above, would be M dyn sin 2 i ≈ 2.6 × 10 10 M within R = 0.85 kpc for the eastern component and M dyn sin 2 i ≈ 2.2 × 10 10 M within R = 0.40 kpc for the western component, again under the assumption that f = 1.0 in Eq. ( 15). The bottom panels of Fig. 18 show the dynamical mass radial distributions, limited to the extension of the first-moment maps.
Assuming a baryon fraction f baryon (e.g., Tacconi et al. 2020), and supposing that at high redshift the gas content of a galaxy A122, page 21 of 33 is dominated by its molecular component and the atomic gas mass M HI can be neglected (see, e.g., Tacconi et al. 2018), the dynamical mass could yield an estimate of the stellar mass M * = f baryon M dyn -M mol -M HI and of the molecular gas fraction f mol = M mol /( f baryon M dyn -M mol ), thereby solving the question of whether or not the galaxy lies on the MS. However, the dynamical mass profiles derived for HerBS-89a do not extend to large enough radii to allow for a determination of the total mass.
For the sake of completeness, we also performed the reconstructions of multiple velocity channels across the H 2 O(2 02 -1 11 ) emission line, adopting the same approach as above, using a coarser velocity re-binning to boost the S/Ns (see Appendix A). As already suggested by the slightly different observed and reconstructed source-plane zeroth moment maps (see Fig. 12), the individual 12 CO(9-8) and H 2 O(2 02 -1 11 ) velocity channels reconstruct to slightly different locations in the source plane (see Figs. A.3 and A.4). However, a first moment map derived from the reconstructed H 2 O(2 02 -1 11 ) channel maps shows a velocity gradient similar to the one shown in the top left panel of Fig . A.3 for the 12 CO(9-8) emission line, albeit at a slightly different position angle, suggesting once more that this aspect of the system's kinematics has been robustly recovered. Due to the higher S/N of the 12 CO(9-8) emission line, we have limited the kinematic 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 12 CO emission lines in this system. Observations with higher S/N and angular resolution that can support a more detailed lens model will be needed to resolve the one vs. two sources question in a definitive way for this system, shed further light on the properties of this galaxy, unveil the finer details of its structure and kinematics, and fully recover the secrets of its starburst nature.
In the following subsections, we present and discuss the various molecular absorption and emission lines other than CO and H 2 O that have been detected in the spectrum of HerBS-89a. The individual molecular lines are described starting with the very dense molecular gas tracers, namely: HCN and/or NH, and NH 2 (Sect. 8.1), followed by the two molecular ions, OH + and CH + , which probe the low density molecular gas with low H 2 fractions (Sect. 8.2).
8.1. Very dense gas tracers: NH 2 , and HCN and/or NH NH 2 . Amidogen (NH 2 ) 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 3 ). It is an asymmetrical molecule with a 2 B 1 ground electronic state characterized by a complex rotational spectrum (similar to H 2 O + ). It was first detected in the interstellar medium (in SgrB2) by van Dishoeck et al. (1993) and Goicoechea et al. (2004). NH 2 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 (2 02 -1 11 )(5/2-3/2) emis-sion line is seen next to the OH + (1 0 -0 1 ) absorption line, and, possibly, the NH 2 (2 02 -1 11 )(3/2-1/2) emission line as well.
In the case of HerBS-89a, two emission lines of NH 2 are detected. The first is the NH 2 (2 02 -1 11 )(5/2-3/2) emission line, which overlaps with the OH + (1 0 -0 1 ) absorption line, such that only the red-shifted part of the NH 2 emission line is apparent in the spectrum. However, the nearby NH 2 (2 02 -1 11 )(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 20 -2 11 )(5/2-3/2) is weaker but uncontaminated by nearby lines.
As the data are of low S/N, a detailed study of NH 2 , first reported here in a high-z galaxy, remains out of reach. However, from the velocity-integrated emission map of the highangular resolution data of the NH 2 (2 20 -2 11 )(5/2-3/2) emission line (Fig. 4), it is clear that the NH 2 emission closely resembles the distributions of the 12 CO(9-8) and water emission lines, with the southern peak being stronger, in accordance with the fact that this molecule also traces the high-density gas. Dedicated future observations can further explore the properties of this molecule in HerBS-89a or in other high-z galaxies to study the properties of their dense gas reservoirs.
HCN(11-10) and/or NH. The absorption feature observed at 246.72 GHz in HerBS-89a may be due to redshifted and/or NH(1 2 -0 1 ). 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 H 2 3 × 10 4 cm -3 ), with critical densities 100 to 1000 times higher than of CO. Although HCN emission lines are typically one to two orders of magnitude fainter than CO emission lines, HCN has been observed in its low-J transitions in local starburst galaxies, revealing a tight correlation between the HCN luminosity and the star formation rate (e.g., Gao & Solomon 2004;Zhang et al. 2014). Emission lines of low-J HCN have also been detected in a few high-z galaxies (Weiß et al. 2007;Riechers et al. 2010;Oteo et al. 2017;Cañameras et al. 2021, and references therein), and HCN is seen in the stacked restframe 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 absorptionfoot_4 . These transitions of HCN have critical densities around 10 9 -10 10 cm -3 and probe very dense (n H 2 ≈ 3 × 10 6 cm -3 ) gas. The high-J lines of HCN are thought to be populated by radiative pumping of infrared photons (rather than by collisional excitation) in an intense high-temperature (>350 K) radiation field Rangwala et al. (2011).
In HerBS-89a, the optical depth of the 246.72 GHz line is estimated to be τ(ν 0 ) ∼ 0.24 ± 0.12 yielding, in the case of HCN, a column density of (0.38 ± 0.1) × 10 14 cm -2 (Table 3). Future observations of both high-J HCN lines (in absorption) and low-J lines (in emission) in HerBS-89a would allow us to confirm if this absorption line is indeed due (or partly due) to HCN, and, once verified, allow us to probe and constrain the properties of the reservoir of extremely dense gas in this high-z starburst galaxy.
-NH (imidogen) is a molecular radical that is linked (as for NH 2 ) to the formation and destruction of NH 3 . It was first detected in the diffuse interstellar medium (Meyer & Roth 1991) and later in the dense molecular gas of SgrB2 (Cernicharo et al. 2000;Goicoechea et al. 2004). The Herschel SPIRE-FTS spectra of Mrk 231 and Arp 220 (González-Alfonso et al. 2018;Rangwala et al. 2011) show the presence of NH with various transitions detected. NH has a spectrum similar to that of OH + , since these molecules are isoelectronic. Many of the transitions of NH overlap with OH + transitions, in particular NH(1 2 -1 0 ) with OH + (1 0 -0 1 ), with similar column densities derived for OH + and NH in Mrk 231 (González-Alfonso et al. 2018). In the case of HerBS-89a, attributing all of the 246.72 GHz absorption line to NH implies a column density of (6.7 ± 1.9) × 10 14 cm -2 , which is also comparable to the column density of OH + (Table 3). Together with NH 2 , the likely detection of NH in HerBS-89a opens up the possibility of exploring the chemistry of N-bearing molecules in a high-z galaxy by measuring other transitions of NH and searching for NH 3 , and exploring whether (as is the case in SgrB2) the derived column densities are compatible with grain-surface chemistry and sputtering by shocks (see, e.g., Goicoechea et al. 2004).
8.2. Low density gas tracers: OH + and CH + OH + . The molecular ion OH + is a reliable and powerful tracer of molecular outflows in DSFGs as has been demonstrated with the Herschel SPIRE-FTS studies in Arp 220 (Rangwala et al. 2011) andMrk 231 (González-Alfonso et al. 2018), and in highz galaxies with NOEMA in HLFS 3 (Riechers et al. 2013) and ALMA (Indriolo et al. 2018, who reported OH + (1 1 -0 1 ) absorption line in two lensed sources, i.e., the Cosmic Eyelash and SDP.17). As discussed in Indriolo et al. (2018), OH + absorption is thought to arise in the extended cool, diffuse, low H 2 fraction gas that surrounds galaxies. OH + has been found to be useful in constraining the cosmic-ray ionization rate of atomic hydrogen, ζ H 2 , particularly in gas with a low H 2 fraction. The three ground state lines of OH + (all are split into hyperfine components, which are blended due to the broad line widths measured in galaxies) are generally observed in absorption. In a few cases, they can be seen in emission as a red wing, as observed in Mrk 231 or Arp 220, or as recently reported for the OH + (1 1 -0 1 ) line in the z = 6.03 quasar SDSS231038+1855 by Li et al. (2020) and the hot dust-obscured galaxy W0410-0913 at z = 3.63 by Stanley et al. (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 1 lines (J = 0, 1, 2) are seen in absorption (Fig. 5 middle panel). The (1 1 -0 1 ) and (1 2 -0 1 ) transitions are the strongest, while the (1 0 -0 1 ) transition is nearly five times weaker (Table 2). None of the OH + lines in HerBS-89a show any indication of emission above the continuum as in Mrk 231 or Arp 220 (González-Alfonso et al. 2018;Rangwala et al. 2011). All three lines have comparable line widths with ∆V ≈ 500 km s -1 and the inferred OH + column density is N OH + ∼ 10 15 cm -2 (Table 3), comparable to the value derived for the Cosmic Eyelash (Indriolo et al. 2018). The OH + spectra show clear indications of red-shifted gas, with all three transitions peaking in opacity at ≈+100 km s -1 relative to the systemic velocity, including the (1 0 -0 1 ) transition after the contamination of the NH 2 emission line is accounted for (see 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.17foot_5 . The CH + (1-0) line is found both in emission and absorption, providing critical information on the mechanisms of mechanical energy dissipation and/or strong UV irradiation in these high-z galaxies. Due to its high critical density for excitation (∼10 7 cm -3 for the J = 1-0 transition), the CH + line is seen in emission in dense (>10 5 cm -3 ) gas, probing shocked regions powered by galactic winds, and in absorption, in large (>10 kpc) reservoirs of turbulent, cold (∼100 K), lowdensity (∼100 cm -3 ) gas (see, e.g., Falgarone et al. 2017).
In HerBS-89a, the CH + (1-0) is detected in absorption with a width of ≈400 km s -1 , comparable to that of the OH + lines, and is also red-shifted to ≈+100 km s -1 (Fig. 5). The inferred CH + column density is estimated to be N CH + ∼ 10 14 cm -2 (Sect. 8), similar to what has been derived in the case of the Cosmic Eyelash (Falgarone et al. 2017).
The low S/N, tentative detection of CH + emission could indicate that in HerBS-89a there is also dense gas (>10 5 cm -3 ) with a significant velocity dispersion (the line has a full width at zero intensity of 1000 km s -1 and is centered at the systemic velocity), as was reported for three of the sources studied by Falgarone et al. (2017). However, higher quality data are required to probe this dense gas component in greater detail.
In HerBS-89a, the three OH + ground-state 1 J -0 1 lines (J = 0, 1, 2) and the CH + (1-0) line are all seen in absorption and red-shifted with respect to the systemic velocity (defined by the CO emission lines) by ≈+100 km s -1 . As shown in the velocity maps of these absorption lines, the red-shifted gas is distributed over the southern and northern dust continuum peaks (in the image plane), uniformly covering the distribution of the blueand red-shifted gas traced in the 12 CO(9-8) and water emission lines (Fig. 4). The kinematics of the absorbing gas in HerBS-89a indicate that this red-shifted gas is therefore not kinematically related to the dense molecular gas. We argue here that we are in fact tracing, in both OH + and CH + , low density gas that is flowing into the central regions of HerBS-89a.
The number of high-z sources for which measurements of absorption by molecular ions such as CH + or OH + are available is still low; however, it is noteworthy that the majority of the sources for which observations of these species have been made do show blue-shifted absorption lines, indicating vigorous outflow activity. Out of the five high-z lensed galaxies where CH + (1-0) has been found in absorption, only one sourcefoot_6 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 -1 with respect to the systemic velocity as defined by the CO emission lines. A similar velocity shift in the Cosmic Eyelash is also observed in the OH + (1 1 -0 1 ) and H 2 O + (1 11 -0 00 ) absorption lines (Indriolo et al. 2018). Recently, Butler et al. (in prep.) report OH + (1 1 -0 1 ) absorption lines 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 In contrast to the above trends, the detection of red-shifted molecular absorption lines in high-z systems seems much less common. Only two cases are known, namely the Cosmic Eyelash (Falgarone et al. 2017;Indriolo et al. 2018) and HerBS-89a (here described), out of the total of 18 high-z sources observed to date for which CH + (1-0) and/or OH + (1 1 -0 1 ) (or other tracers such as OH, H 2 O or H 2 O + ) have been seen in absorption.
Direct observational evidence of inflow is indeed scarce. Signatures of infalling gas are notoriously difficult to observe: the accreting material may have low metallicity, low density, and/or a small covering factor (less than 10% at z ∼ 1.5; see, e.g., Fumagalli et al. 2011), and it might be ionized or obscured by outflows. The need for a good alignment of an accretion structure with a background continuum source along the line of sight lowers the chance of detection even further.
The findings here described for molecular absorption lines in high-z galaxies are in line with the scarcity of infall evidence gathered from previous studies. Using optical spectroscopy, Rubin et al. (2012) report evidence of cool, metal-enriched gas accretion onto galaxies at z ∼ 0.5. Wiseman et al. (2017) found a rare case of inflow of metal-poor gas from the intergalactic medium onto a z = 3.35 galaxy. More recently, Ao et al. (2020) report evidence of infalling gas in a Lyman-α blob at z = 2.3, and Daddi et al. (2020) presented observations revealing cold gas filaments accreting toward the center of a massive group galaxy at z = 2.91. In the local Universe, a few examples of inflow activity have been found through Herschel spectroscopic observations of compact luminous infrared galaxies. These include NGC 4418 and Zw 049.057, where reversed P-Cygni profiles, a clear signature of infall, were revealed in the [OI] 63 µm fine-structure lines (González-Alfonso et al. 2012;Falstad et al. 2015); Arp 299 where the absorption of the ground state of OH at 119 µm is found red-shifted by ∼+175 km s -1 , suggesting a low excitation inflow in this source (Falstad et al. 2017); Circinus, which shows an unambiguous OH 119 µm inverted P-Cygni profile (Stone et al. 2016); and IRAS 10091+4704, the starburstdominated 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 -1 (Herrera-Camus et al. 2020).
Theoretically, gas flowing into less massive (M ≤ 10 12 M ) galaxies should be dynamically and thermally cold, while more massive halos receive most of their baryons as hotter (T > 10 5 K) gas (e.g., Dekel & Woo 2003;Kereš et al. 2005). Thus cold, dense, metal-poor circum-galactic gas is often interpreted as direct evidence of accretion. Metal-poor Lyman-limit systems have been used as tracers of accretion in observations (e.g., Lehner et al. 2013) and simulations (e.g., Fumagalli et al. 2011).
Inflow and outflow rates of galaxies are constrained by cosmological simulations (e.g., Kereš et al. 2005) and analytic models based on the mass assembly history of dark matter halos (Bouché et al. 2010;Lilly et al. 2013). Erb (2008) reproduces the gas mass fraction and metallicity distribution of galaxies with a simple chemical model involving gas inflow and outflow. Inverting their reasoning, Yabe et al. (2015) constrain gas inflow and outflow rates based on the observed stellar masses, gas mass fractions, and metallicities of star-forming galaxies at z ∼ 1.4.
Taken together, the above evidence strongly indicates that galaxies possess large reservoirs of circum-galactic gas eligible for accretion (Tumlinson et al. 2017). Nevertheless, evidence for fuel does not necessarily imply the presence of fueling.
We can describe the gas flow as a cylindrical channel that brings the gas onto the galaxy (see Fig. 19 for a schematic representation). Given the distance d inflow of the infalling clouds from the galaxy (i.e., given the line-of-sight distance that the gas has yet to travel), and assuming that the accreting gas travels at a constant speed V inflow , the mass inflow rate can be computed as follows:
Here X denotes either the OH + or CH + molecular ion; N l is the column density of the inflowing gas; m X is the mass of the molecular ion; A obs is the de-magnified projected area covered by the inflowing gas in the sky; φ is the angle between the line of sight and the inflow direction; and V obs is the measured red-shifted velocity of the line. With these definitions, the term V obs / cos φ represents the de-projected inflow velocity V inflow , and the cross-sectional area of the cylinder perpendicular to the flow is given by A inflow = A obs / sin φ.
The N l (OH + ) value measured from the high angular resolution spectra for OH + (1 1 -0 1 ) and (1 2 -0 1 ) is (9.6-11.2 × 10 14 ) cm -2 (Table 3). The mass of the OH + molecular ion is m OH + = 2.824 × 10 -23 g. The observed projected area can be computed from the Sérsic radius R s of the OH + source plane reconstruction (0.21-0.23 kpc; Table 6) yielding A obs = 0.13-0.17 kpc 2 . The observed velocity of the flow is 100 km s -1 . The distance that the inflowing gas travels from the observed position to the center of the galaxy is one of the most uncertain quantities in Eq. ( 16). Typical distances found in the literature are 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
and
modulo the angle φ between the flow direction and the line of sight, which is unknown.
In order to transform these results into the total inflow rate, we assume that the gas mass is dominated by hydrogen, and we adopt as abundances relative to hydrogen the quantities f OH + = 5.7 × 10 -8 (computed as the Galactic average value from Table 5 of Indriolo et al. 2015) and f CH + = 7.6 × 10 -9 (Godard et al. 2014;Indriolo et al. 2018). By dividing Eq. ( 16) by f X and multiplying by the ratio m H /m X of the molecular ion's mass relative to H, we derive mass inflow rates of Ṁinflow sin φ cos φ ∼ 3-14 M yr -1 from OH + , and ∼6-24 M yr -1 from CH + .
In the so-called "bathtub" or "gas regulator" model (Bouché et al. 2010;Lilly et al. 2013;Sánchez Almeida et al. 2014;Peng & Maiolino 2014;Dekel & Mandelker 2014;Somerville & Davé 2015, among many others), the evolution of the gas content of a galaxy can be formulated as follows:
where Φ is the inflow rate; SFR = M gas , as per the Schmidt-Kennicutt star formation law (Schmidt 1959;Kennicutt 1998b); R is the fraction of mass of newly formed stars quickly returned to the ISM through stellar winds and supernovae; Ψ = λSFR is the outflow rate, linked to the SFR by the mass loading factor λ; finally, η rec is the fraction of gas ejected by outflows falling back onto the galaxy. Additional factors, 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 2 , pockets of dense gas should be present; second, the most important formation channel for OH + (and CH + ) requires high temperatures, which could be produced in an accretion shock in the infalling gas as it enters the halo of the galaxy, although this scenario remains speculative. Yabe et al. (2015) use a simple analytic model of galaxy chemical evolution to constrain the inflow (and outflow) rates of three samples of star forming galaxies at z = 0, 1.4 and 2.2 (Erb et al. 2006;Peeples & Shankar 2011, and their work). They show that the average gas inflow rate of the samples is 2 to 3× their SFRs. Turning into physical units, Ṁinflow is a few solar masses per year at z = 0 growing to ∼100 M yr -1 at z = 2.2. Model predictions reproduce these values for dark matter halos masses of M DH ∼ 2 × 10 12 M (e.g., Bouché et al. 2010;Yabe et al. 2015). If the SFR is fed by the inflow, then it is naturally proportional to the inflow rate. The detailed process of gas accretion and its relation to the fuel of SFR is more complex and still poorly understood. The gas fraction and SFR efficiency (i.e., the reciprocal of τ depl ) come into play. How they relate to
Ṁinflow is yet to be fully investigated by observations, because of the rare occurrence of direct inflow measurements.
In this frame, HerBS-89a would require a very high inflow rate to sustain its large measured star formation rate. Although the angle between the line of sight and the inflow is unknown, the values derived above seem too low to support this scenario.
Alternatively, the flow could be part of a galactic fountain, in which a fraction of the gas expelled by the starburst winds falls back onto the galaxy (the term η M gas in Eq. ( 19)). In this case, an outflow should also be detected, as described by -λ M gas . No evidence of any outflow feature was found in the spectra of HerBS-89a, neither in the form of blue-shifted emission (e.g., inverted P-Cygni profile) nor through absorption components of the detected lines. Nevertheless, it is worth noting that the absorption lines of OH + here reported are not isolated, and have emission or absorption lines of other molecular species on their blue sides, perhaps preventing the identification of a blueshifted 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.
We report new millimeter NOEMA and VLA observations of the lensed starburst galaxy HerBS-89a at z = 2.95 that are complemented with near-infrared and optical data (obtained with the HST and GTC) that reveal the foreground lensing galaxy. The ability to process large instantaneous bandwidths with NOEMA enables an in-depth exploration of a sizeable part of the 1 mm spectrum of HerBS-89a, covering in total 121 GHz in the restframe, leading to the detection of transitions of molecular ions and molecular species that have not been previously observed in a high-z galaxy. In addition to the strong emission lines of 12 CO(9-8) and H 2 O(2 02 -1 11 ), we detect the three OH + groundstate 1 J -0 1 lines (J = 0, 1, 2) and the CH + (1-0) line, all four seen in absorption and red-shifted with respect to the systemic velocity of HerBS-89a, two transitions of NH 2 in emission, and HCN(11-10) and/or NH(1 2 -0 1 ) lines, seen in absorption. modeling the lensing of HerBS-89a, we reconstructed the dust continuum emission and the molecular emission lines in the source plane deriving the morphology and the kinematics of this high-z galaxy.
The main conclusions of this paper are as follows: -The high-angular resolution (0. 3) images of HerBS-89a in dust continuum and the emission lines of 12 CO(9-8) and H 2 O(2 02 -1 11 ) show two distinct components -one to the north and a brighter slightly more extended component to the A122, page 25 of 33 south, which are linked by an arc-like 1. 0 diameter structure reminiscent of an Einstein ring. HerBS-89a is gravitationally amplified by a foreground massive galaxy at z phot = 0.9 +0.3 -0.5 , with an estimated deflector mass of M lens = 9.8 +3.2 -4.9 ×10 10 M . Using lens 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 12 CO(1-0) measurements, we derive a mass for the molecular gas reservoir of M mol = (2.1 ± 0.4) × 10 11 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 H 2 ∼ 500 cm -3 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 12 L , and a dust mass of M dust = (2.6 ± 0.2) × 10 9 M , yielding a dustto-gas ratio δ GDR ≈ 80. The star formation rate is estimated to be SFR = 614 ± 59 M yr -1 and the depletion timescale depl = (3.4 ± 1.0) × 10 8 years. -Lens modeling suggests that the warmest and densest material in HerBS-89a (i.e., as traced by dust, 12 CO(9-8), and H 2 O(2 02 -1 11 ) emission) is distributed across a region a few kpc across. We cannot distinguish definitively between a single-disk and a two-component description of this distribution on the basis of existing observations, but dynamical and molecular gas mass estimates are consistent in either scenario, and the two-component scenario in particular provides a natural explanation for the very large observed CO velocity widths. -A series of molecular emission and absorption lines are detected in HerBS-89a that trace very dense gas (up to n H 2 ∼ 10 9 -10 10 cm -3 ). NH 2 is seen in emission and detected in two transitions, while an absorption feature at 246.72 GHz corresponds to HCN(11-10), NH(1 2 -1 0 ) or a combination of both. Future observations of additional transitions of these species (particularly HCN and NH) would be needed to further explore the densest regions of HerBS-89a. -The three fundamental transitions of OH + and CH + (1-0) are all seen in absorption and all red-shifted by ≈+100 km s -1 with respect to the systemic velocity of HerBS-89a. The molecular ions OH + and CH + both trace low-density (n H 2 ∼ 100 cm -3 ) gas. The fact that the velocity red-shifted absorption lines cover substantially the distribution of the CO and water emission lines in the image (and source) planes indicates that this low-density gas must be kinematically unrelated to the dense gas. We here argue that we are tracing a rare case of inflow of low-density gas, indicating that HerBS-89a is accreting matter from its surroundings. Based on the available data, we have estimated that the mass inflow rate is on the order of ≈10-20 M yr -1 . Various scenarios are possible to account for the observations including inflowing gas from the intergalactic medium, a galactic fountain or a stream of matter between the two merging components. The results here described highlight the importance of processing large instantaneous bandwidth in the study of high-z galaxies, as such observations allow us to efficiently measure multiple lines and explore in an unbiased way the entire available frequency range for the presence of molecular species. In the case of HerBS-89a, we were able to probe this system from the very highest to the lowest molecular gas densities with 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.
The adopted frequencies correspond to the strongest of the hyperfine transitions for the specific ∆J(Indriolo et al.
2015, and references therein).2 The rest frequency of the CH + (1-0) transition is taken from Müller (2010).
http://deneb.astro.warwick.ac.uk/phsaap/hipercam/ docs/html/
The 4000 Å discontinuity is defined as the ratio between the average flux density between 4050 and 4250 Å and that between 3750 and 3950 Å(Bruzual 1983).
The HCN(11-10) transition was not discussed inRangwala et al. (2011).
The first detection of CH + was made in the far-infrared spectrum of the planetary nebula NGC 7027 through rotational lines from J = 2-1 to 6-5 seen in emissionCernicharo et al. (1997).
In the case of the galaxy G09v1.40 reported inFalgarone et al. (2017), adopting the precise redshift of z = 2.0924 ± 0.0001(Yang et al. 2017) shifts the CH + (1-0) absorption line to ∼-200 km s -1 with respect to the systemic velocity of the source as defined by the CO emission lines; the CH + (1-0) absorption line is therefore blue-shifted in this system -see alsoButler et al. (in prep.).