© The Authors 2023

1. Introduction

The interstellar medium (ISM), comprising interstellar gas, dust, and cosmic rays (CRs), is a fundamental part of the ecosystem of galaxies. The cold ISM (gas and dust with T <  10³ K), which we simply refer to as the ISM hereafter unless otherwise specified, exhibits complex physical structures across different scales and plays a crucial role in a range of important physical processes, including star formation (Spitzer 1978). Energetic phenomena in the ISM, including shocks and jets, along with interactions with the electromagnetic radiation field and CRs emitted from various sources such as stars, supernovae, and active galactic nuclei (AGN), foster conditions that enable a vast number of chemical reactions (Tielens 2013, 2021). Consequently, a comprehensive understanding of the ISM in galaxies requires observational studies that delve into its chemical complexity (e.g., Omont 2007; Viti 2016).

Astrochemical studies of nearby galaxies have flourished over the past two decades, driven by advancements in radio/(sub)millimeter instruments. Deep spectral line surveys covering large frequency ranges (see Table 1 of Martín et al. 2021 for a list of extragalactic line surveys and references) allow us to study a wide variety of ISM tracers without preselecting the canonical species, such as ¹²C¹⁶O (CO hereafter), thereby offering an unbiased perspective of the ISM by appreciating its chemical complexity. These line surveys of nearby galaxies have unveiled a complicated picture of how different ISM phases interplay with radiation in the process of star formation and AGN activities and have greatly enriched our understanding of the evolution and ecosystem of galaxies. Interferometers such as the NOrthen Extended Millimeter Array (NOEMA) and the Atacama Large Millimeter/submillimeter Array (ALMA) now enable different gas tracers in nearby galaxies to be mapped down to the scales of giant molecular clouds, unveiling significant spatial variations in the ISM chemistry that are distinct from those found in the Milky Way, potentially arising from super star clusters and AGN (e.g., Meier et al. 2015; Martín et al. 2021; Sakamoto et al. 2021).

However, in contrast to local galaxies, high-redshift galaxies have remained largely unexplored (though there have been several detections in molecular absorbers against bright quasars; Muller et al. 2011, 2014; Tercero et al. 2020), primarily due to sensitivity limitations. The lack of high-sensitivity line surveys at high redshifts prevents us from attaining a comprehensive understanding of the ISM in the early Universe, which often hosts the most extreme ISM conditions and is a key laboratory for testing physical and chemical theories of the ISM.

A breakthrough came with the discovery of large samples of strongly gravitationally lensed dusty star-forming galaxies at high redshifts (e.g., Negrello et al. 2010; Vieira et al. 2013; Cañameras et al. 2015). Lensing magnification enables us to achieve sensitivities that allow for the detection of more complex molecules beyond CO. These strongly lensed infrared-luminous galaxies are among the most active dusty starbursts in the early Universe (e.g., Blain et al. 2002; Casey et al. 2014; Hodge & da Cunha 2020). They are rich in molecular gas and dust and feature bright emission and/or absorption lines, which are often highly excited by intense radiation fields, such as infrared (IR), ultraviolet (UV), and X-ray. Additionally, CRs, stellar winds, and shocks can also serve as excitation sources. The discovery of these galaxies has paved the way for in-depth spectral line surveys, particularly in galaxies with extreme ISM conditions, facilitating an unparalleled insight into their ISM. The spectral richness of such high-redshift dusty galaxies has been demonstrated in recent dedicated observations of the lensed Cloverleaf quasar at z = 2.56 where a variety of molecules have been reported, including HCN, HNC, HCO⁺, CN, and CO isotopologs (see Guélin et al. 2022, and references therein). Stacking together the spectra of high-redshift dusty star-forming galaxies from ∼200 to 800 GHz in the rest frame (Spilker et al. 2014; Reuter et al. 2023; Hagimoto et al. 2023), several molecular lines in addition to CO have also been detected, including H₂O, H₂O⁺, ¹³COCO, HCN, HCO⁺, CH and CN. However, stacked spectra only reflect the averaged quantities of the sample, and results can vary with different weighted average methods, making it difficult to interpret the physical properties.

Therefore, to accurately depict the ISM properties, it is crucial to conduct in-depth spectral line surveys in individual sources. However, individual detections of molecular gas tracers beyond CO lines remain scarce. Apart from the Cloverleaf (Guélin et al. 2022), only a limited number of high-redshift sources have been reported with detections of a few species, predominantly only single transitions per molecule (e.g., Solomon et al. 2003; Vanden Bout et al. 2004; Carilli et al. 2005; Gao et al. 2007; Guélin et al. 2007; Riechers et al. 2006, 2007, 2010; Riechers et al. 2011a; Danielson et al. 2011; Oteo et al. 2017; Béthermin et al. 2018; Cañameras et al. 2021; Rybak et al. 2022, 2023). As a result, our understanding of the high-redshift ISM remains limited.

In this paper, we report the first results of NOEMA line surveys toward two strongly lensed high-redshift galaxies of contrasting extremes, APM 08279+5255 and NCv1.143, from the SUNRISE (Submillimeter molecUlar liNe suRveys in dIstant duSty galaxiEs) project. This survey project includes multiple deep broadband millimeter spectral line surveys toward four high-redshift targets – APM 08279+5255, NCv1.143, G09v1.97, and BR 1202+0725 – using NOEMA (this paper) and ALMA (Yang et al., in prep.). The first results presented in this work highlight the detection of a rich set of emission lines in the rest-frame frequency from about 330 to 550 GHz in two sources with very different ISM environments – the dusty quasar APM 08279+5255 and the starburst galaxy NCv1.143 (throughout this work, we use the term “ISM environments” to refer either AGN- or starburst-dominated ISM). These comprehensive and unbiased line surveys, spanning a continuous bandwidth, are essential for detecting multiple transitions of molecules beyond CO. They allow us to study the ISM by probing a large range of ISM conditions, provide rich and robust ISM diagnostic tool sets, improve our current potentially biased understanding of the physics and chemistry of the ISM, and reveal how various extreme ISM environments influence its properties.

Throughout this work, we assume a Chabrier (2003) stellar initial mass function (IMF) and a definition of infrared luminosity integrated over 8–1000 μm. Based on the starburst evolutionary synthesis models in Leitherer & Heckman (1995) with a Salpeter (1955) IMF, Kennicutt (1998) derived the star formation rate (SFR) to infrared luminosity, SFR-L_IR, calibration for starburst galaxies (with continuous bursts of age ∼100 Myr) as SFR = 1.7 × 10⁻¹⁰ (L_IR/L_⊙) M_⊙ yr⁻¹. After translating this calibration to the Chabrier (2003) IMF using a factor of ∼0.63 for the SFR calibration (Madau & Dickinson 2014), we derive SFR = 1.1 × 10⁻¹⁰ (L_IR/L_⊙) M_⊙ yr⁻¹ (with ≈30% uncertainty due to variations in star formation histories within individual sources; Kennicutt 1998). We note that using a top-heavy stellar IMF, which is found in dusty starburst galaxies, may reduce the SFR/L_IR ratio by a factor of at least ∼2 (Cai et al. 2020) and up to ∼5 (Yan et al. 2017; Zhang et al. 2018). Additionally, if the burst time is shorter by one order of magnitude, the calibration factor will increase by a factor of 3 (Madau & Dickinson 2014), introducing additional uncertainties to the SFR-L_IR calibration.

We adopt a spatially flat Lambda cold dark matter (ΛCDM) cosmology model with H₀ = (67.4 ± 0.5) km s⁻¹ Mpc⁻¹ and Ω_M = 0.315 ± 0.007 (Planck Collaboration VI 2020). This translates to luminosity distances and scales of 35 700 ± 370 Mpc and 7.44 kpc/″ for APM 08279+5255 and 31 970 ± 320 Mpc and 7.17 kpc/″ for NCv1.143.

This paper is organized as follows: Sect. 2 introduces the targets of the surveys, Sect. 3 describes the NOEMA observations and the process of data reduction, and Sect. 4 presents a detailed analysis of the data, including line identification. Section 5 discusses the continuum emission covered by our line survey. Section 6 delves into a detailed discussion about the line survey results and a comparison between the two sources, which is followed by chemical modeling in Sect. 7. The impact of differential lensing is addressed in Sect. 8. We conclude in Sect. 9.

2. Targets – APM 08279+5255 and NCv.1.143

APM 08279+5255 is a z = 3.911 strongly lensed broad absorption line quasar with complex multiple absorption features mainly seen in the C IV line (Irwin et al. 1998) and a prodigious apparent bolometric luminosity of 5 × 10¹⁵ L_⊙ (Egami et al. 2000). It has been suggested that radiation pressure plays an important role in producing quasar outflow (Saez & Chartas 2011). Observations of Hβ and Mg II suggest a mass of the central supermassive black hole (SMBH) of ∼10¹⁰ M_⊙ (Saturni et al. 2018). APM 08279+5255 is one of the brightest dust- and gas-rich quasars in the far-infrared and submillimeter (e.g., Downes et al. 1999; Egami et al. 2000; Weiß et al. 2007). It is extremely dusty with a total infrared luminosity L_IR = (3.4 ± 0.4)×10¹⁴μ⁻¹ L_⊙ (Leung et al. 2019, note the different definition of L_IR¹). The molecular gas is found to be highly excited with extreme conditions by examining the CO spectral line energy distribution (SLED). The CO excitation can be attributed to two prominent excitation components – a cold one with T_kin ∼ 65 K and n_H2 ∼ 10⁵ cm⁻³, and a warm one with T_kin ∼ 220 K and n_H2 ∼ 10⁴ cm⁻³ which is probably heated by the central AGN (Weiß et al. 2007). Observations of multiple H₂O lines in APM 08279+5255 suggest that the submillimeter H₂O emission originates in similar regions as CO with a dust temperature T_d ∼ 220 K and the optical depth at 100 μm, τ₁₀₀ ≈ 0.9, where the X-ray radiation of the quasar nucleus penetrates deeply into the buried regions and contributes to the ISM heating (van der Werf et al. 2010; Lis et al. 2011). From the early optical and near/mid-infrared observations, APM 08279+5255 shows a triple-image structure, likely composed of the lensed images of the quasar (Oya et al. 2013). Based on these observations, APM 08279+5255 is believed to exhibit a very high lensing magnification (μ ∼ 90–120). This magnification is notably sensitive to the size of the source within the source plane, showing a differential lensing factor of up to 10 (Ledoux et al. 1998; Egami et al. 2000). However, later analysis of the 0.3″ resolution CO(1 − 0) Very Large Array (VLA) image suggests a naked cusp lens model with a significantly lower magnification of ∼7 (Lewis et al. 2002). Moreover, a magnification of μ = 2–4 was derived based on a different lens model based on a deeper 0.3″ CO(1 − 0) image (Riechers et al. 2009), suggesting that the CO disk has an apparent (magnified) radius of about 1.3 kpc, consistent with the CO size derived from excitation analysis (Weiß et al. 2007). In this lens model, differential lensing is insignificant (variations ≲20% for source radius smaller than 5 kpc) and is at most a factor of 2 for the most extreme intrinsic source radius that exceeds 20 kpc. Throughout this work, we adopt this latest lens model.

The dusty starburst galaxy NCv1.143 is among the brightest strongly lensed sources discovered from the Herschel-ATLAS survey (or H-ATLAS for short; Eales et al. 2010); its IAU name is H-ATLAS J125632.7+233625 (or HerBS-5, Bakx et al. 2018). At z = 3.565, NCv1.143 is an intrinsically hyper-luminous infrared galaxy with an apparent total infrared luminosity¹ of L_IR = (1.4 ± 0.2)×10¹⁴μ⁻¹ L_⊙. The magnification is derived to be μ = 11.3 ± 0.7 based on the 880 μm continuum from Bussmann et al. (2013) and μ = 12.2 ± 1.2 based on the 2 mm continuum from Yang (2017). Throughout this work, we have adopted the latter value derived from the lens model described in Appendix C. Thus, the intrinsic infrared luminosity is ≳10¹³ L_⊙ with a molecular gas reservoir of (6.2 ± 1.6)×10¹⁰ M_⊙ (Yang et al. 2017). The source has one of the highest L_IR surface densities in the H-ATLAS sample (Bussmann et al. 2013), reaching Σ_LIR ∼ ×10¹³ L_⊙ kpc⁻². A series of extensive follow-up studies (e.g., Yang et al. 2016, 2017) found no clear AGN signatures, indicating that this galaxy is likely dominated by starburst activity with an estimated SFR = 1.3 × 10³ M_⊙ yr⁻¹, although a deeply buried AGN cannot be completely ruled out. Studies of the molecular gas conditions using the CO SLED in NCv1.143 found two prominent excitation components of molecular gas – a lower-excitation one with a kinematic temperature T_kin ∼ 20 K and a molecular gas density n_H2 ∼ 10^4.1 cm⁻³, and a high-excitation molecular gas component with T_kin ∼ 63 K and n_H2 ∼ 10^4.2 cm⁻³ (Yang et al. 2017). Detection of the J_up = 2, 3 and 4 H₂O lines in NCv1.143 suggests the existence of an optically thick (τ₁₀₀ >  1) dusty nucleus with a dust temperature T_d ∼ 75 K, surrounded by an extended cooler disk with T_d ∼ 35 K (Yang et al. 2016; Yang 2017). The similarities between the line profiles of the high-J CO and H₂O lines indicate that they both arise from similarly dense warm gas in star-forming regions across the galaxy.

In addition to CO and H₂O, detections of dense gas tracers such as HCN, HNC, and HCO⁺ in APM 08279+5255 (Wagg et al. 2005; García-Burillo et al. 2006; Weiß et al. 2007; Riechers et al. 2010) and NCv1.143 (Yang 2017) have also been reported. These detections, however, lean toward the brightest lines, focusing on specific frequency windows and often capturing only one transition per molecule. This limited scope presents a challenge for conducting a thorough examination of the ISM and its interaction with radiation fields, as such an analysis necessitates the consideration of complex astrochemistry.

3. Observations and data reduction

The deep 3 mm spectral line surveys of APM 08279+5255 and NCv1.143 were conducted using NOEMA with its PolyFix correlator (Gentaz 2019). The observations were carried out under projects W18EB and S18DC (PIs: C. Yang & A. Omont). The PolyFix correlator can process a total instantaneous bandwidth of 7.744 GHz for up to 12 antennas per polarization (dual polarization) per sideband – namely the upper sideband (USB) and the lower sideband (LSB) – which provides a frequency coverage of 15.488 GHz per tuning in dual-polarization, with a fixed channel spacing of 2 MHz across LSB and USB, separated by a gap of 7.744 GHz. As shown in Fig. 1 and Table 1, we designed the tunings to allow for a continuous frequency coverage from 71 to 109.5 GHz for NCv1.143 and from 71 to 111 GHz for APM 08279+5255 (except for a gap between 101.7 GHz and 103.3 GHz, intentionally placed to cover some crucial lines at the high-frequency end, close to 111 GHz). In the rest frame, the frequency coverage for APM 08279+5255 corresponds to 347.71–500.76 GHz and 507.28–546.26 GHz, and for NCv1.143, it corresponds to 323.15–500.73 GHz. Therefore, the total frequency coverage for both sources extends to approximately 200 GHz in the rest frame. This extensive coverage, as intended, enables us to detect multiple transitions of some key molecular emitters.

Fig. 1. Tuning setups of the NOEMA observations. Top panel: Atmospheric transmission curve at the zenith of the NOEMA site for the frequency range covered by our line survey. The red and blue lines show the transmission curves for PWV = 1 mm and PWV = 5 mm, respectively. Middle and bottom panels: Tuning setups, their sky (observed) frequency coverage, and the RMS noise reached in each spectral window for APM 08279+5255 and NCv1.143. For display purposes, we smoothed the spectral resolution to 50 MHz, resulting in an RMS that is five times smaller than the noise level reported in Table 1 for the 2 MHz bins. As expected, the edges of the sidebands and the conjunction channels of the two spectral windows within each sideband exhibit large RMS values.

As summarized in Table 1, the observations spanned two semesters from June 2018 to March 2019. Since our aim is to detect weak spectral features, we conducted our observations in compact C or D configurations to maximize sensitivity. The number of antennas used ranged from six to ten, with ten being used most of the time (Table 1). The baselines remain compact, with lengths from 24 to 176 m (D configuration) for most of the observations, while for a small portion of the observations, we also used baselines up to 368 m (C configuration). The resulting synthesized beam sizes with natural weightings have modest/low resolutions of ∼1.7″×3.0″ to 4.2″×5.6 ″(about 12×22 kpc to 31×41 kpc). The pointing positions are given in Table 1. We note that, for APM 08279+5255, there is an offset of about 9″ between the phase center and the peak flux position (RA 08h31m41.7s Dec +52°45′17.5′′; see Fig. 2), which does not have any significant impact in our observation, as the position shift can be neglected considering the size of the primary beam (FWHM ∼50″ at 3 mm). Most observations were carried out in good or excellent weather conditions, except for those on November 22 and December 13, 2018, when the weather conditions were moderate. The observations of NCv1.143 were performed mostly in the summer semester and partially in the winter semester, with a mean phase RMS of ∼15° and precipitable water vapor (PWV) below 7mm, while APM 08279+5255 was observed in the winter semester with a mean phase RMS of ∼5° and PWV below 5mm. Consequently, the total data loss due to flagging is below 20%. The total on-source time for all the tunings was ∼22.5 h for NCv1.143 and ∼16.4 h for APM 08279+5255. In the end, the flux RMS levels in the native 2 MHz spectral resolution reach ∼0.6–1.6 mJy beam⁻¹ and ∼0.6–1.0 mJy beam⁻¹ for NCv.143 and APM 08279+5255, respectively. The RMS is also relatively stable across all the frequency coverage, except for the edge of the spectral windows where the values are high (Fig. 1). Details of the setups and conditions of the observations are listed in Table 1 and Fig. 1. The phase and bandpass were calibrated by measuring standard calibrators that are regularly monitored at NOEMA, including 3C279, 3C273, MWC349, and 0923+392. The accuracy of the absolute flux calibration is estimated to be about ∼10% in the 3 mm band. The final uv tables containing the visibilities were produced after calibration using the GILDAS² packages CLIC, with the native 2 MHz spectral resolution. Then imaging, CLEANing, uv fitting, and spectra extraction were performed with GILDAS’s MAPPING on the uv tables.

Fig. 2. NOEMA CLEANed images of APM 08279+5255 and NCv1.143, made from the visibilities obtained by averaging all the USB channels from the highest-frequency tunings (w18eb003 and s18dc003; Table 1). These tunings were chosen since they are expected to offer the highest spatial resolution among all the spectral windows in our tunings. The contour levels begin at 5σ and increase in increments of 10σ. The unit of the color bar and contours is arbitrary, with values representing fluxes integrated over the entire USB of the tuning. The primary objective here is to illustrate that the sources are unresolved rather than provide precise flux measurements. The offset for APM 08279+5255 is explained in Sect. 3.

4. Data analysis

Upon acquiring the calibrated uv tables, we collapsed all the channels from each sideband per tuning. Subsequently, we imaged and CLEANed those collapsed uv tables to validate that the sources are not spatially resolved by the synthesized beams. As an example, we show the CLEANed images from the USB of the highest frequency tunings (s18dc003 and w18eb003, which correspond to the tunings with the highest spatial resolution) for NCv1.143 and APM 08279+5255 in Fig. 2. From the perspective of the all-channel-collapsed images, it is clear that both sources are not significantly resolved by the synthesis beams with dimensions of approximately 3.0″×1.7″–5.6″×4.2″ (using natural-weighting) for the signal surpassing the 5 σ levels.

However, since the all-channel-collapsed images are dominated by the dust continuum, it is important to check if any of our observed line emissions can be significantly more extended than the dust continuum. Notably, some molecular line maps (e.g., low-J CO) show a more extended spatial distribution than the dust continuum in dusty star-forming galaxies (e.g., Ivison et al. 2011; Calistro Rivera et al. 2018; Tsukui et al. 2023). After further examining specific channels where the emission is dominated by lines (e.g., CO), we find no significant difference in the size of the source between the all-channel-collapsed images and the line-dominated channel images (all the differences are below 3σ). We thus conclude that all the line emissions and the dust continuum are not resolved by our beams. Additionally, we do not expect any significant missing flux, given the compactness of the emission with regard to the beam sizes.

Since the sources are not resolved and remain compact for both the emission from the continuum and the lines, with the purpose of extracting the total fluxes and maximizing the signal-to-noise ratios, the emission of our sources can thus be approximated using simple elliptical Gaussian models in the uv plane across all the channels. Such an assumption is also consistent with other long-baseline interferometric data of the sources, namely the 0.5″×0.3″ resolution 1.2 mm dust continuum, CO and H₂O line images of NCv1.143 (see Yang 2017 nad our Appendix C) and the 0.3″ resolution 2.6 mm dust continuum and CO images of APM 08279+5255 (Riechers et al. 2009). These high-angular-resolution data show that the largest structures of the ISM emission are ≲1″ for APM 08279+5255 and ≲2″ for NCv1.143. These structures are characterized by the lensing-produced multiple image components distributed along the Einstein rings. The sizes of the Einstein rings are slightly smaller than the synthesized beam sizes as listed in Table 1. The individual image components are about ≲0.2″ and 0.7″ for APM 08279+5255 and NCv1.143, respectively.

4.1. Extracting spectra

For flux extraction, we assume a model of a single two-dimensional elliptical Gaussians and directly fit the visibilities in the uv plane. We verified our fitting results by checking the residual images and did not find any significant residual fluxes. Nevertheless, we caution that our purpose is only to extract the total fluxes in each channel, and the study of the ISM morphology is beyond the scope of this work. Therefore, while the two-dimensional elliptical Gaussian models are accurate enough as approximations to extract source flux and reduce the ambiguity in CLEANing and aperture selection in the traditional method of spectra extraction, they do not reflect the source morphology (see Appendix C for discussion of the morphology of the sources with high spatial resolution data). Considering the spatial resolution of our line survey data and the main goal of maximizing spectral sensitivity and extracting fluxes, two-dimensional elliptical Gaussian is thus a good approximation.

The two-dimensional elliptical Gaussian profiles are characterized by the central positions relative to the phase centers, δ_RA and δ_Dec, the major and minor axes, r_major and r_minor, the position angle PA, and the total flux. We used the UVFIT package of the GILDAS to fit the visibilities in the uv-plane directly with the model. We fixed the parameter of δ_RA, δ_Dec, r_major, r_minor, and PA to the values derived from the all-channel-collapsed uv tables for each tuning (assuming that the parameters, except for the flux, do not vary significantly between the LSB and USB within the same tuning), and only allow the flux to change across all the channels. Accordingly, the fitted fluxes, along with the errors across different channels, produce the final spectra containing fluxes and their errors.

The central positions, which are part of the fixed parameters derived in each tuning, are listed in Table 2. For NCv1.143, we find an elliptical Gaussian size of ∼(1.3–1.6)″×(1.2–1.4)″, depending on the frequencies. For APM 08279+5255, the values of r_major/r_minor are from 1.41 ± 0.38 to 2.60 ± 1.60, indicating its elongated morphology, with similar PA values derived from 40 ± 10 to 46 ± 21 deg across three tunings. This is consistent with the elongated distribution seen at high-spatial-resolution images (Riechers et al. 2009). For NCv1.143, the aspect ratio (r_major/r_minor) is from 1.03 ± 0.06 to 1.40 ± 0.09, showing little degree of deviation from circular symmetry, consistent with high-angular-resolution images (Appendix C). Both uv-plane models are robust; the fitted parameters have a good agreement across all three tunings for each source.

The positions of APM 08279+5255 are particularly well constrained because of the stable observation conditions and the compactness of the source compared with the synthesis beams (Table 2). Nevertheless, we also find that the values of r_major and r_minor of APM 08279+5255 are decreasing with increasing frequencies, suggesting that the true source size is significantly smaller than the synthesis beam. While for NCv1.143, similar values of r_major and r_minor across three tunings indicate that the source is only slightly smaller than the synthesis beam.

Fig. 3. NOEMA 3 mm band spectral surveys of APM 08279+5255 and NCv1.143, shown as the blue and violet spectra, respectively. For visualization purposes, the spectrum of NCv1.143 has been shifted up by 2 mJy. Both spectra are binned to 50 MHz (about 150 km s−1) to highlight the line detection. Error bars are indicated by thin lines (the same color as the spectra) overlaid on the spectra. Dashed black lines identify the lines detected in both sources, while dashed orange lines highlight the lines that are only detected in one of the two sources. The line names are labeled above the dashed lines. The peak values indicated for CO(3 − 2), CO(4 − 3), and [C I](1-0) emission lines in NCv1.143 are the true peak values before the 2 mJy shift, and they are outside the box of the figure. A zoomed-in view is provided in Fig. A.1.

From the aforementioned method, we obtained six spectra for each source corresponding to all the USB and LSB of the three tunings per source. The spectra are then combined through a simple linear re-grid along the frequency axis, taking the errors of each channel into account. The final combined spectra are displayed in Fig. 3 (a zoom-in view is provided in Fig. A.1). The following analyses throughout the work are then based on these final combined spectra.

4.2. Line identification

To maximize the accuracy of the line identification, we adopt two complementary methods that cross-verify each other to produce the final line catalogs. First, we utilized MADCUBA³ (MAdrid Data CUBe Analysis; Martín et al. 2019a) to identify molecular species by simultaneously fitting all their transitions within our frequency coverage, also accounting the brightness of the lines as predicted by the local thermal equilibrium (LTE) assumption. Second, we input the identified line detection list, as priors, into our matched filter identification code to cross-check and ensure that no line features are overlooked. In the following sections, we provide the details of these processes.

4.2.1. Local thermal equilibrium analysis and line identification

As an initial step to identify the spectral line detections in our observations, we deployed MADCUBA to conduct line fittings of all the lines across the entire spectra, assuming LTE conditions. Fitting the emission of all transitions for all the possible species, rather than identifying individual spectral features by their central observed frequency, allows for a robust identification when multiple transitions are detected. Additionally, this method accounts for line blending based on identified species and their LTE-predicted fluxes. We included a total of 29 species in the model, drawing from previously published extragalactic line surveys (Martín et al. 2021). A list of the molecules and transitions detected is compiled in Table 3. Within the 29-species list, undetected species included in the fit are CH₃CN, CH₃OH, H₂CO, H₂S, N₂D⁺, the isotopologs C¹⁷O, DNC, C¹⁵N, HC¹⁵N, H¹⁵NC, and the vibrational states of HNC. We do not identify any other feature outside this list in our spectra (which is cross-checked in Sect. 4.2.2).

The model adopts a linewidth of 520 and 330 km s⁻¹ for APM 08279+5255 and NCv1.143, respectively, as derived from CO emission from our data. These values are consistent with previous measurements (Weiß et al. 2007; Yang et al. 2017). Here, we assume an apparent source size of 0.13 arcsec², which should be within a factor of ≲3 from the sizes of both sources (with significant uncertainties; see our Appendix C and Riechers et al. 2009). The effect of the source size in the model partially degenerates with the column density, where a smaller source size would require higher column densities to model the observed spectra. Given that the bulk of our analysis is rooted in comparing relative column densities (to that of C¹⁸O, see Sect. 6), the uncertainties associated with size estimations are likely canceled out if the line-emitting regions are similar. If this is the case, the choice of source size will have minimal influence on our overall conclusions. We also acknowledge that some of the lines might have different sizes than that of the C¹⁸O lines, such as more extended low-J CO. However, because of the limited spatial resolution of our data, future high-angular-resolution observations are needed to disentangle the uncertainties in the sizes of the emitting regions for a more accurate accounting of the size variations across all the lines. With the assumed source size, the models are still within the optically thin regime for almost all the lines. A much smaller source size would result in line saturation, mostly on the CO transitions toward NCv1.143, while the rest of the species is still in the optically thin regime. For those species for which the excitation temperature could not be constrained with the observations, mainly those with a single spectral line feature, a T_ex of 36 and 25 K has been derived from HCN, HCO⁺, and HNC and assumed for all other species (except for H₂O and the vibration excited HCN), for APM 08279+5255 and NCv1.143, respectively. It is important to note that H₂O lines are mainly powered by radiative pumping in both sources (van der Werf et al. 2010; Yang et al. 2016); therefore, the fitting here (assuming a molecular-gas temperature of 230 K for both sources) serves only to identify the H₂O features. We assume an ad hoc high temperature of 300 K, which is the typical value found in Arp 220 (Martín et al. 2011) and NGC 4418 (Sakamoto et al. 2010), to fit the vibrational transitions of HCN given its high energy levels. Table 4 lists the fitted column density (N) and temperature of the molecular species under LTE for both sources. We note that the excitation temperature (T_ex) should be treated as a lower limit of the kinetic temperature (T_kin) of the molecular gas (Goldsmith & Langer 1999). Nevertheless, the higher T_ex values found in APM 08279+5255 as compared to NCv1.143 reflect that the molecular gas conditions are generally more extreme in the former.

4.2.2. Matched-filtering line identification and the final fitting results

After identifying the spectral lines using MADCUBA under the LTE assumptions as described in the previous section, we fed in the list of identified lines as priors and performed a second search based on a matched-filtering algorithm, where we convolve the spectra with Gaussian kernels based on the best fit CO linewidths, utilizing the specutils package of AstroPy (similar approach as in, e.g., Boogaard et al. 2023). This process leveraged the well-detected bright CO-line profile to enhance detection and filter out noise for better detecting weak emission features.

First, we conducted a “blind” matched-filter detection using only a 5-σ detection threshold. Then, we further lowered our detection threshold based on the prior information about the line detection from MADCUBA, especially for the case of the blended lines, and performed a second line search. Accordingly, we identified all the emission and absorption lines after this search. Statistical tests are performed to choose the adequate threshold to make sure more than 99.7% of the detected lines are real. In the end, we also verified all the line identifications visually as a triple-check. We do not see any significant unidentified lines beyond the initial MADCUBA list. Thus, our initial list of 29 species is sufficient for line identification.

The final cross-verified detections are summarized in Table 3. During this process, we simultaneously fitted all the line features, along with the underlying continuum emission (Figs. A.2 and A.3). The discussion on the physical aspects of the continuum for both sources is presented in Sect. 5. In the fitting process, we also deployed a Markov chain Monte Carlo (MCMC) technique, where we used samplers that are twice the number of the free parameters with 2000 interactions after 1000 burn-in steps. The derived posterior distribution of the fitted fluxes can better account for the overall uncertainties. The resulting line fluxes are listed in Table 5.

5. Continuum: Thermal dust and free-free emission

Due to the exquisite capability of NOEMA to calibrate the continuum over large frequency ranges, we can constrain the continuum emission of the two sources using the full coverage of the continuum flux densities across the rest-frame frequency range of ∼330–550 GHz.

For both galaxies, we assume two components that contribute to the continuum at this frequency range: the Rayleigh-Jeans tail of the modified black body emission from the dust and the free-free (thermal Bremsstrahlung) emission produced by free electrons scattering off ions commonly in HII regions. Our assumption aligns with the global spectral energy distribution (SED) analysis of these two sources (Leung et al. 2019; Yang et al. 2017), where the contribution from the nonthermal part is negligible around ∼330–550 GHz. The full SEDs of the two galaxies align with the typical synchrotron emission in normal galaxies (Condon 1992), contrasting with radio-loud sources where nonthermal emissions can have significant contributions (e.g., Falkendal et al. 2019).

As shown in Fig. 4, to simplify the fitting, we used the Rayleigh-Jeans approximation for the thermal dust emission while adopting a α = −0.1 spectral index for the free-free; thus, we have the following expression of the total flux, which consists of two power-law components:

$$\begin{array}{c}\hfill C^{\mathrm{R}-\mathrm{J}}{({\nu }_{\mathrm{rest}}/{\nu }_0)}^{(2+\beta )}+C^{\mathrm{free}-\mathrm{free}}{({\nu }_{\mathrm{rest}}/{\nu }_0)}^{-0.1}\end{array}$$(1)

Fig. 4. Fitting of the continuum of APM 08279+5255 and NCv1.143. Fluxes and errors are in black and blue, respectively. The solid red lines represent the total flux from the model. The dotted and dashed orange lines are the free-free and thermal dust components from the model, respectively. It is clear from the fitting that APM 08279+5255 has a flatter spectrum, which can be explained by contributions from free-free and/or corona emission from the SMBH, while the continuum of NCv1.143 is dominated by thermal dust.

where C^R − J and C^{free − free} are the scaling factors of the fluxes that reflect relative contributions from each component, ν is the rest-frame frequency (ν₀ = 300 GHz) and β is the dust emissivity index. The fitting results are shown in Fig. 4. The fitted values for C^R − J and C^{free − free} are 0.38 ± 0.03 mJy and 0.00 ± 0.02 mJy for NCv1.143 and C^R − J = 0.15 ± 0.01 mJy and C^{free − free} = 0.38 ± 0.01 mJy for APM 08279+5255. Therefore, for NCv1.143, we find the continuum covered by our line survey is totally dominated by dust emission, with β = 1.71 ± 0.02, while for APM 08279+5255, the free-free contribution is more prominent with a contribution from about 57% at 350 GHz to 27% at 500 GHz, and the dust continuum contributing the rest 43% to 73% with β = 1.73 ± 0.03. While the thermal dust flux contribution of 0.27 mJy is in excellent agreement with Leung et al. (2019) for APM 08279+5255, our free-free component (0.37 mJy) is at least three times larger than previous studies of the global SED (Stacey et al. 2018; Leung et al. 2019). We note that the cosmic microwave background can affect the fitted spectral index in the Rayleigh-Jeans part depending on redshift and dust temperature. Following da Cunha et al. (2013), assuming the dust temperatures of ≳ 40 K (Yang et al. 2017; Leung et al. 2019), we estimate that the value of β may be slightly underestimated by about 9%, while this value can be even smaller in APM 08279+5255 due to a higher dust temperature (see also Zhang et al. 2016). This will further increase the fraction of the free-free contribution toward the low-frequency end.

This raises the question of what causes such a difference in the ratio between free-free and thermal dust emission. Comparing the two sources, while the dust emission remains similar, the free-free emission from APM 08279+5255 is significantly higher. Both sources have high apparent SFRs around 1.2 − 1.5 × 10⁴ μ⁻¹ M_⊙ yr⁻¹. Using the correlation between free-free luminosity and the SFR (Algera et al. 2022), assuming the typical electron temperature of the HII regions of 10⁴ K and a thermal fraction of 1 (neglecting synchrotron emission around 300–500 GHz), we derive a luminosity of free-free emission of about 1.7 × 10³¹ erg s⁻¹ Hz⁻¹ at 350 GHz, which translates to fluxes of about 0.06 mJy, about six times lower than the observed value. But this value is in agreement with the free-free contribution from the global SED fitting of APM 08279+5255 (Stacey et al. 2018; Leung et al. 2019). One possibility is that the contribution from the synchrotron emission in APM 08279+5255 is significant because of its AGN activities. However, longer wavelength data of APM 08279+5255 showed that the synchrotron emission is almost two orders of magnitude smaller than thermal dust and free-free, which rules out this scenario (Stacey et al. 2018; Leung et al. 2019). Another possibility is an additional contribution from very cold dust (T < 15 K) emission (typically found in nearby star-forming galaxies, Galliano et al. 2018) in APM 08279+5255. However, considering the similarity of the bulk of the dust emission in NCv1.143 and APM 08279+5255, it will be difficult to explain why such very cold dust is not present in NCv1.143. The third possibility is an additional millimeter flux contribution from the non-self-absorbed part of the synchrotron radiation from the hot corona around the SMBH (e.g., Laor & Behar 2008; Inoue & Doi 2014; Kawamuro et al. 2022). Accordingly, using the 2–10 keV X-ray luminosity of APM 08279+5255 (μL_X = 3 × 10⁴⁶ erg s⁻¹; Bertola et al. 2022) and the correlation between millimeter and X-ray corona emission (Kawamuro et al. 2022), we find a corresponding millimeter flux of about 0.05 mJy at 350 GHz (assuming the millimeter corona emission peaks around 300 GHz), with significant uncertainties. Such a millimeter emission from the corona component may explain the observed elevated fluxes of APM 08279+5255 at lower frequencies. If this is the case, the value of β could be even higher to account for the steeper increase with less contribution from the free-free after including this corona component. In the global SED of APM 08279+5255 shown in Stacey et al. (2018), it is evident that there is flux excess, which is about 0.1–0.2 mJy depending on the frequency, that cannot be explained by free-free combining thermal dust around rest frame 100–400 GHz. The excess value is in broad agreement with a 50% contribution from the millimeter corona flux inferred from the X-ray luminosity. However, we will need to further combine data at the rest-frame frequency range of 50–300 GHz, where the corona emission usually peaks, to test this scenario. This is beyond the scope of this work and will be presented in another paper (del Palacio et al., in prep.).

6. A detailed look at the detected molecules and a comparison of the two sources

As displayed in Fig. 3 and Table 3, we have detected a total of 17 species, including two isotopologs of CO and one isotopolog of CS, plus the vibrational line of the HCN. Except for the [C I] line, all of the species have at least two transitions detected in at least one of the sources. Among these lines, the CH, NO, HCN-VIB (rotational transitions within the vibrationally excited state of HCN v₂ = 1f), 380 GHz H₂O maser, and H₃O⁺ lines are the first high-redshift detections in individual sources reported.

Given the existing extensive discussion of the CO SLED in APM 08279+5255 and NCv1.143 (Weiß et al. 2007; Yang et al. 2017), we do not delve into the discussions of the CO lines in this work. Here, we only highlight that the CO SLED of APM 08279+5255 shows a more elevated shape toward the high-J CO lines in Fig. 5, indicating an overall more excited molecular gas condition (Rosenberg et al. 2015).

Fig. 5. Line luminosities of APM 08279+5255 and NCv1.143 (based on the line fluxes listed in Table 5, where we also include several detections from the literature). Error bars are propagated from the measured flux errors. The upper panel shows the line luminosity in K km s−1 pc2, $L_{\mathrm{line}}^{\prime }$, normalized by the value of C18O(4–3). The lower panel shows the line luminosity in L⊙, Lline, normalized by the C18O(4–3) line. The lowest levels of the transitions of each molecule detected are labeled at the top of each panel, except for H2O and H3O+, where the rest-frame frequencies in GHz are given. Each transition is identified in the upper part of the figure. Compared with NCv1.143, the excitation condition is more extreme in APM 08279+5255 as seen from the SLEDs from several molecules.

From the observed integrated fluxes in Table 3, we derive the apparent line luminosities, μL_line (in units of L_⊙) and $\mu L_{\mathrm{line}}^{\prime }$ (in units of K km s⁻¹ pc²) following Solomon et al. (1992) using $L_{\text{line}}=1.04\times {10}^{-3}{I}_{\text{line}}{\nu }_{\text{rest}}{(1+z)}^{-1}{D}_{\text{L}}^2$ and $L_{\mathrm{line}}^{\prime }=3.25\times {10}^7{I}_{\mathrm{line}}{\nu }_{\mathrm{obs}}^{-2}{(1+z)}^{-3}{D}_{\mathrm{L}}^2$, where I_line, ν and D_L are in the units of Jy km s⁻¹, GHz and Mpc, respectively. We note here that because the ¹³CO(4–3) and CS(9–8) lines are blended due to their close frequencies, we extrapolate the integrated flux of CS(9–8) by combining the data of CS(7–6), CS(8–7) and CS(10–9), and we subtracted this flux from the blended line flux for obtaining the integrated flux of ¹³CO(4–3). Considering the smooth CS SLED in both sources, we do not expect a strong overestimation of the CS(9–8) line flux, which can cause an underestimation of the integrated flux of the ¹³CO(4–3) line in APM 08279+5255.

Assuming that the magnification factors do not vary significantly among the emissions (Sect. 8), we remove the magnification factor by normalizing the luminosities using the C¹⁸O(4–3) line (given that the C¹⁸O line is likely optically thin as we argue in Appendix D and the emitting size of the C¹⁸O is very similar to most of the other dense gas tracers; see, e.g., Martín et al. 2021). In Fig. 5, we show the comparison between the line luminosities normalized by C¹⁸O(4–3) (L_line and $L_{\mathrm{line}}^{\prime }$) of APM 08279+5255 and NCv1.143. We caution that the ratio in Fig. 5 depends on the normalization factor, which, however, does not affect the overall shapes of the SLEDs.

As shown in Fig. 5, the luminosity of CO(4–3) and C¹⁸O(4–3) are similar for both sources, despite APM 08279+5255 having an overall elevated CO SLED. APM 08279+5255 has a significantly low luminosity ratio of [C I](1-0)/C¹⁸O(4–3) compared to NCv1.143, while the ¹³ CO(4 − 3)/C¹⁸O(4–3) ratio is higher than NCv1.143. We also see a similarly high ratio of the luminosities of other dense gas tracers, such as CCH, HCN, HCO⁺, HNC, and CS over the value of C¹⁸O(4–3). If there is no deficit of C¹⁸O in APM 08279+5255, then the lines of dense gas tracers are at least a factor of 2 brighter in APM 08279+5255, indicating that the gas conditions in the quasar host are more extreme. This is consistent with molecular gas conditions derived from the CO SLEDs of APM 08279+5255, where the thermal pressure reaches 2 × 10⁶ K cm⁻³ (Weiß et al. 2007). While in NCv1.143, the thermal pressure is ≲1 × 10⁶ K cm⁻³ (Yang et al. 2017). This is also in line with our non-LTE excitation analysis of the dense gas tracers presented in Sect. 6.5.

Fig. 6. Column density ratios to C18O in APM 08279+5255 and NCv1.143 derived by MADCUBA (Table 4), which can be treated as the relative abundances. The plot is ordered, from left to right, according to the relative abundance in APM 08279+5255 from high to low. Some local sources from the literature are also shown: NGC 253 (CMZ; Martín et al. 2021; Rangwala et al. 2014), NGC 1068 (inner ∼2 kpc; Aladro et al. 2015), NGC 4418 (Costagliola et al. 2015), and Arp 220 (Martín et al. 2011).

Before applying any gas excitation and radiative transfer model to the observed SLEDs, it is useful to inspect their shapes, where different heating mechanisms can leave different imprints (Rosenberg et al. 2015). In Fig. 5, the CO SLEDs reveal a clear difference at the high energy levels, where APM 08279+5255 shows almost constant $L_{\mathrm{line}}^{\prime }$ (thermalized), while the luminosity drops by a large factor for higher energy level CO lines in NCv1.143. This suggests that the gas conditions are more extreme in APM 08279+5255 than in NCv1.143. For the dense gas tracer lines of HCN, HCO⁺, HNC, and CS, APM 08279+5255 has slightly more elevated SLEDs than NCv1.143, showing a more extreme heating mechanism in the quasar host. The most striking and complex difference is observed in the H₂O lines, where we find slightly elevated H₂O(2₁₁ − 2₀₂) and H₂O(3₂₁ − 3₁₂) luminosity in NCv1.143, while the 448 GHz H₂O(4₂₃ − 3₃₀) is about nine times brighter in APM 08279+5255. Such a difference suggests that the J_up = 2 and 3 H₂O lines might arise from very different conditions compared to the region that the H₂O(4₂₃ − 3₃₀) traces. The latter might be tightly related to AGN-powered radiative pumping. Thus, this line is much more excited in APM 08279+5255 than in NCv1.143. This is further explained by the detailed analysis of H₂O excitation in Sect. 6.3. Besides H₂O, the H₃O⁺ line in APM 08279+5255 is also much enhanced. Additionally, we identify the first high-redshift detection of the HCN-VIB line – three HCN-VIB lines with J_up ranging from 4 to 6 are detected in APM 08279+5255, while they remain undetected in NCv1.143.

6.1. Comparison of the LTE-derived abundances

To understand the astrochemical differences between the two sources, we compare their relative abundances derived using MADCUBA, together with some local prototypical galaxies: the central molecular zone (CMZ) region of the starburst galaxy NGC 253 (Martín et al. 2021; Rangwala et al. 2014), the central ∼2 kpc nuclear region of the Seyfert 2 galaxy NGC 1068 (Aladro et al. 2015), globally integrated values of the luminous infrared galaxy with a deeply buried nuclei NGC 4418 (Costagliola et al. 2015) and the archetypal ultra-luminous infrared galaxy Arp 220 (Martín et al. 2011). We also note that all the line surveys mentioned above are performed at scales ∼1–2 kpc that is comparable with the size of our galaxies, and this can potentially minimize the chemical variations due to differences of spatial scales (Butterworth et al. 2022).

Fig. 7. Comparison of the abundance ratios of APM 08279+5255 (APM 08279 for short), NCv1.143, and three archetype local galaxies – the CMZ of NGC 253 (Martín et al. 2021; Rangwala et al. 2014), ∼2 kpc nuclear region of NGC 1068 (Aladro et al. 2015), and NGC 4418 (Costagliola et al. 2015) – for the species detected in both APM 08279+5255 and NCv1.143. The abundance ratios of the galaxies are computed by dividing the column densities of each molecule by that of C18O, and these ratios are then compared among the galaxies. The horizontal axis is ordered by the increasing value of the ratio APM 08279/NCv1.143 (in the case of C34S and HCN-VIB, where the value of NCv1.143 is missing, we use the ratio of APM 08279/NGC 253 for ordering). Circles are AGN-dominated (APM 08279, NGC 1068, and NGC 4418) over starburst-dominated (NCv1.143 and NGC 253) ratios, diamonds are comparisons between similar types, and crosses are starburst-dominated over AGN-dominated ratios. The upper panel shows the ratio of the two high-redshift galaxies over local galaxies, while the lower panel shows the ratio of AGN-dominated galaxies over starburst-dominated galaxies only.

Given that we do not have a direct estimate of the absolute abundance of each molecule, we take the column densities listed in Table 4, normalize all the values with that of the C¹⁸O (as a proxy of total H₂, it is optically thin as shown in Appendix D; besides, the size effect mentioned in Sect. 4.2.1 is minimum because the size of the C¹⁸O emission is close to other dense gas tracers as found in, e.g., Martín et al. 2021), as a proxy of the relative abundance ratio (Fig. 6). Such a ratio also removes the uncertainty of the unknown source size (which includes the lensing magnification) if we assume that all the molecules are well-mixed. Nevertheless, we acknowledge the uncertainties as the molecular abundances can vary by a large factor within a galaxy (e.g., García-Burillo et al. 2014; Huang et al. 2022). However, given that our observation is spatially unresolved, such an assumption can be used as a first-order approximation until high-spatial-resolution data are available.

The comparison of the relative abundance (with respect to C¹⁸O) of all the molecules between the two sources, plus four nearby galaxies, is shown in Fig. 6. The molecules are ranked, from left to right, with the decreasing order of the relative molecular abundances of APM 08279+5255. Except for the neutral carbon, almost all the molecules show higher relative abundances in the quasar APM 08279+5255 than in NCv1.143, where the difference becomes most significant in the dense gas tracers like HCN, HCO⁺, HNC, CS, and C³⁴S. Notably, the CN abundance of NGC 4418 reported by Costagliola et al. (2015) is exceptionally high compared with all other galaxies. We also note here that the relative abundances are quite uncertain, given many sources of the uncertainties involved in the LTE analysis.

To better understand the differences in the relative abundances between the two sources and compare these abundances with some local sources, we present the ratios of relative abundances for APM 08279+5255 and NCv1.143 in Fig. 7. In addition to the ratio between APM 08279+5255 and NCv1.143, we include the ratios for our high-redshift targets when compared with – the CMZ of NGC 253, which is driven by pure starburst activity (Martín et al. 2021); NGC 4418, characterized by AGN activity (Costagliola et al. 2015); and NGC 1068, the inner ∼2 kpc nuclear region, where the abundances are expected to be dominated by AGN (Aladro et al. 2015). In Fig. 7 we group the ratios between each two of these galaxies into three categories: (1) AGN/starburst (solid circles), which are ratios of APM 08279+5255/NCv1.143, APM 08279+5255/NGC 253, and NGC 1068/NCv1.143; (2) AGN/AGN and starburst/starburst (solid diamonds), which are ratios of APM 08279+5255/NGC 1068, APM 08279+5255/NGC 4418, and NCv1.143/NGC 253; and (3) starburst/AGN (crosses), which are NCv1.143/NGC 1068, NCv1.143/NGC 4418, and NGC 253/NGC 1068.

From left to right in Fig. 7, we rank the molecules in the increasing order of the relative abundance ratio of APM 08279+5255/NCv1.143. The figure highlights the fact that the dense gas tracers like HNC, HCN, HCO⁺, and CS are considerably more abundant in APM 08279+5255 than in NCv1.143, while the difference is less prominent in CN, ¹³CO and CCH. It is worth noting that in high-temperature environments, where AGN can heat up the dust and gas, providing adequate conditions for high-temperature chemistry, CN can be converted into HCN effectively (Harada et al. 2010), leading to a decrease in CN and an increase in HCN. This is fully consistent with the fact that the relative abundance ratio of HCN/CN is significantly higher in APM 08279+5255 than in NCv1.143. The relative abundance ratio APM 08279+5255/NGC 253 shows a similar trend as APM 08279+5255/NCv1.143 – abundances of the dense gas tracers are enhanced, with the most extreme enhancement seen in C³⁴S. However, when we examine the ratios of APM 08279+5255/NGC 1068 and APM 08279+5255/NGC 4418 (AGN/AGN ratios), we do not see significantly enhanced relative abundances of the dense gas. The values of the ratios of APM 08279+5255/NGC 1068 and APM 08279+5255/NGC 4418 (AGN/AGN) are smaller than APM 08279+5255/NCv1.143 and APM 08279+5255/NGC 253 (AGN/starburst), and they are close to or below unity. As opposed to the ratios of APM 08279+5255/NCv1.143 and APM 08279+5255/NGC 253, we do not see systematic trends of these AGN/AGN ratios.

Conversely, when examining the abundance of NCv1.143 compared to local AGN-dominated sources (starburst/AGN), values lower than unity are consistently observed, which aligns with the higher values noted in the ratios of APM 08279+5255 over NCv1.143 and NGC 253 (AGN/starburst). This is not surprising, as the ratios of NCv1.143/NGC 1068 and NCv1.143/NGC 4418 are the inverse of the AGN/starburst ratios. When comparing NCv1.143 with NGC 253, we see a similar trend as seen in the AGN/AGN ratios (APM 08279+5255/NGC 1068 and APM 08279+5255/NGC 4418), where the values are found distributed around unity, and they are also often higher than the starburst/AGN ratios.

While the absolute abundance ratios for all the molecules in both our targets remain uncertain, distinct trends can be found in Fig. 7. First, when examining each species, the AGN/starburst ratios (APM 08279+5255/NCv1.143 and APM 08279+5255/NGC 253) generally exhibit higher values that are mostly above unity, whereas the starburst/AGN ratios (NCv1.143/NGC 4418 and NCv1.143/NGC 1068) often have lower values that are mostly below unity; the similar pairs, AGN/AGN or starburst/starburst (APM 08279+5255/NGC 4418, APM 08279+5255/NGC 1068, and NCv1.143/NGC 253) ratios, fall in between. This indicates that the molecular abundances in APM 08279+5255 behave like those found in the AGN-dominated sources, while the abundances in NCv1.143 resemble the starburst chemistry. Second, the dense gas tracers such as HCN, HCO⁺, HNC, CS, N₂H⁺, and C³⁴S are generally more abundant among the molecules studied in AGN-dominated sources compared to other molecules as shown in the lower panel of Fig. 7. Similar trends of an enhanced abundance of these dense gas tracers relative to C¹⁸O have also been found in the AGN-dominated CMZ in comparison to the starburst regions in NGC 1068 (Nakajima et al. 2023).

The parallel between APM 08279+5255 and local AGN prototypes, coupled with the similarity in behavior of NCv1.143 to the starburst CMZ of NGC 253, suggests that the ISM chemistry in APM 08279+5255 is influenced by the central AGN, through either the high-temperature chemistry proposed by Harada et al. (2010) or X-ray irradiate chemistry (Meijerink et al. 2007). Correspondingly, the chemistry of the ISM in NCv1.143 likely bears a resemblance to the CMZ of NGC 253, and we do not see any significant evidence of any influence from a buried AGN in this source from the molecular abundance.

Fig. 8. Energy diagram of an H2O molecule. Black lines indicate the H2O lines detected from our line survey, while gray lines indicate the observed water transitions from other works (see Table 5 for references), with their frequencies indicated (in GHz). The H2O data are taken from the JPL Molecular Spectroscopy website https://spec.jpl.nasa.gov.

6.2. The [C I] line

Both galaxies are detected in the neutral carbon [C I] ³P₁–³P₀ fine structure line at 492 GHz ([C I](1-0) hereafter), with luminosities $\mu L_{\mathrm{C}\mathrm{I}]}^{\prime }$ of 3.55 × 10¹⁰ K km s⁻¹ pc² and 10.1 × 10¹⁰ K km s⁻¹ pc² for APM 08279+5255 and NCv1.143, respectively. The luminosity of the [C I] line is weaker in APM 08279+5255 compared with NCv1.143, which can be explained by the X-ray-dominated region (XDR) model when the column density of gas N_H is around 10²⁴ cm⁻² (Meijerink & Spaans 2005).

The [C I] emission lines allow us to derive the properties of the atomic carbon gas in these systems and compare the H₂ masses derived independently from the CO lines. We note that the reported values here are not corrected for magnification. We estimated the neutral carbon masses using Eq. (1) in Weiß et al. (2005), assuming a [C I] excitation temperature equal to T_ex = 30 K, which is close to the value in Walter et al. (2011), ⟨T_ex⟩ = 29.1 ± 6.3 K, and the mean temperature of ⟨T_ex⟩ = 25.6 ± 1.0 K found by Valentino et al. (2020):

$$\begin{array}{c}\hfill M_{[\mathrm{C}\mathrm{I}]}=5.706\times {10}^{-4}Q(T_{\mathrm{ex}})\frac{1}{3}{e}^{23.6/T_{\mathrm{ex}}}{L}_{[\mathrm{C}\mathrm{I}](1-0)}^{\prime },\end{array}$$(2)

where Q(T_ex)=1 + 3e^{−23.6 K/T_ex} + 5e^{−62.5 K/T_ex} is the partition function of [C I] and the result is expressed in units of M_⊙. The derived masses are μM_[C I] = (4.40 ± 0.17) × 10⁷ M_⊙ and (1.25 ± 0.05)×10⁸ M_⊙ for APM 08279+5255 and NCv1.143, respectively.

In the case of NCv1.143, Yang et al. (2017) reported the detection of the [C I](2-1) emission line with a line flux of 4.4 ± 0.9 Jy km s⁻¹. Assuming LTE and, under the condition that both [C I] lines are optically thin, the excitation temperature equals the kinetic temperature:

$$\begin{array}{c}\hfill T_{\mathrm{ex}}/\mathrm{K}=38.8/\ln (2.11/[L_{[\mathrm{C}\mathrm{I}](2-1)}^{\prime }/L_{[\mathrm{C}\mathrm{I}](1-0)}^{\prime }]).\end{array}$$(3)

Using the measured flux of the [C I](2-1) emission line yielding $\mu L_{[\mathrm{C}\mathrm{I}](2-1)}^{\prime }=(4.9\pm 1.0)\times {10}^{10}$ K km s⁻¹ pc² (Yang et al. 2017), we derive an excitation temperature of 26.3 K for NCv1.143, close to the above-adopted value and comparable to the excitation temperatures found in previous studies of high-redshift galaxies (e.g., Nesvadba et al. 2019; Valentino et al. 2020). As pointed out in Sect. 4.2.1, this value of T_ex is likely a lower limit of T_kin. This is consistent with the T_kin of 20–63 K derived from the large velocity gradient (LVG) model of the CO SLED, while the [C I] excitation temperature is lower than the dust temperature T_dust = 40 K in NCv1.143 (Yang et al. 2017). This effect has been pointed out by Papadopoulos et al. (2022) via a sample of 106 galaxies, that the excitation condition of the [C I] lines are strongly sub-thermal.

Adopting an atomic carbon abundance of 8.4 × 10⁻⁵ (Walter et al. 2011) and an excitation temperature of 30 K, we derive total molecular masses of μM_mol = (8.7 ± 0.3)×10¹⁰ M_⊙ and (2.5 ± 0.1)×10¹¹ M_⊙ for APM 08279+5255 and NCv1.143, respectively. In the case of NCv1.143, the value is significantly smaller than the apparent molecular gas mass (5.2 ± 1.1)×10¹¹ M_⊙ derived by Yang et al. (2017) and for APM 08279+5255 the molecular gas mass estimated by Riechers et al. (2009) of 5.3 × 10¹¹ M_⊙ is higher by a factor of ∼2.

Adopting the average $L_{[\mathrm{C}\mathrm{I}]}^{\prime }/M_{\mathrm{mol}}$ conversion factor by Dunne et al. (2022), α_[C I] = 17.0 M_⊙ (K km s⁻¹ pc²)⁻¹, the molecular gas masses would be μM_mol = (6.03 ± 0.24)×10¹¹ M_⊙ and (8.3 ± 1.0)×10¹¹ M_⊙ for APM 08279+5255 and NCv1.143, respectively, which are more comparable with the CO-estimated total molecular gas masses (Weiß et al. 2007; Yang et al. 2017).

6.3. Submillimeter H₂O lines

Thanks to the observations of extragalactic H₂O lines done by Herschel (Pilbratt et al. 2010) and ground-based telescopes (e.g., Omont et al. 2011, 2013; Yang et al. 2013, 2016), it is now established that H₂O is one of the most important interstellar molecules after CO. It is because the submillimeter H₂O lines are bright, capable of several key diagnostics of the ISM, and provide critical complementary information to the CO emission lines (e.g., González-Alfonso et al. 2010; van der Werf et al. 2011; Jarugula et al. 2019; Yang et al. 2020; Stanley et al. 2021; Pensabene et al. 2022; Decarli et al. 2023). This is due to the combination of the high abundance and large dipole of H₂O, as well as its rich energy-level structure and rotational excitation processes (Fig. 8), including infrared pumping in dusty starburst regions (González-Alfonso et al. 2022), and collisional excitation where the heating source is from post-shock gas and/or shock-front radiation heating in young stellar objects (van Dishoeck et al. 2021; Dutkowska & Kristensen 2023), and the warm dense gas regions powered by warm dust heated by massive stars and/or AGN (Liu et al. 2017). H₂O has a central role in the oxygen chemical networks in the ISM, where H₂O chemistry can be dominated by grain desorption, ionic reactions in typical conditions, and neutral-neutral reactions in high-temperature conditions (see the review of van Dishoeck et al. 2013, and references therein).

For APM 08279+5255 and NCv1.143, we explored their H₂O spectra at rest-frame frequencies below 550 GHz, where there are only three relatively strong H₂O lines available – H₂O(4₁₄ − 3₂₁), H₂O(4₂₃ − 3₃₀) and H₂O(5₃₃ − 4₄₀): all three are detected in APM 08279+5255 and two of them, H₂O(4₁₄ − 3₁₂) and H₂O(4₂₃ − 3₃₀), in NCv1.143 (Tables 3 and 5). This number may appear marginal compared to the typical strong submillimeter H₂O lines that have been observed at rest frequencies between 557 and 1410 GHz in local starburst galaxies such as Mrk 231 and Arp 220 (e.g., Pilbratt et al. 2010; Hayward et al. 2011; Yang et al. 2013). However, these water lines with E_up/k_B ∼ 320–730 K in APM 08279+5255 and NCv1.143 can provide great insights into the ISM conditions, as demonstrated in the following subsections.

6.3.1. The 380 GHz H₂O(4₁₄−3₂₁) maser line

Extragalactic water masers are powerful tools for measuring masses of SMBHs and can be used to derive angular distances to the host galaxies with high precision, thus constraining critical cosmological parameters such as the Hubble constant (Herrnstein et al. 1999; Lo 2005; Gao et al. 2016; Pesce et al. 2020). Despite many searches for H₂O megamasers at high redshifts, to date, only one has been reported. The 22 GHz H₂O megamaser was identified in the z = 2.64 lensed quasar MG J0414+0534 (Impellizzeri et al. 2008). The H₂O(4₁₄ − 3₂₁) line has also been detected in the same source (Kuo et al. 2019; Stacey et al. 2020). However, due to the complex broad line profile, it is not definitively confirmed whether this is a maser line or not.

Fig. 9. Spectra of the 380 GHz H2O(414 − 321) line in APM 08279+5255 and NCv1.143 in red (the green error bar indicates the averaged noise level around the spectra per bin), overlaid with the scaled-down profile of the CO(4 − 3) line. The 380 GHz H2O line appears narrower by a factor of a few compared with CO(4 − 3), indicating that it is likely associated with maser origins. Black arrows indicate a ∼2 σ satellite line detection, possibly from the H2O maser, which could originate from accretion disks around the central SMBH.

This ortho-H₂O(4₁₄ − 3₂₁) line at 380.197 GHz (380 GHz hereafter), with its metastable backbone level J_Ka, Kc = 4₁₄ (E_up = 324 K), is one of the most promising H₂O maser lines predicted by radiative transfer models (e.g., Neufeld & Melnick 1991; Gray et al. 2016). The detection of the 380 GHz H₂O in APM 08279+5255 and NCv1.143 (with signal/noise ratio, S/N ∼ 11 and 6, respectively – see Table 5) is one of the most important results of our NOEMA line survey. In both sources, unlike MG J0414+0534 (Kuo et al. 2019), the profiles of this line have about half the width of the CO and all the other emission lines, clearly indicating a maser emission line (Fig. 9). This is the highest-redshift maser detection to date. We note that the 380 GHz H₂O is the only submillimeter H₂O maser detected in both sources in the rest-frame frequency range of 330–550 GHz. Moreover, the 380 GHz H₂O maser is at least four times stronger than its 22 GHz H₂O(6₁₆ − 5₂₃) counterpart in APM 08279+5255 (undetected by Ivison 2006). The flux is significantly bright when juxtaposed with other possible submillimeter masers (at 355, 380, 437, 439, and 471 GHz; see Gray et al. 2016, with the exception of the 475 GHz line as described later in this section), as well as the undetected 22 GHz canonical H₂O maser line. It is thus one of the best extragalactic H₂O masers found so far and should be systematically searched for in high-redshift strongly lensed submillimeter galaxies.

Besides the narrow feature at the central velocity, we also identify tentative (∼2.0-σ) detections of narrow emission features close to the 380 GHz maser line in APM 08279+5255 (Fig. 9). These narrow high-velocity satellite lines of the 380 GHz maser are up to about −2000 km s⁻¹ away from the central velocity, and they are likely originated from maser spots on an accretion disk surrounding the SMBH, similar to what has been observed in local galaxies (Pesce et al. 2023). If this is true, their velocities are more than twice as large compared to the highest velocity maser spots of the 22 GHz H₂O line from a local survey of the H₂O megamasers disks (Kuo et al. 2011), indicating that the accretion disk, in which the 380 GHz H₂O masers are located, is very close to the central AGN or that the SMBH of APM 08279+5255 is much more massive than the local sources from Kuo et al. (2011). Taking scaling relations from local megamaser disks (e.g., Gao et al. 2017), such a high velocity indicates a mass of ∼10¹⁰ M_⊙ of the SMBH in APM 08279+5255, which is consistent with estimates using other methods (e.g., Trevese et al. 2013). However, we caution that the high-velocity satellite lines here are very tentative detections, and thus we have been carrying out ongoing follow-up observations with the Green Bank Telescope (GBT) to confirm if these features are real. A detailed analysis of the maser line, together with the follow-up GBT data, will be presented in a separate paper.

6.3.2. The 448 GHz H₂O(4₂₃ − 3₃₀)

The ortho-H₂O(4₂₃ − 3₃₀) line at 448.001 GHz is predominantly excited through the absorption of 79 and 132 μm photons in the transitions 4₂₃ ← 3₁₂ and 4₂₃ ← 4₁₄ (Fig. 8). The influence of collisional excitation on the 448 GHz H₂O line is generally deemed negligible due to the high energy levels involved (E_up/k_B > 400 K) and the high critical densities of these levels. The 4₂₃ level is thus likely mostly radiatively populated. Recently, this key H₂O line has been detected for the first time in the ISM of both local and high-redshift galaxies, where this optically thin line is found to be a good tracer of the deeply buried nuclei that have intense far-infrared fields (Pereira-Santaella et al. 2017; Yang et al. 2020; González-Alfonso et al. 2021). This makes the 448 GHz H₂O line particularly effective in tracing the elusive population of buried AGN in environments with extreme far-infrared radiation energy densities.

From Fig. 5, we can see that the 448 GHz H₂O line exhibits the most significant difference among all the water lines between APM 08279+5255 and NCv1.143, where this line is about nine times stronger in the quasar. This strongly indicates that there is a deeply buried radiation source with a large column of dust and gas in APM 08279+5255 that is nonetheless visible in the optical. Most importantly, the brighter 448 GHz H₂O emission suggests that the APM 08279+5255 quasar must have a significantly more intense far-infrared field buried in the nuclear region than NCv1.143. It is worth noticing that in APM 08279+5255, we have also detected several transitions of the HCN-VIB line, which is believed to be the tracer of the compact obscured nuclei (Aalto et al. 2015; Falstad et al. 2021), where the column density reaches beyond 10²⁴ cm⁻². Comparing the ratio between the 448 GHz H₂O with the HCN-VIB(3–2) line, we find a comparable flux ratio of about 6.1 in APM 08279+5255, which is close to the value of 7.9 ± 1.6 found in a local compact obscured nucleus of ESO 320-G030 (Falstad et al. 2021; González-Alfonso et al. 2021). Because APM 08279+5255 is visible in the soft X-ray (Bertola et al. 2022) as well as in the rest-frame UV and optical (e.g., Saturni et al. 2018), is it unclear how APM 08279+5255 is linked to the typical compact obscured nuclei in the local Universe featured by strong HCN-VIB emission but mostly no X-ray detection (Falstad et al. 2021). Future high-spatial-resolution data are needed to have a clearer picture of the ISM structure in order to answer this question.

6.3.3. The 475 GHz H₂O(5₃₃ − 4₄₀) line

The para-H₂O(5₃₃ − 4₄₀) at rest frame 474.689 GHz (475 GHz hereafter) has a very high upper energy level of E_up = 725 K, which is even higher than the commonly observed J_up = 5 H₂O line transition 5₂₃ − 5₁₄ at 1411 GHz (E_up = 642 K) in local galaxies (Yang et al. 2013). Observations of evolved stars found the 475 GHz H₂O lines are masers associated with fast outflows (Menten et al. 2008).

In dusty warm dense regions of galaxies, H₂O molecules can also absorb 53 μm far-infrared photons through 5₃₃ ← 4₂₂ (Fig. 8). Such absorption features have been detected in NGC 4418 (González-Alfonso et al. 2012). In this case, the high H₂O column density can efficiently populate the 4₂₂ level and thus sustain this far-infrared pumping process. Such a process may potentially provide a route for the formation of H₂O maser. Nevertheless, the observed line profile of H₂O(5₃₃ − 4₄₀) in APM 08279+5255 is similar to those of other dense gas tracers, suggesting the H₂O emission is from similar regions of the dense gas tracers that are tracing the bulk motion of the gas at large scales. This is contradictory to a maser origin of H₂O emission, where commonly, a specific physical condition has to be met (Gray et al. 2016 found a typical excitation condition with a kinetic temperature of ∼1000 K), leading to very localized emitting regions (thus very narrow linewidths), such as accretion disks surrounding SMBHs.

6.3.4. Model of the H₂O emission lines and its interpretation

Fig. 10. Best-fit model and the likelihood of the parameters from the model of the H2O lines and the dust continuum in APM 08279+5255. Top left: Observed lensing-corrected (black data points) versus modeled H2O integrated flux (green and blue are the extended and compact components, respectively; red is the total). Bottom-left: Lensing-corrected modeled dust SED from the two components (same colors as in the top panel). The photometric data are taken from Leung et al. (2019). Right panels: Likelihood distributions of the key parameters of the H2O model, with the same colors as in the left panels. The arrows indicate the maximum likelihood values (Table 6).

Fig. 11. Same as Fig. 10 but for NCv1.143. The photometric data for the dust SED are from Bussmann et al. (2013) and Yang et al. (2017).

Given that in both APM 08279+5255 and NCv1.143, there are detections of other H₂O lines from the literature (Lis et al. 2011; van der Werf et al. 2011; Yang et al. 2016; Yang 2017), it is worth combining all the detected H₂O lines together with the dust continuum to constrain the excitation conditions. For APM 08279+5255, there are eight H₂O lines from J_up = 2 to 5 included, while for NCv1.143, six H₂O lines are included (Table 5, Figs. 5 and 8). Here, we assume that the H₂O submillimeter excitation is dominated by the absorption of far-infrared dust-emitted photons, with additional contribution from collisional excitations. Then, we used the radiative transfer model carried out by González-Alfonso et al. (2014b); González-Alfonso et al. (2021) to constrain the physical properties of the dusty ISM. Here we adopt a magnification of μ ∼ 4 (Riechers et al. 2009) for APM 08279+5255 and μ ∼ 9 for NCv1.143 (Appendix C), and all the following H₂O analysis are performed after lensing correction.

The model assumes a number N_C of independent components, which are spherically symmetric sources with uniform physical properties, namely: the dust temperature T_dust, the continuum optical depth at 100 μm along a radial path τ₁₀₀, the column density of H₂O along a radial path N_H2O, the H₂ density n_H2, the velocity dispersion ΔV, and the gas temperature T_gas. The model components are classified into groups according to their physical parameters, each group covering a regular grid in the parameter space (T_dust, τ₁₀₀, N_H2O, n_H2). However, it should be noted that the conditions of the gas that the H₂O lines trace can be different from the CO-traced conditions. We keep fixed ΔV = 100 km s⁻¹ and T_gas = 150 K (which is the typical condition and the models are insensitivity to the choice of T_gas = 150 K, González-Alfonso et al. 2021). For both APM 08279+5255 and NCv1.143, N_C = 2 is required to obtain a satisfactory fit to the data, as shown below.

Radiative pumping by far-infrared radiation means that H₂O is sensitive to the dust SED, and thus, our goal is to fit the H₂O SLED and the far-infrared SED simultaneously. For this purpose, we use as input data of our minimization method the fluxes of the H₂O lines (and the upper limit of the H₂O,448 GHz line in NCv1.143) and a number N_cont of continuum flux densities. For APM 08279+5255 N_cont = 2 (at λ_rest ≈ 51 and 677 μm), while we use N_cont = 6 for NCv1.143. These continuum data are highlighted in magenta in the SEDs in Figs. 10 and 11. In our method, χ² is minimized for all possible combinations among the N_C components, giving for each one the squared radius R² of each component. The best-fit (fiducial) model is the combination that yields the minimum χ², and the regular grid enables the estimation of the likelihood of each parameter. The derived parameters – R, the H₂O abundance relative to H nuclei X_H2O, L_IR, and M_gas of each component – are also inferred. More details are given in González-Alfonso et al. (2021).

We first attempted to fit the H₂O and continuum emission with a single-model component (N_C = 1), but results were unreliable, with a best-reduced χ² value ${\chi }_{\text{red}}^2$ ≈ 4. This was indeed expected because the low-lying H₂O J_up = 2–3 lines are expected to arise in more extended regions than the J_up = 4 lines that trace buried regions (e.g., González-Alfonso et al. 2014a; Pereira-Santaella et al. 2017). A better fit was found with N_C = 2 components, with ${\chi }_{\text{red}}^2$ ≈ 0.9; in Figs. 10 and 11, the arrows indicate the best-fit values and solid histograms show their likelihood distributions.

In both galaxies, we find that the J_up ≥ 4 lines are formed in a warm “compact” component (R ∼ 330 pc and ∼150 pc in APM 08279+5255 and NCv1.143, respectively), which we identify with the nuclear region (nuclear core). The nuclear cores are moderately optically thick in the far-infrared with τ₁₀₀ ∼ 2, with warm dust temperatures of over 100 K (Figs. 10 and 11). In NCv1.143, this core has an intrinsic luminosity of L_IR ∼ 4 × 10¹² L_⊙, resulting in an extreme infrared luminosity surface density Σ_IR ≡ L_IR/(4π R²) = 1.4 × 10¹³ L_⊙ kpc⁻². This translates into a surface SFR value of SFR/(π R²)∼ 6.2 × 10³ M_⊙ yr⁻¹ kpc⁻² if the contribution to L_IR by a possible obscured AGN is negligible since no strong evidence of the presence of an AGN has been found in NCv1.143 (Yang et al. 2016). In APM 08279+5255, however, the compact nuclear core has an extreme luminosity of L_IR ∼ 4 × 10¹⁴ L_⊙ that results in Σ_IR = 3 × 10¹⁴ L_⊙ kpc⁻². Obviously, this can only be explained by a significant contribution of radiation from the AGN.

On the other hand, the J_up ≤ 3 lines in both galaxies, pumped by absorption of dust-emitted 101 and 75 μm photons, are formed in a more extended R ∼ 1 kpc starburst region with moderate T_dust ∼ 50 K, most likely the circumnuclear disk of the host. In NCv1.143, this size is comparable to the projected half-light effective radius traced by the H₂O(2₁₁ − 2₀₂), H₂O(3₂₁ − 3₁₂) and CO(10–9) line emission as presented in Appendix C (Fig. C.2), and also to the size of the averaged 870 μm dust continuum of submillimeter galaxies (SMGs, R ∼ 1 kpc; Gullberg et al. 2019). The circumnuclear disk in NCv1.143 dominates the luminosity with L_IR ∼ 6.5 × 10¹² L_⊙ and Σ_IR = 6 × 10¹¹ L_⊙ kpc⁻² (corresponding to ∼230 M_⊙ yr⁻¹ kpc⁻²).

We note that excluding the continuum flux density at 51 μm in the fit of APM 08279+5255 yields a very similar best-fit model, with the far-infrared and mid-infrared continuum emission well fitted. The estimate of very warm T_dust of the nuclear region is mainly driven by the H₂O(5₃₃ − 4₄₀) 475 GHz line. This means that the H₂O lines are indeed tracing the shape of the SED, and high columns of H₂O are present in the inner ∼0.5 kpc in APM 08279+5255 provided that these nuclear regions (still) have large reservoirs of gas of ∼10¹⁰ M_⊙.

It is also worth noting the similarities and differences between the two sources. Despite their different SEDs, the extended circumnuclear disks have similar characteristics in APM 08279+5255 and NCv1.143 (see Table 6). For the compact core, most of the parameters in Table 6 are also similar between APM 08279+5255 and NCv1.143. However, the dust temperature and infrared-luminosity of the core in APM 08279+5255 are significantly higher than those in NCv1.143. As shown in Figs. 10 and 11, the compact nuclear core in APM 08279+5255 is dominating the total infrared luminosity of the entire galaxy, while in NCv1.143, the extended circumnuclear disk has a much larger contribution to the total infrared luminosity. Therefore, as diagnosed by the H₂O lines, the main difference between the two sources is a much more prominent and powerful nuclear core component, indicating the “turn-on” phase of the extreme quasar in APM 08279+5255, which may be a short-lived phenomenon.

6.4. H₃O⁺ lines

As a symmetric top molecule, H₃O⁺ plays an important role in the oxygen chemistry network that is related to the H₂O formation in ion-neutral reactions (Gerin et al. 2016). First detections of the H₃O⁺ lines were achieved in two local starburst galaxy, M 82 and Arp 220, from the transition of p-H₃O⁺($3_2^{+}-2_2^{-}$) at 364.797 GHz, with E_up/k_B = 139 K (van der Tak et al. 2008, see also Sakamoto et al. 2021). H₃O⁺ has also been detected in the Cloverleaf quasar (Guélin et al. 2022). Another study of five local galaxies utilized the abundance ratios of H₃O⁺ with H₂O and HCO⁺ and found the chemistry in XDRs can explain some of the high column density ratios of H₃O⁺ over H₂O and HCO⁺ in some sources (Aalto et al. 2011).

Both the H₃O⁺($3_2^{+}-2_2^{-}$) and H₃O⁺($3_0^{+}-2_0^{-}$) lines are robustly detected in APM 08279+5255. These are the first high-redshift detections of any H₃O⁺ lines. While in NCv1.143, these two lines are below the noise level; therefore, we do not discuss NCv1.143 in detail here until further robust detections are reached. The LTE-derived column density of H₃O⁺ is about N_H3O⁺ ∼ 10¹⁵ cm⁻² after lensing correction. This value is comparable with NGC 253, NGC 1068, and NGC 6240, where the XDR and/or a high CR ionization rate in the photodissociation-dominated region (PDR) can explain the column density (Aalto et al. 2011).

The abundance ratio between H₃O⁺ and H₂O can also be used as a probe of ionization states. Using the H₂O column density derived from Sect. 6.3.4 of 2.5 × 10¹⁷ cm⁻², we find a relative abundance ratio of H₃O⁺/H₂O,∼ 3 × 10⁻³, which can be explained by both the PDR and the XDR (Riechers et al. 2009). Indeed, H₃O⁺ lines have also been detected in the pure-starburst NGC 253 (Martín et al. 2021; Holdship et al. 2022). While the aforementioned diagnostics are inconclusive, the most evident abundance ratio that suggests the XDR plays an important role comes from H₃O⁺/HCO⁺, where in APM 08279+5255 we find a very large abundance ratio of ∼ 40, which is difficult to explain by PDRs (Meijerink et al. 2007). This suggests that the XDR may significantly influence the abundances of HCO⁺ and H₃O⁺ abundances in APM 08279+5255, though the possibility of H₃O⁺ enhancement by CR cannot be entirely ruled out.

Fig. 12. Violin plots for the LVG-MCMC-derived posterior distributions of nH2, Tkin, and Nmol/dV in APM 08279+5255 and NCv1.143 (upper and lower panels, respectively). The white dot represents the median values, while the solid blue and dashed blue lines represent ±1σ and ±2σ, respectively. The violin shapes represent the posterior distributions of each parameter as presented in the corner plots (Figs. B.1 and B.2).

6.5. Typical high dipole moment tracers – HCN, HCO⁺, HNC, CS and C³⁴S – and non-LTE modeling

Because of their large dipole moment, the consensus is that the lines of HCN, HCO⁺, HNC, and CS are probably among the best tracers of dense gas where stars are forming. This is demonstrated by the linear correlation between the line luminosity of the dense gas tracers, especially HCN, and L_IR, which holds tightly over ten orders of magnitude, ranging all the way from dense cores of Giant Molecular Clouds to ultra-luminous infrared galaxies (e.g., Gao & Solomon 2004; Wu & Evans 2005; Gao et al. 2007; Zhang et al. 2014; Chen et al. 2015; Jiménez-Donaire et al. 2019; Zhou et al. 2022). Unlike the super-linear correlation found between the SFR and the total molecular gas traced by CO(1 − 0), the tight relation between the dense gas and the SFR suggests that the dense gas links more tightly with star formation rather than CO(1 − 0) (Gao & Solomon 2004). Therefore, the dense gas tracers such as HCN, HCO⁺, HNC, and CS serve as crucial probes for understanding the most intense star formation in the SMGs. Nevertheless, the detection of the standard dense gas tracers remains meager at high redshifts because the emission lines of these dense gas tracers are very weak compared to CO. HCN is only reported in a handful of individual luminous high-redshift galaxies and only detected in strongly lensed sources (Vanden Bout et al. 2004; Carilli et al. 2005; Gao et al. 2007; Riechers et al. 2007, 2011a; Carilli & Walter 2013; Oteo et al. 2017; Béthermin et al. 2018; Cañameras et al. 2021; Rybak et al. 2022, 2023). These lines provide the essential signature of massive starbursts, and it has been suggested that the dense gas fraction, as derived from the comparison between HCN and low-J CO lines, is higher in the high-redshift dusty star-forming galaxies (Oteo et al. 2017), though Rybak et al. (2022) argue that there is no difference in dense gas fractions.

In our NOEMA line survey, we have acquired the largest amount of multiple transitions of several dense gas tracers to date, including HCN, HCO⁺, HNC, CS, C³⁴S, NO, and N₂H⁺ simultaneously in single objects. This enables us to understand the gas excitation conditions by directly looking into each species; thus, we are able to perform excitation analysis via the LVG method (e.g., Goldreich & Kwan 1974) for APM 08279+5255, where at least three transitions are detected per molecule. For NCv1.143, using some strict priors based on the previous studies of this source, we can also get loose constraints of the excitation conditions (the lack of data constraints here is well reflected by the posteriors of the parameters, Fig. 12 and Appendix B), where only two transitions are detected for most of the dense gas tracers. For the purpose of comparing APM 08279+5255 and NCv1.143, we picked HCN, HCO⁺, HNC, CS, and C³⁴S for a detailed non-LTE analysis for the ISM conditions (N₂H⁺ is excluded in the LVG analysis due to the lack of collisional rates data).

Fig. 13. Fitting of the dense gas SLEDs using LVG models via the MCMC technique. Observed fluxes and their errors are indicated with black symbols. Orange lines show the best fit (Table 7), and light-yellow lines indicate 100 examples within the 15.87% and 84.14% quartiles of the posteriors.

Following the same implementation as described in Yang et al. (2017), we use a one-dimensional (1D) non-LTE radiative transfer code RADEX (van der Tak et al. 2007), assuming an escape probability of β = (1 − e^−τ)/τ from an expanding sphere geometry (Sobolev 1960) to model the observed flux. At the same time, we also adopted the MCMC code emcee as we fit the integrated line flux of the SLED for HCN, HCO⁺, HNC, and CS for both APM 08279+5255 and NCv1.143, as well as C³⁴S for APM 08279+5255. Here, we do not simultaneously fit all the tracers with a single condition because they do not trace the same gas conditions (Tielens 2005). The collisional rates are taken from the LAMDA database (Schöier et al. 2005). After deploying 100 walkers with 1000 interactions after 100 burn-in steps, we derive the posterior distributions of the kinetic temperature of the molecular gas (T_kin), the volume density of H₂ (n_H2), and the column density of each molecule unit velocity difference (N_mol/dV), and the solid angle (Ω_app) of the source where the lensing magnification factor is included.

Table 7 lists the results from the LVG-MCMC analysis, while the fitted SLEDs are shown in Fig. 13. The LVG models can well reproduce the observed integrated fluxes. However, it is worth noticing that the J_up = 5 HNC line is about 40% brighter than the LVG model (over 1 σ level). It is difficult to explain such a flux excess of the J_up = 5 HNC line. One possible cause is the blending by the CN(4–3) line. Nevertheless, our current data limit the ability to assess this possibility thoroughly, and we do not see a similar overestimate of the HNC(5–4) flux in NCv1.143. Thus, the source of the flux excess of the HNC(5–4) line remains elusive.

In order to facilitate a more effective comparison of the results from our LVG-MCMC analysis, we present the posterior distributions of the two sources using the violin plots in Fig. 12. It becomes apparent that all the n_H2 values derived from HCN, HCO⁺, HNC, CS, and C³⁴S significantly exceed that of the CO from previous studies (Weiß et al. 2007 and Yang et al. 2017, n_H2 ∼ 10³–10⁴ cm⁻³), falling within the range of 10⁵–10⁷ cm⁻³. This is not surprising as these mid-J lines of the dense gas tracers probe more extreme conditions of the high-density ISM than the CO lines. Within APM 08279+5255, the gas density traced by HCO⁺ is markedly lower compared to that traced by all other molecules. As for the kinetic temperature, T_kin, all tracers align to suggest a similar value of approximately 100 K. When examining the column density values derived from the LVG analysis, APM 08279+5255 displays well-defined values, with most of the species falling within the 10¹⁴–10¹⁵ cm⁻²/km s⁻¹ range. Notably, we observe that the column density of HCO⁺ is significantly below other species, particularly in APM 08279+5255.

When compared with the LVG-derived column densities and assuming a small dV, the abundances of APM 08279+5255 and NCv1.143 align with the MADCUBA LTE-derived values listed in Table 4. The only exception is the CO abundance of NCv1.143, where the LVG-derived value is one order of magnitude smaller in NCv1.143. The LVG-derived values confirm the trend we found in Fig. 7 that the relative abundance of the dense gas tracers seems all boosted in the quasar where AGN is dominated. In the next section, we take a closer look at how different ISM environments might affect the abundance of the dense gas tracers in the two sources.

6.6. AGN diagnostics with HCN, HCO⁺, HNC and CS

Here, we place the observed line flux ratios from high dipole moment tracers from our line surveys into what has been reported within the local Universe. A number of molecular line ratios have been proposed as tracers of the powering sources and/or physical conditions within the obscured central regions of bright galaxies.

The ratio between the HCN and HCO⁺ is a classic molecular discriminator of nuclear regions dominated by either a starburst or an AGN (e.g., Kohno et al. 2001; Krips et al. 2008; Imanishi et al. 2023). Observations of a statistically significant sample with low spatial resolutions have shown that an enhancement of the HCN/HCO⁺ can occur in galaxies hosting an AGN (Privon et al. 2015; Izumi et al. 2016; Imanishi et al. 2016). However, Privon et al. (2020) found no trends between HCN/HCO⁺ and the AGN fraction from hard X-ray observations. High-resolution imaging of nearby AGN galaxies has shown this enhancement to be located in the circumnuclear disk and not directly toward the AGN position (Viti et al. 2014; Martín et al. 2015), which is claimed to be the result of high-temperature chemistry due to mechanical heating in the surrounding of the AGN (Izumi et al. 2013; Harada et al. 2010). All these local studies show the complex nature of the HCN/HCO⁺ diagnostics.

Izumi et al. (2016) proposed a diagnostic diagram based on the ratios of HCN/HCO⁺(4–3) and HCN(4–3)/CS(7–6), which we cover in our observations. We derive HCN/HCO⁺ ratios of 1.4±0.2 and 1.6±0.2 toward APM 08279+5255 and NCv1.143, respectively. Similar ratios are obtained with the J = 5–4 and 6–5. The HCN/CS ratio is measured to be between 3 and 4⁴ and 6.7±3.3, for APM 08279+5255 and NCv1.143, respectively.

The ratios derived for the quasar APM 08279+5255 are not distinctive from what is observed toward starburst-dominated galaxies. This is not uncommon and matches what is observed in other observations of AGN at a low spatial resolution, which shows a wide range of ratios (Privon et al. 2015; Izumi et al. 2016; Imanishi et al. 2016). Conversely, the ratios observed toward the starburst NCv1.143 are very uncertain and might lie within the range where only galaxies with an AGN are found.

Although we could expect these ratios to be affected by continuum absorption or self-absorption from intervening material, the analysis from Imanishi et al. (2016) showed that line flux ratios are only moderately affected by such effects. Since we do not find any AGN footprint in NCv1.143 (Yang et al. 2016), especially from the absence of detection in radio and in the mid-infrared bands of the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), a buried AGN is unlikely to reside in NCv1.143. However, the sensitivity and the angular resolution of our data on the CS, HCN, and HCO⁺ lines are limited to reaching conclusions on whether they actually point to the presence of such an AGN, although we do see general trends of the abundance difference between AGN- and starburst-dominated sources as discussed in Sect. 6.1, and NCv1.143 is aligned with the non-AGN-dominated abundances.

The abundance ratios we derived for APM 08279+5255 and NCv1.143, expressed as N(HCN)/N(CS), are given in log base ten as $-0.7_{-1.1}^{+1.3}$ and $-0.1_{-2.7}^{+2.5}$, respectively. These correspond to the values of N(HCN)/N(CS) ratios of about 0.2 and 0.8 for APM 08279+5255 and NCv1.143, respectively, with large uncertainties. Compared with a detailed case study of the chemical networks in NGC 1068 (Butterworth et al. 2022), where AGN tends to enhance N(HCN)/N(CS), our values do not show such clear evidence of the enhancement. This is not surprising considering the large uncertainties of the abundance ratios and low spatial resolution of our measurements. The abundance ratios of N(HCN)/N(HCO⁺) for APM 08279+5255 and NCv1.143 are $2.3_{-1.3}^{+1.1}$ and $0.6_{-2.4}^{+2.0}$ in log base ten, which might show a significant abundance enhance of HCN over HCO⁺. This is consistent with what has been found in the AGN-dominated nuclear regions of NGC 1068, at 40 pc to 100 pc scales. If such similarities are valid, it also indicates that the spatial distribution of HCN and HCO⁺ are similar, as the case in NGC 1068 (Butterworth et al. 2022). Thus, their abundance ratio is a better representation of a consistent inference of their ISM condition compared to the abundance ratios of N(HCN)/N(CS).

6.7. CCH

The ethynyl radical, CCH, is abundant in the ISM and is known to be enhanced in PDRs, both from galactic and extragalactic observations, although ionizing processes such as CR ionizing could also be responsible for the CCH production (Martín et al. 2014; Holdship et al. 2021, see the discussion and references therein).

In this work, we observed line ratios of CCH over HCN, HCO⁺, and HNC using their J = 5–4 transitions of 0.29/0.37/0.46 and 0.29/0.46/0.48, toward APM 08279+5255 and NCv1.143, respectively. These ratios fall within the ratios presented by Martín et al. (2014) over a small sample of galaxies with a variety of environments. We do not see any particular flux enhancement of the CCH lines in our targets. However, when we compare the relative abundances from the LTE analysis directly, we find that CCH is more abundant in APM 08279+5255 than in NCv1.143. This is consistent with an enhancement of CCH abundance in the region dominated by AGN and shocks (Nakajima et al. 2023).

6.8. CH lines

The methydilyne radical, CH, is a major interstellar molecule both in the diffuse and dense ISM, as well as in XDRs (e.g., Gerin et al. 2016), typically tracing densities n_H of 100–5000 cm⁻³ and temperatures of 15–100 K (Snow & McCall 2006). It was one of the very first interstellar molecules detected through its visible bands in the optical and later through its Λ-doublet lines in the ground state in radio frequencies. CH is one of the key probes of the diffuse ISM. Its absorption lines have been detected in the foreground galaxy at z = 0.89 toward the quasar PKS 1830-211 (Muller et al. 2014, 2023). However, its study in the ISM of local galaxies is limited by the absence of millimeter lines and the complete atmospheric absorption of its first submillimeter transition at ∼530 GHz (Jacob et al. 2021).

The lowest spin-rotational (560 μm) transition of CH has six hyperfine components grouped at 532.8 and 536.8 GHz. They were detected in four nearby galaxies, NGC 1068, Arp 220, M 82, and NGC 253, using the Herschel (Rangwala et al. 2011, 2014), who found that the CH lines were a factor of ∼20 brighter than the adjacent HCN and HCO⁺ lines.

While this frequency range was not covered for NCv1.143, the two groups of hyperfine components of the CH(4–3) transition were detected with S/N  >   4 in APM 08279+5255 (Table 5). This is the first clear detection of CH in an individual high-redshift galaxy, although Spilker et al. (2014) presented a marginal 2.5-σ detection of these lines in the stacked spectrum of the dusty star-forming galaxies selected by the South Pole Telescope survey.

As shown in Table 8, the ratio of CH lines to adjacent J = 6–5 lines of HCN and HCO⁺ is more than an order of magnitude lower in APM 08279+5255 than that in the local galaxies observed by Rangwala et al. (2014). This holds true even for NGC 1068 and Arp 220, where the ratios are lower by a factor of approximately 20. The origin of such a substantial difference is challenging to comprehend without a detailed comparison to other high-redshift starbursts and AGN galaxies. Additionally, a higher angular resolution of APM 08279+5255 would be beneficial in better understanding the origin of the CH emission, shedding light on the cause of the CH “deficit.” As the maximum factor of differential lensing reaches only up to 3 (Riechers et al. 2009) in APM 08279+5255, it is unlikely the explanation. Interestingly, however, the average CH/HCO⁺ ratio in the stacked spectrum of SPT-discovered dusty galaxies appears intermediate, with an approximate value of 0.5 ± 0.2 (Spilker et al. 2014), a factor of six times smaller than the local values. In contrast, the CH/HCO⁺ flux ratios reach as high as ∼4–5 in the stacked spectrum of the Herschel lensed sources (Hagimoto et al. 2023). In the absence of alternative explanations, the primary cause of the weakness of the CH lines in APM 08279+5255 might be due to the destruction of CH by the strong quasar UV field, which is akin to the weak [C I] emission in APM 08279+5255 compared with NCv1.143.

Taking the LTE-derived column densities, we derive a relative abundance ratio between CO and CH in APM 08279+5255 of CH/CO ∼ 8 × 10⁻⁵. Unlike the typical Milky Way condition where the diffuse gas is dominating the CH emission, in typical fully shielded molecular regions, the CH/CO abundance ratio is expected to be about 10⁻⁴ (Sheffer et al. 2008), which is consistent with our finding in APM 08279+5255, suggesting the CH lines are unlikely to arise from diffuse gas in APM 08279+5255. Such an abundance ratio is higher than the value found in Arp 220, M 82, and NGC 253 (Rangwala et al. 2014) by at least a factor of 2. However, the abundance ratio of APM 08279+5255 is very close to that of NGC 1068, where the XDR is likely driving the high abundance of CH. This is also consistent with our previous finding of the abundance comparison between AGN and non-AGN sources (Sect. 6.1), as well as the explanation of the abundance ratio of H₃O⁺/HCO⁺ (Sect. 6.4). If CH/CO ratios are indeed a robust XDR tracer as proposed in Rangwala et al. (2014), future observations of the CH lines in starburst galaxies will be revealing, and we will be able to compare the abundance ratios between AGN and non-AGN sources to test the assumption. If the test is successful, the CH/CO ratio can be further used as a diagnostic tool for high-redshift ISM observations to detect XDR imprints.

6.9. CO isotopologs: ¹³CO and C¹⁸O

CO isotopologs are conceivably promising tracers of stellar IMF in high-redshift dust-enshrouded starburst galaxies (e.g., Danielson et al. 2013; Romano et al. 2017; Zhang et al. 2018). In those galaxies, the flux ratio of ¹³CO/C¹⁸O may be assumed to be close to the isotope ratio ¹³C/¹⁸O if both lines of ¹³CO and C¹⁸O are optically thin (however, we caution that in very high-density gas, the ¹³CO lines may have a moderate optical depth; see Appendix D). The isotopes, ¹³C and ¹⁸O, predominantly originate from low- to intermediate-mass stars (less than 8 M_⊙) and massive stars (greater than 8 M_⊙), respectively. Therefore, the flux ratio ¹³CO/C¹⁸O, which represents their abundance ratio, can serve as a dust-insensitive stellar IMF probe. Zhang et al. (2018) demonstrated that low flux ratios ¹³CO/C¹⁸O, approaching unity, found in a sample of starburst galaxies can be interpreted as evidence of a top-heavy stellar IMF. This implication can be further substantiated by high flux ratio values of CO/¹³CO, specifically those exceeding 20.

Zhang et al. (2018) present deep observations of the three isotopologs – CO, ¹³CO and C¹⁸O – toward four high-redshift lensed galaxies, including three starbursts and one quasar. They found flux ratios ¹³CO/C¹⁸O and CO/¹³CO ranging from 1.0 to 1.6 and from 20 to 25, respectively, which can be interpreted as evidence for a top-heavy stellar IMF. In the spectral surveys reported here, we detected all three CO isotopologs in NCv1.143 in the J = 3–2 and 4–3 transitions and, in APM 08279+5255 in the J = 4–3 transition (Table 5). However, the ¹³CO(4–3) line is blended with the CS(9–8) line, which is blueshifted with respect to ¹³CO(4–3) by only 26 km s⁻¹. Therefore, the value given for the flux of the blended line in Table 5 should be corrected for the contribution of the CS(9–8) line. The flux of this line may be estimated by interpolation between the observed values of the CS(8–7) and CS(10–9) lines, assuming a smooth excitation SLED shape. We also acknowledge the substantial uncertainties of the integrated fluxes of the CS lines in NCv1.143 (Table 5), which we take into consideration when computing the line ratios conservatively as upper and lower limits. These ratios listed in Table 9 are consistent for the 3–2 and 4–3 lines, and their values for both sources – ¹³CO/C¹⁸O ≲ 1 and CO/¹³CO ≳ 20 – are comparable to the values reported by Zhang et al. (2018). The similar isotopolog ratios can be interpreted as evidence for a top-heavy stellar IMF in both APM 08279+5255 and NCv1.143, with implications for parameters derived using common assumptions about the IMF, SFR, and stellar mass.

6.10. Rotational lines of the vibrationally excited HCN

From our NOEMA line survey, we have detected for the first time emission lines from the vibrationally excited state v₂ = 1f of HCN (HCN-VIB hereafter) at high redshifts. In the APM 08279+5255 quasar, all three transitions of HCN-VIB, 6–5, 5–4, and 4–3, are detected (Table 5).

These HCN-VIB lines have very high energy levels that exceed 1050 K (Aalto et al. 2015), suggesting that our detection of the bright HCN-VIB lines are unlikely caused only by collisional excitation because of the extreme physical conditions needed (Ziurys & Turner 1986). In local galaxies, observations of the HCN-VIB lines found that the most common process powering these lines involve radiative pumping from the mid-infrared (such as 14 μm), where HCN-VIB is absorbing the mid-infrared in the most obscured regions where N_H reach beyond 10²⁴–10²⁵ cm⁻², then releasing the energy via its emission at submillimeter and millimeter bands through the HCN-VIB lines, forming the so-called compact obscured nuclei (e.g., Sakamoto et al. 2010; Imanishi & Nakanishi 2013; Aalto et al. 2015, 2019; González-Alfonso & Sakamoto 2019; Falstad et al. 2019, 2021). This makes the HCN-VIB lines a powerful tool for probing the presence of extremely obscured nuclei.

Based on the observed integrated flux of the HCN-VIB(4–3) line listed in Table 5, we estimate line luminosities of L_{HCN-VIB(4 − 3)} = (2.5 ± 1.0)×10⁷ L_⊙ for APM 08275+5255. This yields a luminosity ratio of L_{HCN-VIB(4 − 3)}/L_IR = (2.9 ± 1.2)×10⁻⁷. This value is at the high end of the ratios found in compact obscured nuclei and an order of magnitude above the threshold of 10⁻⁸ considered to enter the criteria of a compact obscured nucleus (Falstad et al. 2021). Nevertheless, the HCN-VIB line luminosity could be overestimated if the emitting region is much more compact than the dust continuum and is located at the caustic lines, though we argue that the differential lensing is insignificant (Sect. 8). While, for NCv1.143, a 3σ upper limit of the HCN-VIB lines of 0.13 Jy km s⁻¹ results in a limit to its luminosity of L_{HCN-VIB(4 − 3)} <  1.1 × 10⁷ L_⊙, about a factor of 2–3 below the detected line luminosity in APM 08279+5255. Such a big difference in the HCN-VIB line between APM 08279+5255 and NCv1.143 highlights the fact that the nuclear region of the APM 08279+5255 quasar resides in a much more intense infrared radiation field than NCv1.143. This is consistent with the conclusions we derived from our H₂O excitation model in Sect. 6.3.4. Nevertheless, the upper limit of the ratio L_{HCN-VIB(4 − 3)}/L_IR <  9 × 10⁻⁷ in NCv1.143 is still compatible with a compact obscured nucleus, though it is much less powerful than that of APM 08279+5255. Although a better measurement is needed for NCv1.143, this could be consistent with a buried nuclear starburst or AGN as speculated in Sect. 6.6 should the compact obscured nuclei be unequivocally associated with buried AGN. If such a deeply obscured AGN exists in NCv1.143, its influence on the bulk of the ISM is still limited, consistent with most of the diagnostics from other lines presented in this work.

We also note that the surface luminosity of the HCN-VIB line (converting to J = 3–2 using local line ratios between 4–3 and 3–2) is about 50 L_⊙ pc⁻². Such a value is about five times stronger than the brightest compact obscured nuclei in the local Universe, such as NGC 4418 (surface brightness of HCN-VIB is about 9.1 L_⊙ pc⁻², Falstad et al. 2021), making APM 08279+5255 the brightest HCN-VIB source ever discovered. Yet, unlike most of the compact obscured nuclei, where there have been no X-ray detections (Falstad et al. 2021), APM 08279+5255 is rather unique. As argued in Imanishi et al. (2016), the HCN-VIB line can be boosted by AGN; thus, the unusual brightness of the HCN-VIB in APM 08279+5255 is likely powered by the central AGN.

6.11. Other lines: CN, NO, N₂H⁺, and c-C₃H₂

In APM 08279+5255, we have detected the CN(4–3) emission line, confirming the first detection reported by Guélin et al. (2007) with a significantly better S/N. While in NCv1.143, we have detected both the 3–2 and 4–3 transitions of the CN line. CN can be a dense gas tracer with a few factors lower critical density than HCN (Aalto et al. 2002). In PDRs, CN is primarily formed in the ionization fronts where HCN is dissociated. CN radical can also be formed from CH and N via neutral-neutral reactions (Tielens 2013). On the other hand, in hot-temperature environments, CN can also form HCN directly by interacting with H₂ (Harada et al. 2010). Comparing APM 08279+5255 and NCv1.143, we find that the luminosity ratio between CN and HCN is significantly lower in APM 08279+5255. Comparing directly the abundance ratios (Figs. 6 and 7), we do see the same trend of a smaller ratio of CN/HCN in APM 08279+5255 than in NCv1.143. Such a deficiency of CN with regard to HCN in the quasar APM 08279+5255 indicates that rather than the ionization process, hot-temperature chemistry powered by the AGN might dominate the CN formation/destruction (Harada et al. 2013) in APM 08279+5255.

Nitric oxide NO lines have been detected in APM 08279+5255, the first high-redshift detection of this molecule. NO can form from the interaction of N + OH and NH + O, can be dissociated to N⁺ and NO⁺, and can also form CN by interacting with C; it is thus an important part of the network of nitrogen-bearing species (e.g., Neufeld & Dalgarno 1989). NO can also trace dense gas as its critical densities can reach above 10⁵ cm⁻³. The derived column density of NO in APM 08279+5255 is higher than the values found in local starburst galaxies such as M 82 (Aladro et al. 2011) and NGC 253 (Martín et al. 2021), which can indicate a high abundance of nitrogen and oxygen. Future detailed analysis of the nitrogen chemical-reaction network can help us understand the chemistry.

In APM 08279+5255, we have also detected the 4–3 and 5–4 transitions of the classic dense gas tracer, diazenylium N₂H⁺, confirming the first tentative detection of the N₂H⁺(5–4) line report by Feruglio et al. (2017). Also, we report a tentative detection of the cyclopropenylidene c-C₃H₂ line at 485.7 GHz for the first time at high redshifts, which needs further confirmation.

7. Chemical analysis of the dense molecular gas

The unique richness of observational data for both galaxies presented here makes them ideal systems for measuring the heating and cooling processes as well as determining what powers the energy of their ISM via chemical modeling. Using all available lines simultaneously to perform such modeling, however, is an unprecedented and rather challenging task that deserves its own separate and detailed study.

For the purposes of the present work, we focus on the column densities of four species (HCN, HCO⁺, HNC, and CS), which are dense gas tracers and thus can tell much about the conditions prevailing in the molecular gas of both galaxies. In particular, the contribution from the stellar far-ultraviolet (FUV) photons can be neglected since the FUV radiation will be severely extinguished in the dense regions with volume densities predicted by the LVG models (Table 7). For instance, a density of n_H = 10⁵ cm⁻³ is likely to have an effective visual extinction of A_V ≃ 30 mag (Bisbas et al. 2023), and such visual extinction will extinguish the effect of the FUV field even as strong as 10⁵ χ₀⁵, which is at least one order of magnitude higher than the typical value found in dusty high-redshift starburst galaxies (e.g., Cañameras et al. 2018; Rybak et al. 2019). At such high column densities, the chemistry is no longer driven by FUV radiation but rather by CR, which can be quantified by the CR ionization per H-nuclei ratio (ζ_CR/n_H) that describes the ionization degree. We therefore focus only on ζ_CR/n_H hereafter, which plays an important role in triggering chemical networks at high column densities (e.g., Bialy & Sternberg 2015).

For the purposes of this chemical modeling, we use the 3D-PDR⁶ code, which is an astrochemical code treating photodissociation regions (Bisbas et al. 2012). Using various heating and cooling functions and the LVG approximation (Sobolev 1960; Poelman & Spaans 2005), the code performs iterations over thermal balance and terminates once the total heating is approximately equal to the total cooling. In this work, the full UMIST2012 (McElroy et al. 2013) chemical network is used, consisting of 215 species and approximately 3000 reactions. The initial abundances (normalized to total hydrogen) adopted here are [Mg] = 3.981×10⁻⁵, [C⁺] = 10⁻⁴, [S⁺] = 3.236×10⁻⁵, [N] = 6.76×10⁻⁵, [O] = 3×10⁻⁴ [He] = 0.1. A metallicity of Z = 1 Z_⊙ (e.g., Dunne et al. 2022) and a microturbulent velocity of v_turb = 1 km s⁻¹ are further assumed.

In addition, we adopted the following expression of optical depths, which best describes macroturbulent molecular clouds (Papadopoulos & Seaquist 1999):

$$\begin{array}{c}\hfill {\tau }_{\mathit{ij}}=\frac{A_{\mathit{ij}}{c}^3}{8\pi {\nu }_{\mathit{ij}}}{\left(\frac{\mathrm{d}V}{\mathrm{d}r}\right)}_{\mathrm{VIR}}^{-1}\left[\frac{n_j{g}_i-n_i{g}_j}{g_j}\right],\end{array}$$(4)

where

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}V}{\mathrm{d}r}\right)}_{\mathrm{VIR}}=0.65\sqrt{\alpha }{\left[\frac{n_{\mathrm{H}}}{{10}^3{\mathrm{cm}}^{-3}}\right]}^{1/2}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{pc}}^{-1}.\end{array}$$(5)

In the above equations, A_ij is the Einstein coefficient of levels i and j, ν_ij is the corresponding line frequency, n_i and n_j the level populations, and g_i, g_j their statistical weights. The parameter α depends on the density profile (Papadopoulos & Seaquist 1999). Here, we have considered $0.65\sqrt{\alpha }=1$ (see also Bisbas et al. 2015). Equation (4) is used based on the assumption that in macroturbulent clouds, the optical depth is considered as a local property rather than building along the line-of-sight as treated for local clouds.

Fig. 14. Chemical modeling using 3D-DPR to identify the best range of the CR ionization parameter ζCR/nH (x-axis) based on the observed column density ratios. The column densities of HCN (blue), HCO+ (orange), and HNC (green) are normalized with CS. The y-axis shows the ratio of the modeled column density ratio to the observed one. Unity marks agreement between the models and observations. The upper panel shows the results for APM 08279+5255, and the lower panel shows those for NCv1.143. The insets are zoom-in regions of the gray rectangular area in each panel. The thickness of each line represents the uncertainty (light color is for 50% uncertainty, and dark color is for 10% uncertainty).

Absolute column densities strongly rely on the assumptions made on the size of the emitting source, which also introduces uncertainty on the opacity of the emission. In order to minimize this uncertainty and the matching to models, we base our analysis on the normalized column density. In this case, we normalize it to the column density of CS as a reference, which will be less affected by optical thickness than, for example, CO. Here, we assume that the spatial distribution among HCN, HCO⁺, HNC, and CS is similar; however, due to the limit of our current data, this is not certain. Figure 14 shows the results of our chemical modeling in which we plot for both sources the ratio ℛ = R_N,sim/R_N,obs versus ζ_CR/n_H, where ℛ is the ratio between the modeled column densities ratios and the observed ones from Table 7. We adopt two uncertainties for this comparison, a 10% (dark color) and a 50% (light color) of the mean reported in this table, for the purpose of demonstrating how the measurement uncertainties impact our analysis. However, we acknowledge that our estimations of CR ionization rates have large uncertainties.

Tables 10 present the modeled results for NCv1.143 and APM 08279+5255, respectively. It appears that for both cases, no single solution for the column densities considered is found. Using the n_H number density of CS from Table 7 it can be seen that the estimated range of CR ionization rates of ζ_CR ∼ (4 − 40)×10⁻¹⁴ s⁻¹ for APM 08279+5255 (although the high ζ_CR/n_H solution of the HCN/CS ratio suggests a higher value of ζ_CR ∼ 10⁻¹² s⁻¹), and of ζ_CR ≃ (8 − 12)×10⁻¹³ s⁻¹ for the NCv1.143 galaxy. In both galaxies, the HNC/CS ratio shows systematically lower values of ζ_CR, while the results from the HCN/CS and HCO⁺/CS ratios generally agree with each other.

Similar high ζ_CR values, as those found here, have been previously reported in systems from the Galactic Center to high-redshift galaxies – Bremsstrahlung radiation (Yusef-Zadeh et al. 2013) and ${}_3^{+}$ observations (Goto et al. 2014; Le Petit et al. 2016) in the Galactic Center suggest a CR ionization rate ζ_CR ≃ (1 − 100)×10⁻¹⁵ s⁻¹, which is 2–4 orders of magnitude higher than the Milky Way average. The nearby and well-studied NGC 253 starburst is suggested to have a ζ_CR ≳ 10⁻¹⁴ s⁻¹ derived from the observed HCO⁺/HOC⁺ emission ratio (Harada et al. 2021). Similarly, Holdship et al. (2021) report an ionization rate ζ_CR ∼ 10⁻¹¹ s⁻¹ derived from chemical models of C₂H. However, in a subsequent paper and using H₃O⁺ and SO observations, Holdship et al. (2022) suggest a lower value of ζ_CR ≃ (1 − 80)×10⁻¹⁴ s⁻¹ for NGC 253. Very high ζ_CR values have also been derived for other star-forming galaxies. OH⁺, H₂O⁺ and H₃O⁺ observations in Arp 220 find a ζ_CR/n_H ≃ (1 − 2)×10⁻¹⁷ cm³ s⁻¹, which leads to a ζ_CR ≳ 10⁻¹³ s⁻¹ for an estimated n_H ∼ 10⁴ cm⁻³ (González-Alfonso et al. 2013). Similarly, the aforementioned lines observed for Mrk 231 also suggest a ζ_CR/n_H ≃ (1 − 2)×10⁻¹⁶ cm³ s⁻¹ for the outflowing component implying a ζ_CR ∼ (0.5 − 2)×10⁻¹² s⁻¹ (González-Alfonso et al. 2018). Imanishi et al. (2018) studied OH⁺ and H₂O⁺ lines of the extended halos of z ∼ 2.3 galaxies and derived a ζ_CR ∼ (1 − 100)×10⁻¹⁶ s⁻¹, where the ionization rates in the compact star-forming regions are expected to be orders of magnitude higher (Indriolo et al. 2018). Muller et al. (2016) also found higher than Milky Way CR ionization rate in the z = 0.89 molecular absorbers toward PKS1830-211. All these results suggest that the dense molecular gas in these high-redshift star-forming galaxies is subjected to high CR ionization rates. As pointed out by Krumholz et al. (2023), assuming star formation is the dominant source of CR, ζ_CR scales with the gas depletion time. Therefore, such high ζ_CR values of ∼10⁻¹⁴ − 10⁻¹² s⁻¹ indicates a very short gas depletion timescale of possible only a few to a few tens of megayears, taking the prescription from Krumholz et al. (2023), consistent with their intense starburst nature.

Fig. 15. Abundances of HCN, HCO+, HNC and CS versus nH2, for two different CR ionization rates (solid lines for ζCR = 10−13 s−1 and solid-dotted lines for 10−12 s−1). Our model shows that with such a high ionization rate ζCR = 10−13–10−12 s−1, the abundances of HCO+ is expected to be lower than HCN, HNC, and CS in the dense gas regions of nH2 > 105 cm−3. However, with increasing CR ionization rates, at the same H2 densities, the abundance of HCO+ is rising for the dense gas while the abundance of CS is decreasing.

Taking the derived CR ionization rate values in Table 10 for both galaxies, in Fig. 15, we show the variation of the abundances in HCN, HCO⁺, HNC, and CS with gas density. For APM 08279+5255, where the gas densities reach above n_H2 ∼ 10⁵ cm⁻³ and ζ_CR ≃ 10⁻¹³ s⁻¹, we found that our chemical model predicts a low HCO⁺ abundance compared to other species, lower by about two orders of magnitude. This is strikingly consistent with the LVG-derived column densities, as shown in Fig. 12, that the abundance of HCO⁺ is about two orders of magnitude lower than that of HCN, HNC, and CS. For NCv1.143, if we assume a ζ_CR ≃ 10⁻¹² s⁻¹, we see a significant increase in the relative abundance of the HCO⁺ (a detailed discussion of the chemistry of HCO⁺ is given in Appendix E), and a small decrease in the CS abundance while HCN and HNC remain similar values of abundances. For NCv1.143, both a higher gas density and likely a higher ζ_CR lead to similar relative abundances among the dense gas tracers, consistent with Fig. 12. It is also worth noticing that, given such a CR ionization rate, the relative abundance of HCO⁺ drops toward the denser conditions from n_H2 ∼ 10⁵ cm⁻³, while CS abundance increases with increasing gas densities. Around the gas densities of n_H2 ∼ 10^5 − 6 cm⁻³, the abundance ratio of HCO⁺/CS is most sensitive to the change of the CR ionization among the four species, with the HCO⁺ abundance rising and that of CS decreasing as ζ_CR increases.

We note that the discussions above assume that volumetric heating at high column densities arises from the presence of CRs only, neglecting other heating sources that may also contribute to gas heating. For instance, strong X-ray sources can equally contribute to the CR heating and also affect the chemistry in a similar way (e.g., Meijerink et al. 2007; Gallerani et al. 2014), though they are alternated much more quickly than CRs and ∝R⁻². Should such heating sources be accounted for, they will reduce the values of CR ionization rates we find. We therefore note that the ζ_CR values presented above shall be considered as upper limits.

Nevertheless, regardless of whether the heating is caused by CRs only or with contributions from X-rays and shocks, the initial conditions of star formation will be altered under such ISM conditions, causing a boosted Jean mass and a top-heavy stellar IMF (Papadopoulos 2010; Papadopoulos et al. 2011). As predicted by hydrodynamical simulations, such high-temperature gas of ∼100 K and H₂ density of ≳10⁵ will lead to a stellar IMF peaking at ≳15 M_⊙ (Klessen et al. 2007). Moreover, such a top-heavy stellar IMF is consistent with the observed CO isotopolog ratios in both galaxies, which have been interpreted as a result of enrichment by massive stars (Sect. 6.9). A similar scenario is presented in the lensed starburst galaxy SMM J2135-0102, where strong gas heating (likely by CR) leads to an enhanced C¹⁸O abundance (Danielson et al. 2013).

8. Impact of the lensing on the analysis

Given that our observations are based on globally integrated values using unresolved interferometric data, we do not have information on how each line emission is distributed spatially. In both sources, the background galaxy is strongly lensed, while the foreground galaxies are not visible at longer wavelengths and thus do not contaminate our observations. The lensing can bring potential differential lensing effects, where the magnification factor can be different for emissions coming from different parts of the lensed galaxies (e.g., Riechers et al. 2011b; Serjeant 2012; Hezaveh et al. 2012; Fu et al. 2012; Yang et al. 2019; Dye et al. 2022). This can bring more uncertainties to our analysis because we assume minimum differential lensing across all the molecular lines, which is, of course, a first-order approximation, and yet the only assumption we can make until high-angular-resolution imaging of individual lines is obtained.

For NCv1.143, to better understand how lensing magnification impacts our results, we built the lens model based on 0.3″ resolution data of the CO and H₂O, as described in Appendix C. We find a median lensing magnification of about 12 for the dust continuum and about 8–10 for the CO(10–9), H₂O(2₁₁ − 2₀₂) and H₂O(3₂₁ − 3₁₂) emission. Because the line profiles are not strongly distorted by the magnification, and the resolved line profiles of mid/high-J CO and H₂O(2₁₁ − 2₀₂) and H₂O(3₂₁ − 3₁₂) are similar, we can expect that the differential lensing is negligible among the dense gas tracers, at least. This is further supported by the fact that almost all the line profiles from our line survey are similar, except for the 380 GHz H₂O maser line. However, we caution that some of the diffuse gas tracers in NCv1.143 can have a significantly smaller magnification factor, likely leading to an underestimation of the intrinsic fluxes.

For APM 08279+5255, we do not expect significant differential lensing. As noted by Riechers et al. (2009), even under the most extreme scenarios, the magnification fluctuates mildly, from ≲5 with small variations for a source radius up to 3 kpc, down to 3 if the source extends further past 10 kpc. Given that such an extensive expansion is highly improbable for the molecular emissions observed, substantial differential lensing in APM 08279+5255 is unlikely.

9. Conclusions

Using NOEMA, we conducted deep 3 mm spectral line surveys toward a dusty quasar, APM 08279+5255 at z = 3.911, and a dusty star-forming galaxy, NCv1.143 at z = 3.565, with total a bandwidth coverage of about 200 GHz in the rest frame, from around 330 to 550 GHz. The line survey data allow us to study the high-redshift ISM with an unparalleled level of detail. Our main findings are:

• Utilizing MADCUBA (Martín et al. 2019a) plus a matched-filtering method, we detect 38 emission lines in APM 08279+5255 and 25 emission lines in NCv1.143, from [C I], CO, ¹³CO, C¹⁸O, CN, CCH, HCN, HCO⁺, HNC, CS, C³⁴S, H₂O, H₃O⁺, NO, N₂H⁺, CH, and c-C₃H₂, as well as the rotational levels of the vibrationally excited HCN (HCN-VIB). The CH, CCH, c-C₃H₂, N₂H⁺, HCN-VIB, and H₃O⁺ lines are the first-ever individual high-redshift detections of any transitions of the species. The luminosities of the lines between the two sources behave differently: the SLEDs of most of the tracers are more elevated in APM 08279+5255, indicating more extreme excitation conditions. The line luminosities are higher in APM 08279+5255, except for the [C I](1–0) line. The most significant difference is seen in 448 GHz H₂O line, where APM 08279+5255 is about nine times brighter than NCv1.143 in flux.

• Under the LTE assumption, we derived the relative abundances of the molecules for both sources. By comparing them with some prototype local galaxies, we find that NCv1.143 is more similar to that of the CMZ of the starburst-dominated galaxy, NGC 253, while APM 08279+5255 behaves similarly to some local galaxies where the ISM is expected to be dominated by AGN, such as the central ∼2 kpc regions of NGC 1068 and NGC 4418. We have also found that the abundances of dense gas tracers are more enhanced in the AGN-dominated sources. Additionally, we find a lower CN/HCN ratio in APM 08279+5255 than in NCv1.143, suggesting that the AGN-powered hot-temperature chemistry that converts CN into HCN, is important in APM 08279+5255.

• We performed non-LTE LVG modeling of HCN, HCO⁺, HNC, CS, and C³⁴S, which are trace much denser gas than CO: in the range n_H2 ∼ 10⁵–10⁷ cm⁻³. The kinetic temperatures derived from all the species are consistent with T_kin ∼ 50–130 K. Interestingly, the LVG-derived column densities are consistent with the LTE-derived values, and we highlight a generally low abundance of HCO⁺ compared with other dense gas tracers.

• In APM 08279+5255, we detect the 380 GHz H₂O maser lines, which contain high-velocity components that are likely linked to the maser spots located on the fast-rotating accretion disk around the central SMBH. This 380 GHz feature is the only identified maser line in APM 08279+5255 from the whole rest-frame frequency of ∼330–550 GHz. The flux of this line is at least three times stronger than the undetected 22 GHz, making it a very promising submillimeter H₂O megamaser candidate in high-redshift galaxies. Simply scaling from local galaxies, we crudely derived the SMBH mass, which is about 10¹⁰ M_⊙, consistent with previous literature values.

• After including the two water emission lines detected here for the first time, we built a H₂O excitation model using multiple transitions of H₂O and the dust continuum. We find the H₂O excitation conditions in both sources are best explained by two components, the more extended of which are about 1.3 and 0.9 kpc in APM 08279+5255 and NCv1.143, respectively, and both have a dust temperature of ∼50 K. The compact component of APM 08279+5255 (R ∼ 320 pc) has a high dust temperature of about 211 K that dominates the total infrared luminosity output, and that of NCv1.143 (R ∼ 170 pc) has a dust temperature of about 98 K with a minor contribution to the infrared luminosity. Such a difference in their compact component shows that APM 08279+5255 has a highly prominent compact core that exhibits extreme conditions and dominates the total power output at the far infrared, indicating that the AGN in APM 08279+5255 is in the “turn-on” phase, in contrast with NCv1.143. This showcases the power of H₂O lines as a probe of ISM structures in dusty galaxies.

• While the line flux ratios among HCN, HCO⁺, and CS are inclusive of the imprint of the XDRs versus the PDRs, the abundance ratios from the LVG analysis suggest that the AGN plays an important role in regulating the abundance of these dense gas species in APM 08279+5255. This is further supported by the abundance ratios between H₃O⁺ (detected in both APM 08279+5255 and NCv1.143, the first high-redshift detections) and HCO⁺, where we find that XDRs play an important role.

• We have also detected the isotopologs of CO (¹³CO and C¹⁸O) in both sources. We find similar flux ratios of ¹³CO/C¹⁸O and ¹²CO/¹³CO, as in other starburst galaxies. These ratios are interpreted as evidence for a top-heavy stellar IMF as proposed by Zhang et al. (2018), which can significantly impact estimations of the SFRs.

• In APM 08279+5255, we have detected the HCN-VIB lines (rotational transition from J_up = 4 to 6) at high redshifts for the first time. These lines probe deep into the compact obscured nucleus of the quasar. In NCv.1.143, the upper limit of the HCN-VIB line does not rule out the possibility of a deeply buried AGN. Our calculation of the HCN-VIB surface brightness makes APM 08279+5255 the brightest HCN-VIB emitter discovered to date.

• Using a three-dimensional PDR model, we explored the properties of the ISM and the radiation fields by taking the chemistry into account. Using the LVG-derived abundance ratios, we find that both our sources are extremely dust-obscured and that FUV radiation plays a limited role because of the high A_V. However, we do find high ionization rates in both sources at a level of 10⁻¹³ to 10⁻¹² s⁻¹, which is likely caused by high densities of CRs, although we do not rule out the possibility of contributions from X-rays. Such an extreme condition is likely altering the initial conditions of star formation, leading to a top-heavy stellar IMF, consistent with the enhanced C¹⁸O/¹³CO ratio observed. Additionally, our chemical model predicts a low abundance of HCO⁺ compared to HCN, HNC, and CS in such high-density and high-CR-ionization conditions. This provides a possible explanation of the relatively low abundance of HCO⁺ derived from our LVG modeling.

• We built detailed models taking chemistry and radiative transfer into account. However, we caution that our line fluxes are derived from spatially unresolved integrated values. In reality, the spatial distribution of each molecular emission can vary significantly across a single galaxy, depending on the ISM conditions in different regions. This has been observed in the pure-starburst prototype galaxy NGC 253 (Walter et al. 2017; Martín et al. 2021), where the CMZ has a very distinct chemical process and thus different spatial distributions of molecular emission than its galactic disk. The variation is even more pronounced when a central AGN is present, as shown in the case study of NGC 1068 (García-Burillo et al. 2014). In our line survey data, the intrinsic variations found in NGC 253 and NGC 1068 add more uncertainties and lead to difficulties when interpreting the data, although one can group several tracers according to the typical conditions they trace in common and assume that their emission distributions are coincident. Nevertheless, high-resolution imaging of the individual line emission is needed for a more comprehensive picture of the ISM in APM 08279+5255 and NCv1.143.

With the drastic improvement of the bandwidths and sensitivities of (sub)millimeter interferometers, at high redshifts, we can now detect not only CO, the traditional low-density gas tracer, but also the emission of a score of high-excitation and more complex molecules, including some rare isotopologs. With the richness of lines detected from these molecules, we are able to understand the physical conditions, structures, and chemical processes of the ISM at high redshifts in details never before achieved. Future high-angular-resolution line survey data will help us further break the degeneracies in the models of gas excitation and ISM chemistry.

Acknowledgments

We thank the anonymous referee for many helpful comments and suggestions. This work is based on observations carried out under project numbers W15EP, S18DC, and W18EB with the IRAM NOEMA Interferometer. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). The authors are grateful to the IRAM staff for their support. Chentao Yang and Susanne Aalto acknowledge support from ERC Advanced Grant 789410. Chentao Yang gratefully thanks Serena Viti for the enlightening discussions during several conferences. Chentao Yang thanks Melanie Krips for her involvement in the observation proposals. Chentao Yang is in debt of the insightful discussion with Santiago Del Palacio about X-ray corona emission, as well as the wonderful discussions with Jim Braatz, Dom Pesce, Violette Impellizzeri, Christian Henkel, Elizabeth Humphreys, and Chengyu Kuo on H₂O masers. Chentao Yang also thanks Masatoshi Imanishi and Matus Rybak for the nice discussions. EG-A thanks the Spanish MICINN for support under projects PID2019-105552RB-C41 and PID2022-137779OB-C41. Kirsten Knudsen acknowledges support from the Swedish Research Council and the Knut and Alice Wallenberg Foundation. Rob J. Ivison acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. We acknowledge the usage of the following Python packages: Astropy (Astropy Collaboration 2018), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), PyRadex (https://github.com/keflavich/pyradex) , Emcee (Foreman-Mackey et al. 2013), and Corner (Foreman-Mackey 2016). And the following Julia package: Cosmology.jl (Tomasi & Giordano 2019), and Measurements.jl (Giordano 2016).

References

Appendix A: Zoomed-in view of the spectra and the fittings

Fig. A.1. Same as Fig 3, but zoomed in to show details of the line profiles.

Fig. A.2. Fitting results of the spectra of APM 08279+5255, with the x-axis as the rest-frame frequency and the y-axis as the observed flux density. The top panel shows the observed spectrum in blue with the uncertainties in the green line. The red solid lines indicate the fit to the entire spectrum, and the dashed red lines are the position of the identified lines. The lower panel shows the residual.

Fig. A.3. Same as Fig. A.2 but for NCv1.143.

Appendix B: Posterior distribution of the LVG parameters

This section shows the posterior probability distributions of molecular gas density n_H2, gas temperature T_kin and column density per velocity N_mol/dv of HCN, HCO⁺, HNC, CS and C³⁴S in APM 08279+5255 and NCv1.143.

Fig. B.1. Posterior probability distributions of molecular gas density (nH2), gas temperature (Tkin) and column density per velocity (Nmol/dv) of each species in APM 08279+5255 and NCv1.143, with the maximum posterior possibility point in the parameter space shown in orange lines and points. The contours are in steps of 0.5 σ starting from the center.

Fig. B.2. Same as Fig. B.1 but for the CS and C34S lines.

Appendix C: Lens model

Our lens model for NCv1.143 was based on high-spatial-resolution data gathered from NOEMA under project W15EP, which has improved sensitivity and spatial resolution compared with the Submillimeter Array (SMA) observations from Bussmann et al. (2013). Additionally, the NOEMA data also covers the multiple molecular line emissions, including two H₂O, one H₂O⁺, and the CO(10–9) lines.

NOEMA observations were conducted under excellent atmospheric conditions (seeing of 0.3″–1.0″and PWV≤1 mm) from February to March 2016, using A-configuration (baselines are from 93 to 760 m). The nominal angular resolution reached is approximately 0.5″×0.3″at 1.2 mm and 0.7″×0.4″ at 1.8 mm. Seven antennas were used with a total on-sources time of 5.7 hours for the entire series of observations. Standard NOEMA calibrators, including 3C279, 3C273, MWC349, and 0923+392, were used for phase and bandpass calibration. Calibration, imaging, CLEANing, and spectra extraction were executed within the GILDAS packages CLIC and MAPPING.

In Fig. C.1 we present the 1.8 mm and 1.2 mm dust continuum of NCv1.143. These images reveal dust continuum emissions at approximately 390 μm and 260 μm in the rest frame. We also show the line emission of the CO(10–9), H₂O(2₁₁ − 2₀₂), and H₂O(3₂₁ − 3₁₂) lines and the sum of the emissions from the four H₂O⁺ lines. The images show arc-like emission components, which are typically found in galaxy-galaxy lensing systems.

Utilizing a simple lens model such as a singular isothermal ellipsoid (SIE), as used in NCv1.143 and outlined in Kormann et al. (1994), we can reproduce lensed signatures and determine deflector properties based on mass profile assumptions. Following the procedures outlined by Bussmann et al. (2013) and Gavazzi et al. (2011), we used the SIE model and the GRAVLENS software to analyze the source plane emission and the NOEMA CLEANed images. These techniques yield the parameters of both the deflector’s SIE profiles and the lensed source Sérsic profiles.

Our images of NCv1.143, with a resolution of 0.3′′ − 0.4′′, surpass the resolution of the SMA images from Bussmann et al. (2013) and represent the highest spatial resolution data of NCv1.143 to date. In addition to dust continuum images, we obtained data cubes, including line emission images at various velocity channels. We assumed that the deflector comprises a single SIE mass distribution, allowing the center to vary within 0.5′′ of the optical center of the foreground galaxies. The emitted light was approximated by an elliptical exponential light profile, as supported by Bussmann et al. (2013). The modeling was accomplished with the sl_fit code, using an MCMC method to explore the parameter space.

Using the integrated continuum dust emission for mass distribution and dust emission constraining (simultaneously fitting the 1.2 mm and 2 mm continuum images), the best-fit model was then applied to retrieve intrinsic gas properties. We believe this approach adequately captures the key information of the data. The results of our lensing model are detailed in Table . The positions of line emissions are reported for the channel closest to the frequency center, and magnification is integrated over the entire line.

The model aligns well with the observed continuum images (Fig. C.1), allowing us to accurately reconstruct the image of the dust continuum emission in the source plane. The average magnification factor of the dust continuum, μ_2, mm, was found to be 12.2±1.2, and its radius of extension, R_eff, was determined to be (0.8±0.1) kpc. These results are consistent with, but more precise than, the μ₈₈₀ = 11.3±1.7 value reported by Bussmann et al. (2013).

Fig. C.1. NOEMA images of NCv1.143. Color images are the ∼1.8 mm and ∼1.2 mm dust continuum. Yellow, red, and blue contours are the CO(10–9), H2O, and H2O+ lines. The synthesized beams are 0.76″×0.41″and 0.54″×0.30″for the 1.8 and 1.2 mm bands, respectively. The contour levels are: in steps of 1 σ started from ± 3 σ for the H2O(211 − 202) and the H2O+ lines; and in steps of 3 σ started from ± 3 σ for the H2O(321 − 312) and the CO(10–9) lines.

Fig. C.2. Posterior possibility distributions of the Reff and μ from the MCMC sampling. The contours are from 1 σ to 3 σ from inside out. The solid blue line, solid red line, dashed purple line, solid yellow line, and dashed black line indicate the results for H2O(211 − 202), H2O(321 − 312), H2O+, CO(10–9), and the dust continuum, respectively.

Fig. C.3. Lens model of the 254 GHz (top row) and 164 GHz (bottom row) dust continuum emission in NCv1.143. From left to right: (1) Observed image of the 2 mm continuum. (2) Model images with the model parameters of Table . (3) Residuals from the difference between observed and model images. (4) Reconstructed image of the dust continuum emission in the source plane. The critical curves and caustics are shown in red lines in the first three columns and the fourth columns, respectively. The position of the deflector is indicated by red points.

Fig. C.4. From top to bottom, each row: Lens model of the CO(10–9), H2O(321 − 312) and H2O(211 − 202) lines in NCv1.143. From left to right: (1) Observed CLEANed images of the emission line. (2) Model images with the model parameters of Table . (3) Residuals from the difference between observed and model images. (4) Reconstructed images of the emission line in the source plane. The critical curves and caustics are shown in cyan lines in the first three columns and the fourth column, respectively. The position of the deflector is indicated by the cyan points. One should note that the position of the background source in the source plane is not well-constrained. Thus, the differences seen in the positions of the three lines are not significant, considering the uncertainties.

Appendix D: Optical depths

Fig. D.1. 3D-PDR calculations using a density distribution representing the AV, eff − nH relationship (Bisbas et al. 2019). The figure shows the optical depth of 13CO(4–3) and C18O(4–3) lines as a function of the H2 volume density. We find that both lines remain optically thin for nH2 ≲ 106 cm−3), while the 13CO(4–3) starts to have moderate optical depth with increasing H2 densities.

We performed a PDR calculation to explore how the optical depth of ¹³CO(4–3) and C¹⁸O(4–3) builds along with the cloud depth. We constructed a density distribution resulting from the A_V,  eff − n_H relation presented in Bisbas et al. (2019) and Bisbas et al. (2023). This empirical relation connects the local visual extinction with the most probable local number density. The calculation was performed using an FUV intensity of χ/χ₀ = 10 (normalized to the spectral shape of Draine 1978), a CR ionization rate of ζ_CR = 10⁻¹³ s⁻¹ (as derived from Sect. 7), and a metallicity of 1 Z_⊙. The microturbulent velocity controlling the turbulent heating was taken to be v_turb = 10 km s⁻¹, and the velocity dispersion was taken to be v_disp = 180 km s⁻¹ to be in accordance to the observed linewidth. Furthermore, we assumed a ¹²C/¹³C isotope ratio of 30 and ¹²CO/C¹⁸O of 200 (e.g., Martín et al. 2019b; Tang et al. 2019).

Figure D.1 shows the optical depths of the ¹³CO(4–3) and C¹⁸O(4–3) lines as a function of the H₂ density, n_H2. For values of n_H2 ≲ 10⁶ cm⁻³, both lines remain optically thin. However, when n_H2 reach beyond 10⁶ cm⁻³, the ¹³CO(4–3) line starts to have a moderate optical depth while C¹⁸O(4–3) remains optical thin until the density getting close to 10⁷ cm⁻³. Both lines start to get optically thick if n_H2 ≳ 10⁷ cm⁻³. In the conditions of the n_H2 between 10⁶ and 10⁷ cm⁻³, the optical depth of the ¹³CO(4–3) will reach about 10. Such an optical depth will lead to a low flux ratio of ¹³CO/C¹⁸O close to unity and a very low flux ratio of ¹²CO/¹³CO, as pointed out by Zhang et al. (2018). However, the observed ¹²CO/¹³CO ratios in both our sources (Table 9) are about 20–30, inconsistent with a optically thick scenario of the ¹³CO.

Appendix E: Detailed chemical modeling of HCO⁺

Fig. E.1. Chemical models for HCO+. Left panel: HCO+ abundance versus the local total H-nucleus number density, nH, for different CR ionization rates. In all cases, the abundance of HCO+ peaks for densities 8 × 103 ≲ nH ≲ 3 × 105 cm−3. Right panel: Temperature profile, with the iso-(ζ/n) lines overplotted in orange, for different CR ionization rates (ζCR = 10−16 − 10−12 s−1, with the same color code as the left panel).

The left panel of Fig. E.1 displays the HCO⁺ abundance versus the local number density, n_H, for five different CR ionization rates (ζ_CR = 10⁻¹⁶ − 10⁻¹² s⁻¹). The FUV radiation is neglected in these models, as mentioned in the main text, and the macroturbulent optical depth expression (Eq. 4) is adopted. Overall, the abundance of HCO⁺ peaks at densities in the range of n_H ≃ 10⁴ to a few 10⁵ cm⁻³, with the effect to be better seen when ζ_CR ≳ 10⁻¹³ s⁻¹. In the low end of the above ζ_CR values, HCO⁺ is formed mainly via the reactions

$$\begin{array}{cc}& {\mathrm{H}}_2+{\mathrm{HOC}}^{+}\to {\mathrm{HCO}}^{+}+{\mathrm{H}}_2\hfill \\ \hfill & {\mathrm{H}}_2+{\mathrm{CO}}^{+}\to {\mathrm{HCO}}^{+}+\mathrm{H}\hfill \end{array}$$(E.1)

E.1 and E.2 together form approximately 80% of the total HCO⁺ abundance, with E.2 becoming more effective for densities n_H ≲ 10² cm⁻³. For n_H ≳ 10³ cm⁻³, the reaction

$$\begin{array}{c}\hfill {\mathrm{H}}_3^{+}+\mathrm{CO}\to {\mathrm{HCO}}^{+}+{\mathrm{H}}_2\end{array}$$(E.2)

dominates the formation route of HCO⁺. For very high ζ_CR values, HCO⁺ abundance builds rapidly when n_H≳ a few 10³ cm⁻³, below which is in practice negligible (Fig. E.1). Contrary to the aforementioned case, E.2 remains, here, the main reaction that is able to form more than 50% of the total HCO⁺ abundance. In all cases, HCO⁺ is destroyed with an efficiency of almost 100% due to its reaction with free electrons, leading to the formation of CO and H.

The right panel of Fig. E.1 shows the temperature profile of HCO⁺. The temperature profile varies as a function of the CR ionization rate since it inserts heating. We performed various runs with different CR ionization rates, and we have overplotted iso-(ζ/n) lines, the range of which corresponds to the solutions we found for the two galaxies. Assuming the typical density for the dense gas, we find that their gas at high column densities is relatively warm (∼ 50–100 K), consistent with the LVG-derived gas temperatures.

For low densities, an interesting behavior occurs in which the gas temperature decreases as the CR ionization ray increases. This is counter-intuitive, given that CRs heat up the gas. The decrease in gas temperature occurs because the production of C⁺ abundance increases, making the emission line of [C II] increase. The latter is a strong coolant, which in turn decreases the overall gas temperature at low densities (Bisbas et al. 2017).

All Tables

All Figures

Fig. 1. Tuning setups of the NOEMA observations. Top panel: Atmospheric transmission curve at the zenith of the NOEMA site for the frequency range covered by our line survey. The red and blue lines show the transmission curves for PWV = 1 mm and PWV = 5 mm, respectively. Middle and bottom panels: Tuning setups, their sky (observed) frequency coverage, and the RMS noise reached in each spectral window for APM 08279+5255 and NCv1.143. For display purposes, we smoothed the spectral resolution to 50 MHz, resulting in an RMS that is five times smaller than the noise level reported in Table 1 for the 2 MHz bins. As expected, the edges of the sidebands and the conjunction channels of the two spectral windows within each sideband exhibit large RMS values.

In the text

Fig. 2. NOEMA CLEANed images of APM 08279+5255 and NCv1.143, made from the visibilities obtained by averaging all the USB channels from the highest-frequency tunings (w18eb003 and s18dc003; Table 1). These tunings were chosen since they are expected to offer the highest spatial resolution among all the spectral windows in our tunings. The contour levels begin at 5σ and increase in increments of 10σ. The unit of the color bar and contours is arbitrary, with values representing fluxes integrated over the entire USB of the tuning. The primary objective here is to illustrate that the sources are unresolved rather than provide precise flux measurements. The offset for APM 08279+5255 is explained in Sect. 3.

In the text

Fig. 3. NOEMA 3 mm band spectral surveys of APM 08279+5255 and NCv1.143, shown as the blue and violet spectra, respectively. For visualization purposes, the spectrum of NCv1.143 has been shifted up by 2 mJy. Both spectra are binned to 50 MHz (about 150 km s−1) to highlight the line detection. Error bars are indicated by thin lines (the same color as the spectra) overlaid on the spectra. Dashed black lines identify the lines detected in both sources, while dashed orange lines highlight the lines that are only detected in one of the two sources. The line names are labeled above the dashed lines. The peak values indicated for CO(3 − 2), CO(4 − 3), and [C I](1-0) emission lines in NCv1.143 are the true peak values before the 2 mJy shift, and they are outside the box of the figure. A zoomed-in view is provided in Fig. A.1.

In the text

Fig. 4. Fitting of the continuum of APM 08279+5255 and NCv1.143. Fluxes and errors are in black and blue, respectively. The solid red lines represent the total flux from the model. The dotted and dashed orange lines are the free-free and thermal dust components from the model, respectively. It is clear from the fitting that APM 08279+5255 has a flatter spectrum, which can be explained by contributions from free-free and/or corona emission from the SMBH, while the continuum of NCv1.143 is dominated by thermal dust.

In the text

Fig. 5. Line luminosities of APM 08279+5255 and NCv1.143 (based on the line fluxes listed in Table 5, where we also include several detections from the literature). Error bars are propagated from the measured flux errors. The upper panel shows the line luminosity in K km s−1 pc2, $L_{\mathrm{line}}^{\prime }$, normalized by the value of C18O(4–3). The lower panel shows the line luminosity in L⊙, Lline, normalized by the C18O(4–3) line. The lowest levels of the transitions of each molecule detected are labeled at the top of each panel, except for H2O and H3O+, where the rest-frame frequencies in GHz are given. Each transition is identified in the upper part of the figure. Compared with NCv1.143, the excitation condition is more extreme in APM 08279+5255 as seen from the SLEDs from several molecules.

In the text

Fig. 6. Column density ratios to C18O in APM 08279+5255 and NCv1.143 derived by MADCUBA (Table 4), which can be treated as the relative abundances. The plot is ordered, from left to right, according to the relative abundance in APM 08279+5255 from high to low. Some local sources from the literature are also shown: NGC 253 (CMZ; Martín et al. 2021; Rangwala et al. 2014), NGC 1068 (inner ∼2 kpc; Aladro et al. 2015), NGC 4418 (Costagliola et al. 2015), and Arp 220 (Martín et al. 2011).

In the text

Fig. 7. Comparison of the abundance ratios of APM 08279+5255 (APM 08279 for short), NCv1.143, and three archetype local galaxies – the CMZ of NGC 253 (Martín et al. 2021; Rangwala et al. 2014), ∼2 kpc nuclear region of NGC 1068 (Aladro et al. 2015), and NGC 4418 (Costagliola et al. 2015) – for the species detected in both APM 08279+5255 and NCv1.143. The abundance ratios of the galaxies are computed by dividing the column densities of each molecule by that of C18O, and these ratios are then compared among the galaxies. The horizontal axis is ordered by the increasing value of the ratio APM 08279/NCv1.143 (in the case of C34S and HCN-VIB, where the value of NCv1.143 is missing, we use the ratio of APM 08279/NGC 253 for ordering). Circles are AGN-dominated (APM 08279, NGC 1068, and NGC 4418) over starburst-dominated (NCv1.143 and NGC 253) ratios, diamonds are comparisons between similar types, and crosses are starburst-dominated over AGN-dominated ratios. The upper panel shows the ratio of the two high-redshift galaxies over local galaxies, while the lower panel shows the ratio of AGN-dominated galaxies over starburst-dominated galaxies only.

In the text

	Fig. 8. Energy diagram of an H2O molecule. Black lines indicate the H2O lines detected from our line survey, while gray lines indicate the observed water transitions from other works (see Table 5 for references), with their frequencies indicated (in GHz). The H2O data are taken from the JPL Molecular Spectroscopy website https://spec.jpl.nasa.gov.
In the text

Fig. 9. Spectra of the 380 GHz H2O(414 − 321) line in APM 08279+5255 and NCv1.143 in red (the green error bar indicates the averaged noise level around the spectra per bin), overlaid with the scaled-down profile of the CO(4 − 3) line. The 380 GHz H2O line appears narrower by a factor of a few compared with CO(4 − 3), indicating that it is likely associated with maser origins. Black arrows indicate a ∼2 σ satellite line detection, possibly from the H2O maser, which could originate from accretion disks around the central SMBH.

In the text

Fig. 10. Best-fit model and the likelihood of the parameters from the model of the H2O lines and the dust continuum in APM 08279+5255. Top left: Observed lensing-corrected (black data points) versus modeled H2O integrated flux (green and blue are the extended and compact components, respectively; red is the total). Bottom-left: Lensing-corrected modeled dust SED from the two components (same colors as in the top panel). The photometric data are taken from Leung et al. (2019). Right panels: Likelihood distributions of the key parameters of the H2O model, with the same colors as in the left panels. The arrows indicate the maximum likelihood values (Table 6).

In the text

	Fig. 11. Same as Fig. 10 but for NCv1.143. The photometric data for the dust SED are from Bussmann et al. (2013) and Yang et al. (2017).
In the text

Fig. 12. Violin plots for the LVG-MCMC-derived posterior distributions of nH2, Tkin, and Nmol/dV in APM 08279+5255 and NCv1.143 (upper and lower panels, respectively). The white dot represents the median values, while the solid blue and dashed blue lines represent ±1σ and ±2σ, respectively. The violin shapes represent the posterior distributions of each parameter as presented in the corner plots (Figs. B.1 and B.2).

In the text

	Fig. 13. Fitting of the dense gas SLEDs using LVG models via the MCMC technique. Observed fluxes and their errors are indicated with black symbols. Orange lines show the best fit (Table 7), and light-yellow lines indicate 100 examples within the 15.87% and 84.14% quartiles of the posteriors.
In the text

Fig. 14. Chemical modeling using 3D-DPR to identify the best range of the CR ionization parameter ζCR/nH (x-axis) based on the observed column density ratios. The column densities of HCN (blue), HCO+ (orange), and HNC (green) are normalized with CS. The y-axis shows the ratio of the modeled column density ratio to the observed one. Unity marks agreement between the models and observations. The upper panel shows the results for APM 08279+5255, and the lower panel shows those for NCv1.143. The insets are zoom-in regions of the gray rectangular area in each panel. The thickness of each line represents the uncertainty (light color is for 50% uncertainty, and dark color is for 10% uncertainty).

In the text

Fig. 15. Abundances of HCN, HCO+, HNC and CS versus nH2, for two different CR ionization rates (solid lines for ζCR = 10−13 s−1 and solid-dotted lines for 10−12 s−1). Our model shows that with such a high ionization rate ζCR = 10−13–10−12 s−1, the abundances of HCO+ is expected to be lower than HCN, HNC, and CS in the dense gas regions of nH2 > 105 cm−3. However, with increasing CR ionization rates, at the same H2 densities, the abundance of HCO+ is rising for the dense gas while the abundance of CS is decreasing.

In the text

	Fig. A.1. Same as Fig 3, but zoomed in to show details of the line profiles.
In the text

Fig. A.2. Fitting results of the spectra of APM 08279+5255, with the x-axis as the rest-frame frequency and the y-axis as the observed flux density. The top panel shows the observed spectrum in blue with the uncertainties in the green line. The red solid lines indicate the fit to the entire spectrum, and the dashed red lines are the position of the identified lines. The lower panel shows the residual.

In the text

	Fig. A.3. Same as Fig. A.2 but for NCv1.143.
In the text

	Fig. B.1. Posterior probability distributions of molecular gas density (nH2), gas temperature (Tkin) and column density per velocity (Nmol/dv) of each species in APM 08279+5255 and NCv1.143, with the maximum posterior possibility point in the parameter space shown in orange lines and points. The contours are in steps of 0.5 σ starting from the center.
In the text

	Fig. B.2. Same as Fig. B.1 but for the CS and C34S lines.
In the text

Fig. C.1. NOEMA images of NCv1.143. Color images are the ∼1.8 mm and ∼1.2 mm dust continuum. Yellow, red, and blue contours are the CO(10–9), H2O, and H2O+ lines. The synthesized beams are 0.76″×0.41″and 0.54″×0.30″for the 1.8 and 1.2 mm bands, respectively. The contour levels are: in steps of 1 σ started from ± 3 σ for the H2O(211 − 202) and the H2O+ lines; and in steps of 3 σ started from ± 3 σ for the H2O(321 − 312) and the CO(10–9) lines.

In the text

	Fig. C.2. Posterior possibility distributions of the Reff and μ from the MCMC sampling. The contours are from 1 σ to 3 σ from inside out. The solid blue line, solid red line, dashed purple line, solid yellow line, and dashed black line indicate the results for H2O(211 − 202), H2O(321 − 312), H2O+, CO(10–9), and the dust continuum, respectively.
In the text

Fig. C.3. Lens model of the 254 GHz (top row) and 164 GHz (bottom row) dust continuum emission in NCv1.143. From left to right: (1) Observed image of the 2 mm continuum. (2) Model images with the model parameters of Table . (3) Residuals from the difference between observed and model images. (4) Reconstructed image of the dust continuum emission in the source plane. The critical curves and caustics are shown in red lines in the first three columns and the fourth columns, respectively. The position of the deflector is indicated by red points.

In the text

Fig. C.4. From top to bottom, each row: Lens model of the CO(10–9), H2O(321 − 312) and H2O(211 − 202) lines in NCv1.143. From left to right: (1) Observed CLEANed images of the emission line. (2) Model images with the model parameters of Table . (3) Residuals from the difference between observed and model images. (4) Reconstructed images of the emission line in the source plane. The critical curves and caustics are shown in cyan lines in the first three columns and the fourth column, respectively. The position of the deflector is indicated by the cyan points. One should note that the position of the background source in the source plane is not well-constrained. Thus, the differences seen in the positions of the three lines are not significant, considering the uncertainties.

In the text

	Fig. D.1. 3D-PDR calculations using a density distribution representing the AV, eff − nH relationship (Bisbas et al. 2019). The figure shows the optical depth of 13CO(4–3) and C18O(4–3) lines as a function of the H2 volume density. We find that both lines remain optically thin for nH2 ≲ 106 cm−3), while the 13CO(4–3) starts to have moderate optical depth with increasing H2 densities.
In the text

Fig. E.1. Chemical models for HCO+. Left panel: HCO+ abundance versus the local total H-nucleus number density, nH, for different CR ionization rates. In all cases, the abundance of HCO+ peaks for densities 8 × 103 ≲ nH ≲ 3 × 105 cm−3. Right panel: Temperature profile, with the iso-(ζ/n) lines overplotted in orange, for different CR ionization rates (ζCR = 10−16 − 10−12 s−1, with the same color code as the left panel).

In the text