SUNRISE: The rich molecular inventory of high-redshift dusty galaxies revealed by broadband spectral line surveys

Understanding the nature of high-redshift dusty galaxies requires a comprehensive view of their interstellar medium (ISM) and molecular complexity. However, the molecular ISM at high redshifts is commonly studied using only a few species beyond 12 C 16 O, limiting our understanding. In this paper, we present the results of deep 3mm spectral line surveys using the NOrthern Extended Millimeter Array (NOEMA) targeting two strongly lensed dusty galaxies observed when the Universe was less than 1.8Gyr old: APM08279 + 5255, a quasar at redshift z = 3.911, and NCv1.143 ( H -ATLAS J125632.7 + 233625), a z = 3.565 starburst galaxy. The spectral line surveys cover rest-frame frequencies from about 330 to 550GHz for both galaxies. We report the detection of 38 and 25 emission lines in APM08279 + 5255 and NCv1.143, respectively. These lines originate from 17 species, namely CO, 13 CO, C 18 O, CN, CCH, HCN, HCO + , HNC, CS, C 34 S, H 2 O, H 3 O + , NO, N 2 H + , CH, c-C 3 H 2 , and the vibrationally excited HCN and neutral carbon. The spectra reveal the chemical richness and the complexity of the physical properties of the ISM. By comparing the spectra of the two sources and combining the analysis of the molecular gas excitation, we ﬁnd that the physical properties and the chemical imprints of the ISM are di ﬀ erent: the molecular gas is more excited in APM08279 + 5255, which exhibits higher molecular gas temperatures and densities compared to NCv1.143; the molecular abundances in APM08279 + 5255 are akin to the values of local active galactic nuclei (AGN), showing boosted relative abundances of the dense gas tracers that might be related to high-temperature chemistry and / or the X-ray-dominated regions, while NCv1.143 more closely resembles local starburst galaxies. The most signiﬁcant di ﬀ erences between the two sources are found in H 2 O: the 448GHz ortho-H 2 O(4 23 − 3 30 ) line is signiﬁcantly brighter in APM08279 + 5255, which is likely linked to the intense far-infrared radiation from the dust powered by AGN. Our astrochemical model suggests that, at such high column densities, far-ultraviolet radiation is less important in regulating the ISM, while cosmic rays (and / or X-rays and shocks) are the key players in shaping the molecular abundances and the initial conditions of star formation. Both our observed CO isotopologs line ratios and the derived extreme ISM conditions (high gas temperatures, densities, and cosmic-ray ionization rates) suggest the presence of a top-heavy stellar initial mass function. From the ∼ 330–550GHz continuum, we also ﬁnd evidence of nonthermal millimeter ﬂux excess in APM08279 + 5255 that might be related to the central supermassive black hole. Such deep spectral line surveys open a new window into the physics and chemistry of


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 3 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(Tielens , 2021)).Consequently, a comprehensive understanding of the ISM in galaxies requires ob-⋆ We dedicate this paper to the memory of our coauthor and friend, Yu Gao, who passed away in May 2022.⋆⋆ The final data products of the tables derived from UVFIT will be available in electronic form at the CDS via anonymous ftp to cdsarc.cds.unistra.fr(130.79.128.5) or via https://cdsarc.cds.unistra.fr/cgi-bin/qcat?J/A+A/.
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 12 C 16 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. 2011Muller et al. , 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 2 O, H 2 O + , 13 CO, 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. 2006Riechers et al. , 2007Riechers et al. , 2010Riechers et al. , 2011b;;Danielson et al. 2011;Oteo et al. 2017;Béthermin et al. 2018;Cañameras et al. 2021;Rybak et al. 2022Rybak et al. , 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 −10 (L IR /L ⊙ ) M ⊙ yr −1 .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 −10 (L IR /L ⊙ ) M ⊙ yr −1 (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.
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.

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 15 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 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 14 µ −1 L ⊙ (Leung et al. 2019, note the different definition of L IR1 ).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 Article number, page 2 of 39 C. Yang et al.: SUNRISE: Submillimeter molecUlar liNe suRveys in dIstant duSty galaxiEs be attributed to two prominent excitation components -a cold one with T kin ∼ 65 K and n H 2 ∼ 10 5 cm −3 , and a warm one with T kin ∼ 220 K and n H 2 ∼ 10 4 cm −3 which is probably heated by the central AGN (Weiß et al. 2007).Observations of multiple H 2 O lines in APM 08279+5255 suggest that the submillimeter H 2 O emission originates in similar regions as CO with a dust temperature T d ∼ 220 K and the optical depth at 100 µm, τ 100 ≈ 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 1 of L IR = (1.4 ± 0.2)×10 14 µ −1 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 (2017b).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 13 L ⊙ with a molecular gas reservoir of (6.2 ± 1.6)×10 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 Σ L IR ∼ ×10 13 L ⊙ kpc −2 .A series of extensive followup studies (e.g., Yang et al. 2016Yang et al. , 2017) ) found no clear AGN signatures, indicating that this galaxy is likely dominated by starburst activity with an estimated SFR = 1.3 × 10 3 M ⊙ yr −1 , 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 H 2 ∼ 10 4.1 cm −3 , and a high-excitation molecular gas component with T kin ∼ 63 K and n H 2 ∼ 10 4.2 cm −3 (Yang et al. 2017).Detection of the J up = 2, 3 and 4 H 2 O lines in NCv1.143 suggests the existence of an optically thick (τ 100 > 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 2017b).The similarities between the line profiles of the high-J CO and H 2 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 2 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 2017b) 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.

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. and for NCv1.143,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.
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.7sDEC +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  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.
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 arbituary unit 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. the tunings was ∼ 22.5 hours for NCv1.143 and ∼ 16.4 hours 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.0mJy/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 GILDAS2 packages CLIC, with the native 2 MHz spectral resolution.Then imaging, CLEANing, uv fitting,  1).The ∼ 9 ′′ shift in δ RA of APM 08279+5255 was due to a shifted observational phase center.The δ RA and δ DEC positions are the center of the Gaussian profiles relative to the phase center.The rest of the fitted parameters -the major and minor axis of the ellipse, r major and r minor and the position angle PA are all consistent across different spectral windows within the uncertainties.As they are irrelevant to the purpose of this work, we do not list them here.and spectra extraction were performed with GILDAS's MAPPING on the uv tables.

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 signalto-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 2 O line images of NCv1.143 (see Yang 2017b 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.

Extracting spectra
For flux extraction, we assume a model of a single twodimensional 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 twodimensional 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.
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.

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 MADCUBA3 (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.

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 The model adopts a linewidth of 520 and 330 km s −1 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 2 , 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 18 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 18 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 spec-tral 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 2 O and the vibration excited HCN), for APM 08279+5255 and NCv1.143, respectively.It is important to note that H 2 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 2 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.

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

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.Note: LTE-derived column densities and excitation temperatures of APM 08279+5255 (APM for short) and NCv1.143 (NC for short).We note that the values are not corrected for lensing magnification.Temperature values without errors indicate species for which this parameter could not be constrained and was fixed to the value derived from HCN and HCO + .Although most of the lines are well modeled with a single temperature, for APM 0879+5255, the vibrationally excited HCN (v 2 = 1 f ) requires significantly higher kinetic temperatures.The same applies to the H 2 O lines, where a pure collisional excitation cannot explain the fluxes well.Values in brackets are errors. (a) The reported column density log 10 N = A(B) corresponds to N = 10 A ± 10 B .
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: 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 (ν 0 = 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    For the CN lines, we do not distinguish between the hyper-fine splitting caused by the nuclear spin of nitrogen (described by quantum number F) due to limited spectral resolution.For the CCH lines, because the frequencies of the fine and hyper-fine splitting within the same N +1→N are very close (considering the four brightest transitions for a given N J,F , we have the brightest four lines that are within 0.06 GHz), we do not distinguish the hyper-fine lines.For the NO and CH, we only list the fine structure and Λ-doublet because the hyper-fine lines are too close to be separated.
For the aforementioned molecules, we list the E up /k B and ν rest−frame with the highest Einstein A coefficient as a representative (for the CH lines, we list 3/2, 1 + →1/2, 0 − and 3/2, 2 − →1/2, 1 + ).The data of C 34 S, HCN-VIB (v 2 = 1 f ) and c-C 3 H 2 are taken from CDMS (Endres et al. 2016).For the HCN and the N 2 H + lines, we do not consider their hyper-fine structure lines considering the spectral resolution and the intrinsic linewidths of the galaxies.
0.00 ± 0.02 mJy for NCv1.143 and C R−J = 0.15 ± 0.01 mJy and C free−free = 0.38±0.01mJy 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.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.
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 4 µ −1 M ⊙ yr −1 .Using the correlation between freefree luminosity and the SFR (Algera et al. 2022), assuming the typical electron temperature of the HII regions of 10 4 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 31 erg s −1 Hz −1 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 46 erg s −1 ; 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.2mJy 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.).

A detailed look at the detected molecules and a comparison of the two sources
As displayed in Fig. 3 and )   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).
From the observed integrated fluxes in Table 3, we derive the apparent line luminosities, µL line (in units of L ⊙ ) and µL ′ line (in units of K km s −1 pc 2 ) following Solomon et al. (1992) where I line , ν and D L are in the units of Jy km s −1 , GHz and Mpc, respectively.We note here that because the 13 CO(4-3) and CS  lines are blended due to their close frequencies, we extrapolate the integrated flux of CS  by combining the data of CS(7-6), CS  and CS(10-9), and we subtracted this flux from the blended line flux for obtaining the integrated flux of 13 CO .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 13 CO(4-3) line in APM 08279+5255.(Costagliola et al. 2015), and Arp 220 (Martín et al. 2011).
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 18 O(4-3) line (given that the C 18 O line is likely optically thin as we argue in Appendix D and the emitting size of the C 18 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 18 O(4-3) (L line and L ′ line ) 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 18 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 18 O(4-3) compared to NCv1.143, while the 13 CO(4-3)/C 18 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 18 O(4-3).If there is no deficit of C 18 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 6 K cm −3 (Weiß et al. 2007).While in NCv1.143, the thermal pressure is ≲ 1 × 10 6 K cm −3 (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.
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 ′ line (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 ex-treme 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 2 O lines, where we find slightly elevated H 2 O(2 11 -2 02 ) and H 2 O(3 21 -3 12 ) luminosity in NCv1.143, while the 448 GHz H 2 O(4 23 -3 30 ) is about nine times brighter in APM 08279+5255.Such a difference suggests that the J up = 2 and 3 H 2 O lines might arise from very different conditions compared to the region that the H 2 O(4 23 -3 30 ) traces.The latter might be tightly related to AGNpowered 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 2 O excitation in Sect.6.3.Besides H 2 O, the H 3 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.

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 ) 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.
minimize the chemical variations due to differences of spatial scales (Butterworth et al. 2022).
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 18 O (as a proxy of total H 2 , 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 18 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 wellmixed.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 highspatial-resolution data are available.
The comparison of the relative abundance (with respect to C 18 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 34 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 com-pared 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.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, 13 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. 2010a), 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 34 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.143and 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.143and 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.143and APM 08279+5255/NGC 253) generally exhibit higher values that are mostly above unity, whereas the starburst/AGN ratios (NCv1.143/NGC4418 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 2 H + , and C 34 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 18 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. (2010a) 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.

The [C I] line
Both galaxies are detected in the neutral carbon [C I] 3 P 1 -3 P 0 fine structure line at 492 GHz ([C I](1-0) hereafter), with luminosities µL ′ [CI] of 3.55×10 10 K km s −1 pc 2 and 10.1×10 10 K km s −1 pc 2 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 24 cm −2 (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 2 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): 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 ⊙ .

Submillimeter H 2 O lines
Thanks to the observations of extragalactic H 2 O lines done by Herschel (Pilbratt et al. 2010) and ground-based telescopes (e.g., Omont et al. 2011Omont et al. , 2013;;Yang et al. 2013Yang et al. , 2016)), it is now established that H 2 O is one of the most important interstellar molecules after CO.It is because the submillimeter H 2 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 2 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 2 O has a central role in the oxygen chemical networks in the ISM, where H 2 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 2 O spectra at rest-frame frequencies below 550 GHz, where there are only three relatively strong H 2 O lines available -H 2 O(4 14 -3 21 ), H 2 O(4 23 -3 30 ) and H 2 O(5 33 -4 40 ): all three are detected in APM 08279+5255 and two of them, H 2 O(4 14 -3 12 ) and H 2 O(4 23 -3 30 ), in NCv1.143 (Tables 3 and 5).This number may appear marginal compared to the typical strong submillimeter H 2 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.( Impellizzeri et al. 2008).The H 2 O(4 14 -3 21 ) 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.
This ortho-H 2 O(4 14 -3 21 ) line at 380.197 GHz (380 GHz hereafter), with its metastable backbone level J K a ,K c = 4 14 (E up = 324 K), is one of the most promising H 2 O maser lines predicted by radiative transfer models (e.g., Neufeld & Melnick 1991;Gray et al. 2016).The detection of the 380 GHz H 2 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 2 O is the only submillimeter H 2 O maser detected in both sources in the rest-frame frequency range of 330-550 GHz.Moreover, the 380 GHz H 2 O maser is at least four times stronger than its 22 GHz H 2 O(6 16 -5 23 ) 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 2 O maser line.It is thus one of the best extragalactic H 2 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 −1 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 2 O line from a local survey of the H 2 O megamasers disks (Kuo et al. 2011), indicating that the accretion disk, in which the 380 GHz H 2 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 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.

The 448 GHz H 2 O(4 23 -3 30 )
The ortho-H 2 O(4 23 -3 30 ) line at 448.001 GHz is predominantly excited through the absorption of 79 and 132 µm photons in the transitions 4 23 ← 3 12 and 4 23 ← 4 14 (Fig. 8).The influence of collisional excitation on the 448 GHz H 2 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 23 level is thus likely mostly radiatively populated.Recently, this key H 2 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 2 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 2 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 2 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 24 cm −2 .Comparing the ratio between the 448 GHz H 2 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 Xray 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.The para-H 2 O(5 33 -4 40 ) at rest frame 474.689GHz (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 2 O line transition 5 23 -5 14 at 1411 GHz (E up = 642 K) in local galaxies (Yang et al. 2013).Observations of evolved stars found the 475 GHz H 2 O lines are masers associated with fast outflows (Menten et al. 2008).
In dusty warm dense regions of galaxies, H 2 O molecules can also absorb 53 µm far-infrared photons through 5 33 ← 4 22 (Fig. 8).Such absorption features have been detected in NGC 4418 (González-Alfonso et al. 2012).In this case, the high H 2 O column density can efficiently populate the 4 22 level and thus sustain this far-infrared pumping process.Such a process may potentially provide a route for the formation of H 2 O maser.Nevertheless, the observed line profile of H 2 O(5 33 -4 40 ) in APM 08279+5255 is similar to those of other dense gas tracers, suggesting the H 2 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 2 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.

Model of the H 2 O emission lines and its interpretation
Given that in both APM 08279+5255 and NCv1.143, there are detections of other H 2 O lines from the literature (Lis et al. 2011;van der Werf et al. 2011;Yang et al. 2016;Yang 2017b), it is worth combining all the detected H 2 O lines together with the dust continuum to constrain the excitation conditions.For APM 08279+5255, there are eight H 2 O lines from J up = 2 to 5 included, while for NCv1.143, six H 2 O lines are included (Table 5, Figs. 5 and 8).Here, we assume that the H 2 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) and 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 2 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 τ 100 , the col-umn density of H 2 O along a radial path N H 2 O , the H 2 density n H 2 , 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 , τ 100 , N H 2 O , n H 2 ).However, it should be noted that the conditions of the gas that the H 2 O lines trace can be different from the CO-traced conditions.We keep fixed ∆V = 100 km s −1 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 2 O is sensitive to the dust SED, and thus, our goal is to fit the H 2 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 2 O lines (and the upper limit of the H 2 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, χ 2 is minimized for all possible combinations among the N C components, giving for each one the squared radius R 2 of each component.The best-fit (fiducial) model is the combination that yields the minimum χ 2 , and the regular grid enables the estimation of the likelihood of each parameter.The derived parameters -R, the H 2 O abundance relative to H nuclei X H 2 O , 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 2 O and continuum emission with a single-model component (N C = 1), but results were unreliable, with a best-reduced χ 2 value χ 2 red ≈ 4.This was indeed expected because the low-lying H 2 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 χ 2 red ≈ 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 τ 100 ∼ 2, with warm dust temperatures of over 100 K (Figs. 10 and Fig. 11).In NCv1.143, this core has an intrinsic luminosity of L IR ∼ 4 × 10 12 L ⊙ , resulting in an extreme infrared luminosity surface density Σ IR ≡ L IR /(4π R 2 ) = 1.4 × 10 13 L ⊙ kpc −2 .This translates into a surface SFR value of SFR/(π R 2 ) ∼ 6.2 × 10 3 M ⊙ yr −1 kpc −2 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 14 L ⊙ that results in Σ IR = 3 × 10 14 L ⊙ kpc −2 .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 2 O( 2 6).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 12 L ⊙ and 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 2 O(5 33 -4 40 ) 475 GHz line.This means that the H 2 O lines are indeed tracing the shape of the SED, and high columns of H 2 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 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 2 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.

H 3 O + lines
As a symmetric top molecule, H 3 O + plays an important role in the oxygen chemistry network that is related to the H 2 O formation in ion-neutral reactions (Gerin et al. 2016) Sakamoto et al. 2021).H 3 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 3 O + with H 2 O and HCO + and found the chemistry in XDRs can explain some of the high column density ratios of H 3 O + over H 2 O and HCO + in some sources (Aalto et al. 2011).
lines are robustly detected in APM 08279+5255.These are the first high-redshift detections of any H 3 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 3 O + is about N H 3 O + ∼ 10 15 cm −2 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 photodissociationdominated region (PDR) can explain the column density (Aalto et al. 2011).
The abundance ratio between H 3 O + and H 2 O can also be used as a probe of ionization states.Using the H 2 O column density derived from Sect.6.3.4 of 2.5 × 10 17 cm −2 , we find a relative abundance ratio of H 3 O + /H 2 O ∼ 3 × 10 −3 , which can be explained by both the PDR and the XDR (Riechers et al. 2009).Indeed, H 3 O + lines have also been detected in the pure-starburst NGC 253 (Martín et al. 2021;Holdship et al. 2022a).While the aforementioned diagnostics are inconclusive, the most evident abundance ratio that suggests the XDR plays an important role comes from H 3 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 3 O + abundances in APM 08279+5255, though the possibility of H 3 O + enhancement by CR cannot be entirely ruled out.6.5.Typical high dipole moment tracers -HCN, HCO + , HNC, CS and C 34 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 et al. 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. 2007Riechers et al. , 2011b;;Carilli & Walter 2013;Oteo et al. 2017;Béthermin et al. 2018;Cañameras et al. 2021;Rybak et al. 2022Rybak et al. , 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 34 S, NO, and N 2 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. 13 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 34 S for a detailed non-LTE analysis for the ISM conditions (N 2 H + is excluded in the LVG analysis due to the lack of collisional rates data).
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 34 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 2 (n H 2 ), 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: Physical parameters derived from the LVG models.Table 7 lists the results from the LVG-MCMC analysis, while the fitted SLEDs are shown in Fig. 12.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  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 distributions of the two sources using the violin plots in Fig. 13.It becomes apparent that all the n H 2 values derived from HCN, HCO + , HNC, CS, and C 34 S significantly exceed that of the CO from previous studies (Weiß et al. 2007 andYang et al. 2017, n H 2 ∼ 10 3 -10 4 cm −3 ), falling within the range of 10 5 -10 7 cm −3 .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 14 -10 15 cm −2 / km s −1 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. 2016a).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. 2010b).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.2toward 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 44 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. 2016a).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. (2016a) 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.3  −1.1 and −0.1 +2.5 −2.7 , 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.1 −1.3 and 0.6 +2.0 −2.4 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).

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

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 −3 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(Muller et al. , 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(Rangwala et al. , 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) tran-sition 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 −5 .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 −4 (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 3 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 assump-tion.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: 13 CO and C 18 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 13 CO/C 18 O may be assumed to be close to the isotope ratio 13 C/ 18 O if both lines of 13 CO and C 18 O are optically thin (however, we caution that in very high-density gas, the 13 CO lines may have a moderate optical depth; see Appendix D).The isotopes, 13 C and 18 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 13 CO/C 18 O, which represents their abundance ratio, can serve as a dust-insensitive stellar IMF probe.Zhang et al. (2018) demonstrated that low flux ratios 13 CO/C 18 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/ 13 CO, specifically those exceeding 20.Note: APM stands for APM 08279+5255 and NC for NCv1.143.Zhang et al. (2018) present deep observations of the three isotopologs -CO, 13 CO and C 18 O -toward four high-redshift lensed galaxies, including three starbursts and one quasar.They found flux ratios 13 CO/C 18 O and CO/ 13 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 13 CO(4-3) line is blended with the CS(9-8) line, which is blueshifted with respect to 13 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  line.The flux of this line may be estimated by interpolation between the observed values of the CS  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 -13 CO/C 18 O < ∼ 1 and CO/ 13 CO > ∼ 20are 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 2 = 1 f of HCN (HCN-VIB hereafter) at high redshifts.In the APM 08279+5255 quasar, all three transitions of 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 24 -10 25 cm −2 , 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;González-Alfonso & Sakamoto 2019;Falstad et al. 2019;Aalto et al. 2019;Falstad et al. 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 7 L ⊙ for APM 08275+5255.This yields a luminosity ratio of L HCN-VIB(4-3) /L IR = (2.9±1.2)×10−7 .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 −8 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 −1 results in a limit to its luminosity of L HCN-VIB(4-3) < 1.1 × 10 7 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 2 O excitation model in Sect.6.3.4.Nevertheless, the upper limit of the ratio L HCN-VIB(4-3) /L IR < 9 × 10 −7 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 −2 .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 −2 , 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. (2016b), 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 2 H + , and c-C 3 H 2 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 2 (Harada et al. 2010a).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 Fig. 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 5 cm −3 .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 2 H + , confirming the first tentative detection of the N 2 H + (5-4) line report by Feruglio et al. (2017).Also, we report a tentative detection of the cyclopropenylidene c-C 3 H 2 line at 485.7 GHz for the first time at high redshifts, which needs further confirmation.

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 = 105 cm −3 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 5 χ 0 5 , 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-PDR6 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 Dunne et al. 2022) and a microturbulent velocity of v turb = 1 km s −1 are further assumed.
In addition, we adopted the following expression of optical depths, which best describes macroturbulent molecular clouds (Papadopoulos & Seaquist 1999): where In the above equations, A i j is the Einstein coefficient of levels i and j, ν i j 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 √ α = 1 (see also Bisbas et al. 2015).Equation 4is 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.
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 = R N,sim /R N,obs versus ζ CR /n H , where R 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) Fig. 14: Chemical modeling using 3D-DPR to identify the best range of the CR ionization parameter ζ CR /n H (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).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 −14 s −1 for APM 08279+5255 (although the high ζ CR /n H solution of the HCN/CS ratio suggests a higher value of ζ CR ∼ 10 −12 s −1 ), and of ζ CR ≃ (8 − 12) × 10 −13 s −1 for the NCv1.143galaxy.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 highredshift galaxies -Bremsstrahlung radiation (Yusef-Zadeh et al. 2013) Imanishi et al. (2018) studied OH + and H 2 O + lines of the extended halos of z ∼ 2.3 galaxies and derived a ζ CR ∼ (1 − 100) × 10 −16 s −1 , 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 −14 -10 −12 s −1 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.
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 H 2 ∼ 10 5 cm −3 and ζ CR ≃ 10 −13 s −1 , 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. 13, 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 −12 s −1 , 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 simi- 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 n H 2 > 10 5 cm −3 .However, with increasing CR ionization rates, at the same H 2 densities, the abundance of HCO + is rising for the dense gas while the abundance of CS is decreasing.
lar relative abundances among the dense gas tracers, consistent with Fig. 13.It is also worth noticing that, given such a CR ionization rate, the relative abundance of HCO + drops toward the denser conditions from n H 2 ∼ 10 5 cm −3 , while CS abundance increases with increasing gas densities.Around the gas densities of n H 2 ∼ 10 5-6 cm −3 , 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 −2 .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 2 density of ≳ 10 5 will lead to a stellar IMF peaking at ≳ 15 M ⊙ (Klessen et al. 2007).Moreover, such a topheavy 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 18 O abundance (Danielson et al. 2013).

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. 2011a;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-angularresolution 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 2 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 2 O(2 11 -2 02 ) and H 2 O(3 21 -3 12 ) emission.Because the line profiles are not strongly distorted by the magnification, and the resolved line profiles of mid/high-J CO and H 2 O(2 11 -2 02 ) and H 2 O(3 21 -3 12 ) 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 2 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.

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)  -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 starburstdominated 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 hottemperature 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 34 S, which are trace much denser gas than CO: in the range n H 2 ∼ 10 5 -10 7 cm −3 .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 2 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 2 O megamaser candidate in high-redshift galaxies.Simply scaling from local galaxies, we crudely derived the SMBH mass, which is about 10 10 M ⊙ , consistent with previous literature values.
-After including the two water emission lines detected here for the first time, we built a H 2 O excitation model using multiple transitions of H 2 O and the dust continuum.We find the H 2 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 2 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 3 O + (detected in both APM 08279+5255 and NCv1.143, the first highredshift detections) and HCO + , where we find that XDRs play an important role.
-We have also detected the isotopologs of CO ( 13 CO and C 18 O) in both sources.We find similar flux ratios of 13 CO/C 18 O and 12 CO/ 13 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 dustobscured 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 −13 to 10 −12 s −1 , 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 18 O/ 13 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.
. The coordinates and redshifts of APM 08279+5255 and NCv1.143 are taken fromYang et al. (2017) andWeiß et al. (2007), respectively.The tuning names also correspond to the Prog name that can be found on the Science Data Archive of IRAM (http://vizier.u-strasbg.fr/viz-bin/VizieR-3?-source=B/iram/noema).The synthesized beam sizes correspond to natural weighting.The RMS /2 MHz is the root mean square noise per 2 MHz bandwidth of the spectra for each sideband of each tuning.

Fig. 3 :
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.

Fig. 4 :
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.

Fig. 5 :
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 pc 2 , L ′ line , normalized by the value of C 18 O(4-3).The lower panel shows the line luminosity in L ⊙ , L line , normalized by the C 18 O(4-3) line.The lowest levels of the transitions of each molecule detected are labeled at the top of each panel, except for H 2 O and H 3 O + , 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.

Fig. 7 :
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 C 18 O, 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 C 34 S 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.143and 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.

6. 3 Fig. 9 :
Fig. 9: Spectra of the 380 GHz H 2 O(4 14 -3 21 ) 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 H 2 O 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 H 2 O maser, which could originate from accretion disks around the central SMBH.

Fig. 10 :
Fig.10: Best-fit model and the likelihood of the parameters from the model of the H 2 O lines and the dust continuum in APM 08279+5255.Top left: Observed lensing-corrected (black data points) versus modeled H 2 O 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 fromLeung et al. (2019).Right panels: Likelihood distributions of the key parameters of the H 2 O model, with the same colors as in the left panels.The arrows indicate the maximum likelihood values (Table6).

Fig. 12 :
Fig.12: 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 (Table7), and light-yellow lines indicate 100 examples within the 15.87% and 84.14% quartiles of the posteriors.

Fig. 13 :
Fig. 13: Violin plots for the LVG-MCMC-derived posterior distributions of n H 2 , T kin , and N mol /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).

3 -
Note: LVG-MCMC-derived parameters; the median values and the ± 1 σ uncertainties are noted.The results where only two lines were detected in NCv1.143 are also included despite their large uncertainties, which are well captured by the posterior distributions shown in the corner plots (Figs.B.1 and B.2).For comparison, we also include the LVG model results of CO, using the line fluxes fromWeiß et al. (2007) andYang et al. (2017), which trace much less extreme physical conditions compared to the dense gas tracers.

Fig. 15 :
Fig.15: Abundances of HCN, HCO + , HNC and CS versus n H 2 , 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 n H 2 > 10 5 cm −3 .However, with increasing CR ionization rates, at the same H 2 densities, the abundance of HCO + is rising for the dense gas while the abundance of CS is decreasing.

Fig. A. 2 :
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.

Table 2 :
Parameter of the elliptical Gaussian model used in fitting the uv data.

Table 3 :
Detected species and transitions in the two sources from our NOEMA surveys.

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 2 = 1 f ), 380 GHz H 2 O maser, and H 3 O + lines are the first high-redshift detections in individual sources reported.

Table 4
Aladro et al. 2015)Rangwala et al. 2014)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 Papadopoulos et al. (2022) 2 O molecule.Black lines indicate the H 2 O lines detected from our line survey, while gray lines indicate the observed water transitions from other works (see Table5for references), with their frequencies indicated (in GHz).The H 2 O data are taken from the JPL Molecular Spectroscopy website https://spec.jpl.nasa.gov.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 byPapadopoulos 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 −5

Table 6 :
Physical parameters derived from H 2 O modeling of APM 08279+5255 and NCv1.143.Note: a : 90% confidence intervals; b : Values for the fiducial best-fit model, selected for detailed comparison with data, while the uncertainties of the derived parameters are well-characterized in Figs. 10 and 11.
. First detections of the H 3 O + lines were achieved in two local starburst galaxy, M 82 and Arp 220, from the transition of p-H 3 O + (3 + 2 -2 − 2 ) at 364.797 GHz, with E up /k B = 139 K (van der Tak et al. 2008, see also

Table 8 :
Flux ratios between CH versus HCN and HCO + .

Table 10 :
3D-PDR model results for APM 08279+5255 and NCv1.143.The ζ CR value of the third column was calculated using the n H density of CS from the LVG results shown in Table7.
Holdship et al. (2021)etit et al. 2016);Le Petit et al. 2016)in the Galactic Center suggest a CR ionization rate ζ CR ≃ (1 − 100) × 10 −15 s −1 , 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 −14 s −1 derived from the observed HCO + /HOC + emission ratio(Harada et al. 2021).Similarly,Holdship et al. (2021)report an ionization rate ζ CR ∼ 10 −11 s −1 derived from chemical models of C 2 H.However, in a subsequent paper and using H 3 O + and SO observations, Holdship et al. (2022b) suggest a lower value of ζ CR ≃ (1 − 80) × 10 −14 s −1 for NGC 253.Very high ζ CR values have also been derived for other star-forming galaxies.OH + , H 2 O + and H 3 O + observations in Arp 220 find a ζ CR /n H ≃ (1 − 2) × 10 −17 cm 3 s −1 , which leads to a ζ CR ≳ 10 −13 s −1 for an estimated n H ∼ 10 4 cm −3 (González-Alfonso et al. 2013).Similarly, the aforementioned lines observed for Mrk 231 also suggest a ζ CR plus a matchedfiltering method, we detect 38 emission lines in APM 08279+5255 and 25 emission lines in NCv1.143, from [C I], CO, 13 CO, C 18 O, CN, CCH, HCN, HCO + , HNC, CS, C 34 S, H 2 O, H 3 O + , NO, N 2 H + , CH, and c-C 3 H 2 , as well as the rotational levels of the vibrationally excited HCN (HCN-VIB).The CH, CCH, c-C 3 H 2 , N 2 H + , HCN-VIB, and H 3 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 2 O line, where APM 08279+5255 is about nine times brighter than NCv1.143 in flux.