# Diverse Molecular Gas Excitations in Quasar Host Galaxies at $z \sim 6$

Jianan Li, Ran Wang, Antonio Pensabene, Fabian Walter, Bram P.

Venemans, Roberto Decarli, Eduardo Bañados, Pierre Cox, Roberto Neri, Alain Omont, et al.

► To cite this version:

Jianan Li, Ran Wang, Antonio Pensabene, Fabian Walter, Bram P. Venemans, et al. Diverse Molecular Gas Excitations in Quasar Host Galaxies at $z \sim 6$. The Astrophysical Journal, 2024, 962, 10.3847/1538-4357/ad1754. insu-04851593

HAL Id: insu-04851593

https://insu.hal.science/insu-04851593v1

Submitted on 20 Dec 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Diverse Molecular Gas Excitations in Quasar Host Galaxies at $z \sim 6$

Jianan Li<sup>1,2</sup>①, Ran Wang<sup>2</sup>①, Antonio Pensabene<sup>3,4</sup>①, Fabian Walter<sup>4,5</sup>①, Bram P. Venemans<sup>6</sup>①, Roberto Decarli<sup>7</sup>①, Eduardo Bañados<sup>4</sup>①, Pierre Cox<sup>8</sup>①, Roberto Neri<sup>9</sup>①, Alain Omont<sup>10</sup>①, Zheng Cai<sup>1</sup>①, Yana Khusanova<sup>4</sup>①, Fuxiang Xu<sup>2,11</sup>①, Dominik Riechers<sup>12</sup>①, Jeff Wagg<sup>13</sup>①, Yali Shao<sup>14</sup>①, Yuanqi Liu<sup>15</sup>①, Karl M. Menten<sup>16</sup>①, Qiong Li<sup>17</sup>①, and Xiaohui Fan<sup>18</sup>①

<sup>1</sup>Department of Astronomy, Tsinghua University, Beijing 100084, People's Republic of China; rwangkiaa@pku.edu.cn

<sup>2</sup>Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, People's Republic of China; jianan1@mail.tsinghua.edu.cn

<sup>3</sup>Dipartimento di Fisica “G. Occhialini,” Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126, Milano, Italy

<sup>4</sup>Max Planck Institut für Astronomie, Königstuhl 17, D-69117, Heidelberg, Germany

<sup>5</sup>National Radio Astronomy Observatory, Pete V. Domenici Array Science Center, P.O. Box O, Socorro, NM 87801, USA

<sup>6</sup>Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands

<sup>7</sup>INAF—Osservatorio di Astrofisica e Scienza dello Spazio, via Gobetti 93/3, I-40129 Bologna, Italy

<sup>8</sup>Sorbonne Université, UPMC Université Paris 6 and CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98bis Boulevard Arago, F-75014 Paris, France

<sup>9</sup>Institut de Radioastronomie Millimétrique (IRAM), 300 Rue de la Piscine, F-38400 Saint-Martin-d'Hères, France

<sup>10</sup>Institut d'Astrophysique de Paris, Sorbonne Université, CNRS, UMR 7095, 98 bis bd Arago, F-75014 Paris, France

<sup>11</sup>Department of Astronomy, School of Physics, Peking University, Beijing 100871, People's Republic of China

<sup>12</sup>I. Physikalisches Institut, Universität zu Köln, Zülpicher Strasse 77, D-50937 Köln, Germany

<sup>13</sup>PIFI Visiting Scientist, Purple Mountain Observatory, No. 8 Yuanhua Road, Qixia District, Nanjing 210034, People's Republic of China

<sup>14</sup>School of Space and Environment, Beihang University, Beijing, People's Republic of China

<sup>15</sup>Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030, People's Republic of China

<sup>16</sup>Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany

<sup>17</sup>Jodrell Bank Centre for Astrophysics, University of Manchester, Oxford Road, Manchester M13 9PL, UK

<sup>18</sup>Steward Observatory, University of Arizona, 933 N Cherry Avenue, Tucson, AZ 85721, USA

Received 2023 September 15; revised 2023 December 13; accepted 2023 December 18; published 2024 February 13

## Abstract

We present observations using the Northern Extended Millimetre Array (NOEMA) of CO and H₂O emission lines and the underlying dust continuum in two quasars at $z \sim 6$, i.e., P215–16 at $z = 5.78$ and J1429+5447 at $z = 6.18$. Notably, among all published CO spectral line energy distributions (SLEDs) of quasars at $z \sim 6$, the two systems reveal the highest and the lowest CO levels of excitation, respectively. Our radiative transfer modeling of the CO SLED of P215–16 suggests that the molecular gas heated by active galactic nuclei could be a plausible origin for the high CO excitation. For J1429+5447, we obtain the first well-sampled CO SLED (from transitions from 2–1 to 10–9) of a radio-loud quasar at $z \gtrsim 6$. Analysis of the CO SLED suggests that a component from a single photodissociation region could explain the CO excitation in the radio-loud quasar J1429+5447. This work highlights the utility of the CO SLED in uncovering the properties of the interstellar medium in these young quasar-starburst systems at the highest redshift. The diversity of the CO SLEDs reveals the complexities in gas conditions and excitation mechanisms at their early evolutionary stage.

*Unified Astronomy Thesaurus concepts: Submillimeter astronomy (1647); AGN host galaxies (2017); Molecular gas (1073); Early universe (435)*

## 1. Introduction

Active star formation has been detected in the host galaxies of quasars from the local to the distant Universe, revealing coevolution of the supermassive black holes (SMBHs) and galaxies (e.g., Kormendy & Ho 2013; Wang et al. 2013, 2016; Willott et al. 2015; Venemans et al. 2016; Decarli et al. 2017; Pensabene et al. 2020). The cold (~100 K) interstellar medium (ISM) in the quasar host galaxies fuels both the active galactic nucleus (AGN) and the nuclear star formation. Observations of the molecular and atomic lines at submillimeter/millimeter wavelengths provide rich information about the physical conditions and heating power of the cold ISM, which are the keys to understanding the evolution of the young quasar hosts (e.g., Wang et al. 2011, 2013, 2016; Decarli et al. 2017, 2018, 2023; Venemans et al. 2018; Shao et al. 2019, 2022; Li et al. 2020a, 2020b; Meyer et al. 2022).

The CO molecule, as the second most abundant species of the molecular ISM, gives the brightest molecular line emission in the (sub)millimeter regime. The CO spectral line energy distributions (SLEDs), i.e., the CO line fluxes as a function of rotational transition quantum number *J*, are sensitive diagnostics of the physical properties of the molecular gas (e.g., temperature, density, and radiation field strength). The low-*J* (*J* ≤ 3) CO emission lines directly constrain the total molecular gas mass of galaxies (e.g., Walter et al. 2003; Riechers et al. 2006, 2011; Wang et al. 2011; Shao et al. 2019; Ramos Almeida et al. 2022; Montoya Arroyave et al. 2023). Low-*J* CO lines are found to be sublinearly correlated with far-infrared (FIR) luminosity (e.g., Riechers et al. 2006; Riechers 2011; Kamenetzky et al. 2016). The mid-*J* (4 ≤ *J* ≤ 8) CO lines are reported to be linearly correlated with FIR luminosity (e.g., Greve et al. 2014; Lu et al. 2014; Liu et al. 2015; Kamenetzky et al. 2016), which is consistent with the gas heated by the far-ultraviolet (FUV; 6 eV < *h*ν < 13.6 eV) photons from young massive stars (Carilli & Walter 2013). The high-*J* (*J* ≥ 9) CO lines require high temperature and density to be excited, and processes such as X-ray (~1–100 keV) heating from the AGN, mechanical heating from shocks, and cosmic-ray heating are

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
frequently proposed to explain their excitation (e.g., Weiß et al. 2007; Spinoglio et al. 2012; Meijerink & Kristensen 2013; Li et al. 2020b). Thus, the CO SLEDs also serve as probes for discriminating between different gas heating scenarios (Meijerink & Spaans 2005; Meijerink et al. 2007; Spaans 2008), e.g., the photodissociation region (PDR), where the gas is heated by the FUV photons from young massive stars and the AGN-heated X-ray-dominated region (XDR).

Searching for direct evidence of the influence of the AGN activity on the ISM excitation is observationally challenging even in well-studied local/low-$z$ samples. Valentino et al. (2021) find that for a sample of normal $z \sim 1-1.7$ galaxies that host AGNs, marginal correlations have been found between the presence and strength of AGNs and the CO excitation on galactic scales. On nuclear scales of 250 pc, Esposito et al. (2022) report weak correlations between the CO excitation and either the X-ray or FUV flux in a sample of 35 local active galaxies, which suggests that neither the AGNs nor star formation has a dominant effect on the gas excitation in the centers of these galaxies. Studies of CO SLEDs of individual local AGN host galaxies suggest a mix of a PDR component that contributes to the low-$z$ mid-$J$ CO fluxes and an XDR that dominates the high-$J$ CO emission (e.g., van der Werf et al. 2010; Spinoglio et al. 2012; Pozzi et al. 2017; Mingozzi et al. 2018). Other mechanisms, such as mechanical heating by shocks, dense PDRs, and cosmic-ray heating, are also reported to dominate the high-$J$ CO excitation in some local AGN hosts and starburst galaxies (e.g., van der Werf et al. 2010; Spinoglio et al. 2012; Pellegrini et al. 2013).

Luminous quasars at high $z$ represent the most luminous AGNs across redshifts. They are reported to be responsible for their high CO excitation at high $J$. The lensed quasar APM 08279 at $z=3.9$ is currently the highest excited CO SLED ever published across redshifts. Its CO emission is detected within the central $\sim 550$ pc scale (Riechers et al. 2009), and the high CO excitation is best explained by the X-ray heating from the AGN (Weiß et al. 2007; Bradford et al. 2011). Another example is the Cloverleaf quasar at $z=2.6$, for which the CO SLED is best fitted with a combination of a PDR and an XDR component (Uzgil et al. 2016).

Observations of quasars at $z \sim 6$ from optical to (sub) millimeter wavelengths reveal that $\sim 1/3$ of them reside in extremely active environments with rapid SMBH accretion at rates of $M_{BH} > 10 M_{\odot} \text{yr}^{-1}$ and intense star formation at rates up to a few times $1000 M_{\odot} \text{yr}^{-1}$ in the host galaxy (e.g., Wang et al. 2008; Carilli & Walter 2013; De Rosa et al. 2014; Venemans et al. 2020). In such systems, both the star formation and the SMBH activity lead to the injection of vast amounts of energy into the ISM of the host galaxy, making its ISM an ideal target in which to study the early coevolution between the SMBHs and their host galaxies. Well-sampled CO SLEDs (with at least four transitions detected) have been obtained for five quasars at $z \sim 6$ (Bertoldi et al. 2003; Walter et al. 2003; Riechers et al. 2009; Gallerani et al. 2014; Carniani et al. 2019; Wang et al. 2019; Yang et al. 2019; Li et al. 2020b; Pensabene et al. 2021). Notably, all of them show highly excited CO SLEDs with much brighter high-$J$ CO lines than that of the local starburst galaxies and AGNs, suggesting that the high-$J$ CO emission is likely dominated by an XDR component (Gallerani et al. 2014; Wang et al. 2019; Yang et al. 2019; Li et al. 2020b; Pensabene et al. 2021). However, the shapes of

the CO SLEDs vary from object to object. For example, the CO SLED of J0100+2802 is double-peaked, suggesting more than one gas component (Wang et al. 2019). In contrast, the CO SLED of J2310+1855 reveals one dominant gas component for the mid- and high-$J$ CO transitions (Li et al. 2020b).

In this paper, we use the Northern Extended Millimeter Array (NOEMA) operated by the Institut de Radioastronomie Millimétrique (IRAM) to sample the CO SLEDs of the quasar P215-16 at $z=5.78$ and the quasar J1429+5447 at $z=6.18$, each with at least five transitions ranging from CO(5-4) to CO(10-9). P215-16 is among the optically brightest quasars at $z \sim 6$, first identified in the Pan-STARRS1 survey (Morganson et al. 2012). It is also the source with the brightest 450 and 850 $\mu$m flux densities in our SCUBA 2 High Redshift Bright Quasar Survey (SHERRY) of $5.6 < z < 6.9$ quasars (Li et al. 2020). J1429+5447 is among the most distant radio-loud quasars at $z > 6$ and was discovered in the Canadian-French-Hawaii Quasar Survey (Willott et al. 2010). It is also among the most X-ray-luminous quasars observed by eROSITA (Medvedev et al. 2020; Khorunzhev et al. 2021; Migliori et al. 2023), and the second (sub)millimeter brightest and the $[C II]_{158\mu m}$ brightest radio-loud quasar published at $z > 6$ (Khusanova et al. 2022).

The structure of this paper is as follows. We describe observations and data reduction in Section 2 and observational results in Section 3. In Section 4, we derive the physical conditions of the dust and the gaseous ISM using dust spectral energy distribution models and radiative transfer models, respectively. Discussions of the gas properties and comparison with other galaxy samples are presented in Section 5. We adopt a flat $\Lambda$CDM cosmology model with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$ and $\Omega_M = 0.3$, so that $1''$ corresponds to $5.8$ and $5.6$ kpc at $z=5.78$ and $z=6.18$, respectively.

## 2. Observations and Data Reduction

We obtained NOEMA 3 and 2 mm observations of the quasars P215-16 at $z=5.78$ and J1429+5447 at $z=6.18$, targeting the CO(5-4), CO(6-5), CO(7-6), CO(8-7), CO(9-8), CO(10-9), $OH^+(1_1-0_1)$, $[CII](2-1)$, $H_2O(3_{12}-2_{21})$, and $H_2O(3_{21}-3_{12})$ emission lines, together with the underlying dust continuum emission. The observations were carried out in three separate projects (program ID: W0C3, W18EE, and S20CW), starting on 2012 December 24 and finishing on 2020 December 27, a period where the number of antennas increased from 6 to 11.

We used three separate executions in the compact C/D array configurations to observe the CO(5-4), CO(6-5), CO(8-7), CO(9-8), CO(10-9), $OH^+(1_1-0_1)$, $H_2O(3_{12}-2_{21})$, and $H_2O(3_{21}-3_{12})$ emission lines, as well as the dust continuum emission of the quasar P215-16. Observations of the CO(5-4) and CO(6-5) emission lines were performed with an on-source time of 2.3 hr. We tuned the PolyFix correlator to cover the CO(5-4) line in the lower sideband (LSB) and the CO(6-5) line in the upper sideband (USB), each with a 7.744 GHz bandwidth. The CO(8-7), CO(9-8), and $OH^+(1_1-0_1)$ lines were observed in one frequency setup for 1.9 hr on source, with the CO(8-7) and CO(9-8)/$OH^+(1_1-0_1)$ lines observed in the LSB and USB, respectively. We simultaneously observed the CO(9-8), CO(10-9), $OH^+(1_1-0_1)$, $H_2O(3_{12}-2_{21})$, and $H_2O(3_{21}-3_{12})$ lines with one tuning for 1.7 hr on source, with the CO(9-8) and $OH^+(1_1-0_1)$ covered in the LSB and the CO(10-9), $H_2O(3_{12}-2_{21})$, and $H_2O(3_{21}-3_{12})$ lines covered in
<table><caption>**Table 1**
Measured Emission Line and Continuum Properties</caption><thead><tr><th>Source</th><th>Line</th><th>ν<sub>rest</sub><br>(GHz)</th><th>z</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>S δν<br>(Jy km s<sup>-1</sup>)</th><th>Luminosity<br>(×10<sup>8</sup> L<sub>☉</sub>)</th><th>Beam Size</th><th>Channel Width<br>(km s<sup>-1</sup>)</th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th><th>(5)</th><th>(6)</th><th>(7)</th><th>(8)</th><th>(9)</th></tr></thead><tbody><tr><td rowspan="7">P215-16</td><td>CO(5-4)</td><td>576.268</td><td>5.7831 ± 0.0006</td><td>444 ± 63</td><td>0.90 ± 0.17</td><td>2.4 ± 0.5</td><td>5.774 × 1.73, PA = 5°</td><td>56</td></tr><tr><td>CO(6-5)</td><td>691.473</td><td>5.7837 ± 0.0004</td><td>404 ± 42</td><td>1.01 ± 0.14</td><td>3.3 ± 0.5</td><td>4.882 × 1.45, PA = 4°</td><td>46</td></tr><tr><td>CO(8-7)</td><td>921.800</td><td>5.7834 ± 0.0003</td><td>443 ± 39</td><td>1.54 ± 0.18</td><td>6.7 ± 0.8</td><td>4.445 × 0.92, PA = 8°</td><td>70</td></tr><tr><td>OH<sup>+</sup>(1<sub>1</sub>-0<sub>1</sub>)</td><td>1033.118</td><td>5.7824</td><td>503</td><td>&lt;0.57</td><td>&lt;2.7</td><td>3.14 × 1.09, PA = 5°</td><td>78</td></tr><tr><td>CO(9-8)</td><td>1036.912</td><td>5.7824 ± 0.0004</td><td>503 ± 37</td><td>1.93 ± 0.19</td><td>9.4 ± 0.9</td><td>3.14 × 1.09, PA = 5°</td><td>78</td></tr><tr><td>CO(10-9)</td><td>1151.985</td><td>5.7824</td><td>491 ± 41</td><td>2.00 ± 0.24</td><td>10.8 ± 3.6</td><td>4.41 × 1.09, PA = 6°</td><td>71</td></tr><tr><td>H<sub>2</sub>O(3<sub>12</sub>-2<sub>21</sub>)</td><td>1153.127</td><td>5.7824</td><td>491 ± 41</td><td>0.42 ± 0.18</td><td>2.2 ± 0.9</td><td>4.41 × 1.09, PA = 6°</td><td>71</td></tr><tr><td>H<sub>2</sub>O(3<sub>21</sub>-3<sub>12</sub>)</td><td>1162.912</td><td>5.7824</td><td>491 ± 41</td><td>1.32 ± 0.18</td><td>7.2 ± 2.0</td><td>4.41 × 1.09, PA = 6°</td><td>71</td></tr><tr><td rowspan="7">J1429+5447</td><td>CO(5-4)</td><td>576.268</td><td>6.1824 ± 0.0005</td><td>267 ± 50</td><td>0.39 ± 0.09</td><td>1.2 ± 0.3</td><td>6.32 × 5.11, PA = 46°</td><td>75</td></tr><tr><td>CO(6-5)</td><td>691.473</td><td>6.1841 ± 0.0007</td><td>372 ± 70</td><td>0.53 ± 0.13</td><td>1.9 ± 0.5</td><td>5.33 × 4.33, PA = 48°</td><td>62</td></tr><tr><td>CO(7-6)</td><td>806.652</td><td>6.1841</td><td>372</td><td>0.48 ± 0.10</td><td>2.0 ± 0.4</td><td>2.45 × 2.14, PA = 38°</td><td>62</td></tr><tr><td>[C I](2-1)</td><td>809.342</td><td>6.1841</td><td>372</td><td>&lt;0.31</td><td>&lt;1.3</td><td>2.45 × 2.14, PA = 38°</td><td>62</td></tr><tr><td>OH<sup>+</sup>(1<sub>1</sub>-0<sub>1</sub>)</td><td>1033.118</td><td>6.1860</td><td>271</td><td>&lt;0.21</td><td>&lt;1.1</td><td>3.85 × 2.71, PA = 51°</td><td>50</td></tr><tr><td>CO(9-8)</td><td>1036.912</td><td>6.1860 ± 0.0012</td><td>271 ± 69</td><td>0.29 ± 0.07</td><td>1.3 ± 0.4</td><td>3.85 × 2.71, PA = 51°</td><td>50</td></tr><tr><td>CO(10-9)</td><td>1151.985</td><td>6.1860</td><td>313 ± 119</td><td>0.21 ± 0.06</td><td>1.2 ± 0.3</td><td>3.49 × 2.44, PA = 51°</td><td>45</td></tr><tr><td>H<sub>2</sub>O(3<sub>12</sub>-2<sub>21</sub>)</td><td>1153.127</td><td>6.1860</td><td>313</td><td>&lt;0.18</td><td>&lt;1.1</td><td>3.49 × 2.44, PA = 51°</td><td>45</td></tr><tr><td>H<sub>2</sub>O(3<sub>21</sub>-3<sub>12</sub>)</td><td>1162.912</td><td>6.1860</td><td>313 ± 119</td><td>0.15 ± 0.06</td><td>0.9 ± 0.3</td><td>3.49 × 2.44, PA = 51°</td><td>45</td></tr></tbody></table>

Dust continuum properties

P215-16

J1429+5447

<table><thead><tr><th>ν<sub>cont</sub><br>(GHz)<br>(10)</th><th>Flux density<br>(mJy)<br>(11)</th><th>rms<br>(μJy beam<sup>-1</sup>)<br>(12)</th><th>Ref<br>(13)</th><th>ν<sub>cont</sub><br>(GHz)<br>(14)</th><th>Flux density<br>(mJy)<br>(15)</th><th>rms<br>(μJy beam<sup>-1</sup>)<br>(16)</th><th>Ref<br>(17)</th></tr></thead><tbody><tr><td>83.2</td><td>0.21 ± 0.02</td><td>21</td><td>[1]</td><td>82.2</td><td>0.15 ± 0.02</td><td>19</td><td>[1]</td></tr><tr><td>98.6</td><td>0.27 ± 0.02</td><td>20</td><td>[1]</td><td>98.4</td><td>0.09 ± 0.02</td><td>20</td><td>[1]</td></tr><tr><td>138.3</td><td>1.38 ± 0.03</td><td>29</td><td>[1]</td><td>112.4</td><td>0.34 ± 0.05</td><td>52</td><td>[1]</td></tr><tr><td>153.2</td><td>1.84 ± 0.03</td><td>33</td><td>[1]</td><td>144.2</td><td>0.32 ± 0.02</td><td>23</td><td>[1]</td></tr><tr><td>168.4</td><td>2.37 ± 0.05</td><td>48</td><td>[1]</td><td>159.8</td><td>0.41 ± 0.04</td><td>35</td><td>[1]</td></tr><tr><td>352.7</td><td>16.85 ± 1.10</td><td>...</td><td>[2]</td><td>1.4</td><td>2.68 ± 0.16</td><td>...</td><td>[3]</td></tr><tr><td rowspan="3">666.2</td><td rowspan="3">26.93 ± 7.78</td><td rowspan="3">...</td><td rowspan="3">[2]</td><td>3</td><td>1.81 ± 0.11</td><td>...</td><td>[4]</td></tr><tr><td>32</td><td>0.26 ± 0.02</td><td>...</td><td>[5]</td></tr><tr><td>250</td><td>3.05 ± 0.11</td><td>...</td><td>[6]</td></tr></tbody></table>

**Note.** Columns (1), (2): source name and line ID. Columns (3)–(7): line rest frequency, redshift, line width in FWHM, flux, and luminosity. Columns (8), (9): beam size and channel width. Columns (10)–(13) and (14)–(17): continuum frequency/wavelength, flux, rms noise, and references for P215–16 and J1429+5447, respectively. References: [1] this work, [2] Li et al. (2020), [3] Villarreal Hernández & Andernach (2018), [4] Becker et al. (1995), [5] Wang et al. (2011), [6] Khusanova et al. (2022).

the USB. The flux calibrator was MWC 349 and the phase/amplitude calibrators were 1334-127, 1437-153, 1406-076, and 1435-218.

As for the quasar J1429+5447, we used three separate executions to cover the CO(5-4), CO(6-5), CO(7-6), [C I] (2-1), CO(9-8), OH<sup>+</sup>(1<sub>1</sub>-0<sub>1</sub>), CO(10-9), H<sub>2</sub>O(3<sub>12</sub>-2<sub>21</sub>), and H<sub>2</sub>O(3<sub>21</sub>-3<sub>12</sub>) emission lines and the underlying continuum emission. The CO(5-4) and the CO(6-5) lines were covered in the LSB and USB, respectively, with one frequency setup, for an on-source time of 3.2 hr. The CO(9-8), OH<sup>+</sup>(1<sub>1</sub>-0<sub>1</sub>), CO(10-9), H<sub>2</sub>O(3<sub>12</sub>-2<sub>21</sub>), and H<sub>2</sub>O(3<sub>21</sub>-3<sub>12</sub>) emission lines were observed simultaneously in one tuning in the 2 mm band for 2.1 hr on source. Observations of the CO(7-6) and [C I] (2-1) lines were conducted in one frequency setup for 2.2 hr on source. We observed the continuum emission using all line-free channels in the USBs and LSBs in each observation. For all the observations, 1418+546 was used as phase/amplitude calibrator. The fluxes were calibrated with MWC 349, 2010+723, and 3C 273. The typical calibration uncertainty is <10% in Band 1 and <15% in Band 2.

We reduced the data using the CLIC and MAPPING packages that are part of the Grenoble Image and Line Data Analysis System (GILDAS) software (Guilloteau & Lucas 2000). We obtained the continuum uv tables by averaging all the line-free channels using the UV_AVERAGE task in MAPPING. To generate the continuum-subtracted spectral line uv tables, we adopted the continuum uv table obtained above as a continuum model, and used the UV_SUBTRACT task to subtract it from the original uv tables. The uv tables of the spectral lines and the continuum were cleaned and imaged with natural weighting to maximize the signal-to-noise ratio (S/N) for a point source. The host galaxies are spatially unresolved in beam sizes of ~3"-6" × 1"-2" for P215-16 and 2"-7" × 2"-5" for J1429+5447. The spectra were finally binned to 46-78 km s⁻¹ for the quasar P215-16 and 45-75 km s⁻¹ for the quasar J1429+5447. The rms noise for the continuum emission in all bands is 0.02-0.05 mJy beam⁻¹ for both sources. Detailed information on the targeted lines and derived parameters can be found in Table 1.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Spectra of the CO, H₂O, OH⁺(1₁-0₁), and [C I](2-1) lines for P215-16 (left) and J1429+5447 (right). The continuum emission has been subtracted for each of the spectra. The data are shown as yellow histograms. Gaussian profiles fitted to the line emission are shown as solid or dashed lines. The expected positions of the undetected lines are indicated as dashed blue vertical lines. In the spectra where more than one line is detected, the solid black lines represent the sum of different Gaussian components.</figcaption></figure>

### 3. Results

We detect the CO(5–4), CO(6–5), CO(8–7), CO(9–8), CO(10–9), and the H₂O(3₂₁–3₁₂) lines from the host galaxy of P215–16 at S/N > 5, and report a nondetection of the OH⁺(1₁–0₁) line (Figure 1). A tentative detection (~2.3σ) of the H₂O(3₁₂–2₂₁) line blended with the CO(10–9) line is also

obtained for this object. The continuum emission is detected at S/N > 10 in all the observed bands (Table 1). We extract the spectra from the peak pixels and fit the spectra of CO(5–4), CO(6–5), CO(8–7), and CO(9–8) lines separately, each with a Gaussian profile. The redshifts and line widths derived for different CO lines are consistent with each other within the
uncertainties. The CO(10–9), H₂O(3₁₂–2₁₂), and H₂O(3₂₁–3₁₂) lines are covered in the same sideband. We fit these three lines simultaneously, using a Gaussian model, with the line centers fixed to the redshift derived from the CO(9–8) line and an identical line width for all of the three spectral lines (assuming that these three lines originate from the same gas component). The spectra and the Gaussian fitting results are shown in Figure 1. Adopting the line width of the CO(9–8) line, we obtain a velocity-integrated map centered at the OH⁺(1₁–0₁) line frequency and measure the 3σ detection upper limit for the OH⁺(1₁–0₁) line. The measurements for both spectral lines and dust continuum are presented in Table 1.

For the quasar J1429+5447, we detect the CO(5–4), CO(6–5), CO(7–6), and CO(9–8) lines at S/N > 4, and tentatively detect the CO(10–9) and H₂O(3₂₁–3₁₂) lines at the 3.5σ and 2.5σ levels, respectively. The [C I](2–1), H₂O(3₁₂–2₁₂), and OH⁺(1₁–0₁) lines are undetected (Figure 1). The CO lines detected in J1429+5447 have lower S/N than those detected in P215–16. The S/N of the CO(7–6) line spectrum is low. We fix the Gaussian line width and redshift to that of the CO(6–5) line and fit the line flux. The resulting line widths, redshifts, fluxes, and continuum fluxes are listed in Table 1. The 3σ upper limits of the [C I](2–1), OH⁺(1₁–0₁), and H₂O(3₁₂–2₁₂) lines are estimated based on the rms of the moment 0 maps, which are integrated over channels expected for these lines assuming the same line widths and redshifts as those of their neighboring CO lines. The derived redshifts and line widths of the CO lines are consistent with each other within the uncertainties. These are also consistent with the fitting results of the W component of the CO(2–1) line (the western peak that coincides with the optical quasar) reported in Wang et al. (2011). The continuum emissions are all detected at a >4σ level. More details about the measured dust and spectral line properties are listed in Table 1.

# 4. Dust and Gas Properties

## 4.1. Dust SED Fitting

To derive the dust properties of the two quasars, we analyze their dust spectral energy distributions (SEDs). An optically thick dust model has been suggested to explain the dust emission in a few high-z starburst systems, especially at rest-frame wavelengths ≲200 µm (e.g., Riechers et al. 2013; Wang et al. 2019). Since the quasars in this work are all far-infrared bright sources, it is reasonable to consider that the dust can be optically thick close to the peak of the dust SED. We thus employ the general opacity dust formula published in da Cunha et al. (2021) to fit the dust SEDs of the quasars in this work. In the general opacity regime, the dust emission can be described by

$$S_{\nu} \propto [1 - \exp(-\tau_{\nu})](B_{\nu}[T_{\text{dust}}(z)] - B_{\nu}[T_{\text{CMB}}(z)]), \quad (1)$$

taking into account the effect of the cosmic microwave background (CMB). Here $B_{\nu}[T_{\text{dust}}(z)]$ or $B_{\nu}[T_{\text{CMB}}(z)]$ is the Planck function at the temperature of the dust or CMB at the quasar’s redshift, and the dust optical depth $\tau_{\nu}$ is proportional to the dust surface density $\Sigma_{\text{dust}}$ via

$$\tau_{\nu} = \kappa_{\nu} \Sigma_{\text{dust}}, \quad (2)$$

where $\kappa_{\nu}$ is the dust opacity. Its dependence on frequency is

$$\kappa_{\nu} = \kappa_{0} \left( \frac{\nu}{\nu_{0}} \right)^{\beta}, \quad (3)$$

where $\beta$ is the dust emissivity index and $\kappa_0$ is the emissivity of dust grains per unit mass at frequency $\nu_0$. If we assume a simple spherical geometry, then the dust mass is

$$M_{\text{dust}} = 4\pi R^2 \Sigma_{\text{dust}}, \quad (4)$$

where R is the dust radius of a galaxy. In da Cunha et al. (2021), an additional parameter $\lambda_{\text{thick}}$ is introduced, which is defined when the optical depth reaches 1,

$$\begin{aligned} \tau_{\lambda_{\text{thick}}} &= 1 \implies \lambda_{\text{thick}} \\ &= \lambda_{0} \left( \kappa_{0} \frac{M_{\text{dust}}}{4\pi R^{2}} \right)^{1/\beta}. \end{aligned} \quad (5)$$

Accordingly, the dust emission can be rewritten as

$$S_{\nu} \propto \left[ 1 - \exp\left(-\left(\frac{\nu\lambda_{\text{thick}}}{c}\right)^{\beta}\right) \right] \times (B_{\nu}[T_{\text{dust}}(z)] - B_{\nu}[T_{\text{CMB}}(z)]), \quad (6)$$

where c is the speed of light.

Complemented with SCUBA 2 observations of the dust continuum at 450 µm (666.2 GHz) and 850 µm (352.7 GHz) from Li et al. (2020), we use the general opacity model (Equation (6)) to fit the observed dust SED of P215–16, normalized to the continuum flux density at 666.2 GHz. Under the general opacity assumption, the dust SED is determined by three parameters, namely $T_{\text{dust}}$, $\beta$, and $\lambda_{\text{thick}}$. We consider three different dust models. To break the degeneracy between different parameters in dust SED fitting, we fix $\lambda_{\text{thick}}$ to (1) the average value for submillimeter galaxies (SMGs; $\lambda_{\text{thick}} = 100$ µm; da Cunha et al. 2021; Model 1), and (2) the average value for QSOs at $z \sim 6$ ($\lambda_{\text{thick}} = 264$ µm; derived from Gilli et al. 2022; Model 2), while enabling $T_{\text{dust}}$ and $\beta$ to vary in these two models. In a third model, we enable all three parameters $\lambda_{\text{thick}}$, $T_{\text{dust}}$, and $\beta$ to vary (Model 3). To explore the posterior distributions of the parameters, we use the Markov Chain Monte Carlo (MCMC) Python package emcee to fit the data to the model (Foreman-Mackey et al. 2013), adopting the best-fit values obtained for the least-squares estimation as educated initial guesses of the parameters.

The inclusion of the high-frequency SCUBA 2 observations is critical in constraining the dust SED parameters for P215–16. The least-squares estimation finds best-fit values of $T_{\text{dust}} = 48$ K and $\beta = 2.05$ for Model 1, $T_{\text{dust}} = 111$ K and $\beta = 2.62$ for Model 2, and $T_{\text{dust}} = 80$ K, $\beta = 2.33$, and $\lambda_{\text{thick}} = 199$ µm for Model 3. In some studies, $T_{\text{dust}}$, $\beta$, and $\tau_{1900 \text{ GHz}}$ are adopted to describe the general opacity dust SED model (e.g., Walter et al. 2022; Decarli et al. 2023), where $\tau_{1900 \text{ GHz}}$ is the dust optical depth at 1900 GHz. The derived $\lambda_{\text{thick}}$ values for our model translate into $\tau_{1900 \text{ GHz}}$ of 0.39, 3.84, and 1.71 for Models 1, 2, and 3, respectively, which is consistent with values derived for other quasars at $z \sim 6$ (Walter et al. 2022; Decarli et al. 2023). By using the MCMC method, we obtain fitting results of $T_{\text{dust}} = 50^{+6}_{-5}$ K and $\beta = 1.99^{+0.22}_{-0.22}$, $T_{\text{dust}} = 116^{+18}_{-11}$ K and $\beta = 2.53^{+0.37}_{-0.34}$, and $T_{\text{dust}} = 108^{+47}_{-33}$ K, $\beta = 2.49^{+0.57}_{-0.44}$, and $\lambda_{\text{thick}} = 252^{+57}_{-65}$ µm for
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) Dust SED fitting results for P215−16. The data are shown as black squares. The top and bottom axes show the observed and rest continuum frequencies. The dust SED is fitted with three general opacity models, and the solid lines represent the models with best-fit values (see Section 4.1 for further details). (b)–(d) Posterior distributions of the model parameters for Model 3, 1, and 2, with the same color as the best-fit models in (a). The 2D contours show the [1, 2]σ confidence intervals in the marginalized distributions.</figcaption></figure>

Models 1, 2, and 3, respectively (see Figure 2). The derived dust temperature and emissivity index are in agreement with results obtained for other quasars at $z \sim 6$ when adopting an optically thick assumption (e.g., Wang et al. 2019; Pensabene et al. 2022), as well as results obtained from 3D radiative transfer calculations (Di Mascia et al. 2021). The resulting infrared (IR, 8 µm–1000 µm) luminosities ($L_{IR}$) are $1.80_{-0.20}^{+0.20} \times 10^{13} L_{\odot}$, $4.61_{-0.25}^{+0.25} \times 10^{13} L_{\odot}$, and $2.70_{-0.40}^{+0.40} \times 10^{13} L_{\odot}$ for the three models, respectively. If we adopt the

previously widely used optically thin modified blackbody dust model for $z \sim 6$ quasars with parameters of $T_{\text{dust}} = 47$ K and $\beta = 1.6$ (Beelen et al. 2006), the derived infrared luminosity is $1.9 \times 10^{13} L_{\odot}$, which is close to the result of Model 1.

J1429+5447 is a radio-loud quasar, so we also include a power-law radio emission component $f_{\nu} \propto \nu^{\alpha}$, with a power-law index of $\alpha$ to fit the data in addition to the thermal dust SED model (Equation (6)). Besides the continuum emission detected in this paper, we include published continuum data at
1.4 and 3 GHz from the Very Large Array Sky Survey (VLASS) and Faint Images of the Radio Sky at Twenty Centimeters (FIRST; Becker et al. 1995; Villarreal Hernández & Andernach 2018), and at 32 and 250 GHz obtained using the Karl Jansky Very Large Array (VLA) and NOEMA (Wang et al. 2011; Khusanova et al. 2022), respectively. The observational data cannot constrain all the model parameters for the radio power-law + general opacity model. We thus fix $\lambda_{\text{thick}} = 100$ µm (typical of SMGs) and explore two cases with different assumptions for $T_{\text{dust}}$: namely, 50 K (Model 1) and 100 K (Model 2), representing typical “cold” and “warm” dust, respectively.<fnref n="19" />

The least-squares estimation suggests best-fit values of $\alpha = -0.76$ and $\beta = 2.95$ for Model 1, and $\alpha = -0.76$ and $\beta = 2.66$ for Model 2. The resulting values of $\tau_{1900\,\text{GHz}}$ are 0.26 and 0.30 for Models 1 and 2, respectively, which are consistent with values derived for other quasars at $z \sim 6$ (Walter et al. 2022; Decarli et al. 2023). The MCMC method finds $\beta = 3.02_{-0.53}^{+0.50}$ and $\alpha = -0.76_{-0.05}^{+0.04}$, and $\beta = 2.74_{-0.55}^{+0.58}$ and $\alpha = -0.76_{-0.05}^{+0.04}$, for Models 1 and 2, respectively (Figure 3). The derived infrared luminosities are $6.89_{-0.18}^{+0.18} \times 10^{12} L_{\odot}$ and $3.68_{-0.42}^{+0.42} \times 10^{13} L_{\odot}$ for Models 1 and 2, respectively. If we adopt the widely used modified blackbody dust SED model with parameters $T_{\text{dust}} = 47$ K and $\beta = 1.6$ instead, the infrared luminosity is $9.5 \times 10^{12} L_{\odot}$.

The derived dust properties of P215-16 and J1429+5447 are presented in Table 2. Future high-frequency continuum observations, possibly sampling the peak of the dust SED and/or the Wien tail, will be critical to further constrain the dust temperature as well as the IR luminosity.

## 4.2. ISM Properties

CO SLEDs are sensitive probes of the physical conditions (e.g., temperature, density, and radiation field) of the molecular gas. To derive the physical properties of the molecular gas in the host galaxies of these two quasars, we employ the CO SLED grid models published in Pensabene et al. (2021). These models consider two different gas heating mechanisms, namely a PDR heated by the FUV photons from young massive stars and an XDR heated by the X-ray photons from the AGNs. In the PDR model, the shape of the CO SLED is described by two physical parameters: the gas volume density $n(H_2)$ and the FUV radiation field strength $G_0$ in Habing field units (1 Habing field corresponds to $1.6 \times 10^{-3} \text{ erg s}^{-1} \text{ cm}^{-2}$; Habing 1968). The CO SLED in the XDR model is determined by the gas volume density $n(H_2)$ and the strength of the X-ray radiation field $F_x$ in units of erg s$^{-1}$ cm$^{-2}$. To fully sample the molecular region, which has a typical column density of $>2 \times 10^{22} \text{ cm}^{-2}$ (McKee & Ostriker 2007), we fix the column density of molecular hydrogen to $10^{23} \text{ cm}^{-2}$ for both sets of the PDR and XDR models. The adopted column density is also consistent with the values derived/adopted for other $z \sim 6$ quasars (e.g., Meyer et al. 2022; Decarli et al. 2023). Ranges of $n(H_2)$, $G_0$, and $F_x$ for typical molecular clouds explored in the model grids are $10^2-10^6 \text{ cm}^{-3}$, $10^1-10^6$, and $10^{-2}-10^2 \text{ erg s}^{-1} \text{ cm}^{-2}$ (Pensabene et al. 2021). We consider two cases for the CO SLED fitting: (1) a single PDR component and (2) an XDR component. The fitting procedure for each case is as follows. First, we use the least-squares estimation for a first-order

estimation of the best-fit values of the parameters. Second, we employ the MCMC package emcee to sample the posterior distribution of the parameters while using the best-fit values derived above as educated initial guesses for the parameters (Foreman-Mackey et al. 2013).

We show the fitting results in Figure 4. The CO SLED of P215-16 cannot be fitted with any of our PDR models with $n(H_2)$ and $G_0$ in the ranges $10^2-10^6 \text{ cm}^{-3}$ and $10^1-10^6$ (Figure 4(a)), which are typical ranges for the PDR regions in star-forming galaxies, starbursts, and AGNs (e.g., Gallerani et al. 2014; Pereira-Santaella et al. 2014; Pensabene et al. 2021; Esposito et al. 2022). To further explore whether the data can be fitted with an extreme PDR model, we extend the range of $n(H_2)$ and $G_0$ parameters of Pensabene et al. (2021) to $10^2-10^8 \text{ cm}^{-3}$ and $10^1-10^8$, respectively, while using the same set of assumptions in CLOUDY. The results suggest that the data can be fitted with an extreme PDR of $n(H_2) = 10^{7.1} \text{ cm}^{-3}$ and $G_0 = 10^{7.8}$ (Figure 4(b)). The CO SLED of P215-16 can also be fitted with an XDR model. The least-squares approximation finds a dense XDR component of $n(H_2) = 10^{4.3} \text{ cm}^{-3}$ illuminated by a radiation field of $F_x = 10^{1.7} \text{ erg s}^{-1} \text{ cm}^{-2}$. The MCMC analysis finds $n(H_2) = 10^{4.7 \pm 0.8}_{-0.6}$ and $F_x = 10^{0.87 \pm 0.6}_{-0.4}$ erg s$^{-1} \text{ cm}^{-2}$ for the XDR model, where values and uncertainties are estimated based on the 50th, 16th, and 84th percentiles of the samples in the marginalized distributions of the parameters. It is also possible that the CO SLED of P215-16 is contributed by both PDR and XDR components. However, the current CO SLED measurements are insufficient to constrain the parameters for a PDR + XDR model. Future observations of the high-$J$ ($J \ge 11$) CO lines will be critical to discriminate between different models.

For J1429+5447, the combination of available [C II]$_{158\,\mu m}$ data with the information on the nondetection of the [C I](2–1) line from this work enables us to constrain the physical conditions of the atomic gas. The [C II]$_{158\,\mu m}$/$[\text{C I}](2–1)$ luminosity ratio of $>23$ argues for a PDR origin of the line excitation, as the XDR scenario usually results in a lower line ratio of $\lesssim 10$ (Pensabene et al. 2021). Predictions of the [C II]$_{158\,\mu m}$/$[\text{C I}](2–1)$ ratio in the PDR model are shown in Figure 5 (Pensabene et al. 2021). The [C II]$_{158\,\mu m}$/$[\text{C I}](2–1)$, [C II]$_{158\,\mu m}/L_{\text{IR}}$, and [C I](2–1)/$L_{\text{IR}}$ ratios indicate a PDR component with parameters in the region of $n(H_2) = 10^{3.6} \text{ cm}^{-3}$ and $G_0 = 10^{3.7-4.7}$, adopting an $L_{\text{IR}}$ value of $6.89_{-0.18}^{+0.18} \times 10^{12} L_{\odot}-3.68_{-0.42}^{+0.42} \times 10^{13} L_{\odot}$. Detections of the CO emission lines complemented with previous CO(2–1) observations (Wang et al. 2011) enable us to constrain, for the first time, the physical conditions of the molecular ISM of a radio-loud quasar at $z > 6$.

We analyze the CO SLED of J1429+5447 following a similar procedure to that of P215-16. A single “dense” PDR model illuminated by an intense FUV radiation field is able to reproduce the observational data, with parameters of $n(H_2) = 10^{6.0} \text{ cm}^{-3}$ and $G_0 = 10^{5.4}$. A single XDR model is also able to reproduce the CO SLED with best-fit parameters of $n(H_2) = 10^{3.7} \text{ cm}^{-3}$ and $F_x = 10^{1.5} \text{ erg s}^{-1} \text{ cm}^{-2}$. Using the gas conditions derived from [C II]$_{158\,\mu m}$, [C I](2–1) and $L_{\text{IR}}$ as a constraint on the PDR component, we also try to fit the CO SLED of J1429+5447 with a PDR + XDR model. This leads to a best-fit result of $n(H_2) = 10^{4.3} \text{ cm}^{-3}$ and $F_x = 10^{2.0} \text{ erg s}^{-1} \text{ cm}^{-2}$ for the XDR component, and $n(H_2) = 10^{5.1} \text{ cm}^{-3}$ and $G_0 = 10^{4.2}$ for the PDR component. The PDR component contributes to 96% of the total molecular gas mass

<footnote n="19">¹⁹ We also considered fixing $\lambda_{\text{thick}} = 264$ µm (typical of QSOs at $z \sim 6$), but the models could not fit the SED of J1429+5447 at 250 GHz.</footnote>
<figure><img src="image_6.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) Model 1</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) Model 2</figcaption></figure>

**Figure 3.** (a) Dust SED fitting results for J1429+5447. The data are shown as black squares. The top and bottom axes show the observed and rest continuum frequencies. The dust SED is fitted with the general opacity dust SED model + radio power-law emission (see Section 4.1 for further details). The solid/dashed blue (red) lines represent the best-fit general opacity/radio power-law models for Model 1 (Model 2). (b), (c) Posterior distributions of the model parameters for J1429+5447, respectively. The color scheme for different models is the same as in panel (a). The 2D contours show the [1, 2]σ confidence intervals in the marginalized distributions.

in this case. The MCMC method finds best-fit values in the ranges $n(\text{H}_2) = 10^{5.8^{+0.7}_{-0.3}} \text{ cm}^{-3}$ and $G_0 = 10^{5.6^{+0.7}_{-0.7}}$ (PDR model), $n(\text{H}_2) = 10^{3.8^{+1.1}_{-0.3}} \text{ cm}^{-3}$ and $F_\chi = 10^{0.33^{+1.4}_{-0.7}} \text{ erg s}^{-1} \text{ cm}^{-2}$ (XDR model), and $n(\text{H}_2) = 10^{5.1^{+1.1}_{-0.7}} \text{ cm}^{-3}$ and $G_0 = 10^{4.2^{+0.5}_{-0.5}}$, and $n(\text{H}_2) = 10^{4.4^{+0.9}_{-0.7}} \text{ cm}^{-3}$ and $F_\chi = 10^{1.3^{+1.3}_{-1.3}} \text{ erg s}^{-1} \text{ cm}^{-2}$ (PDR + XDR model). Analysis of the gas properties using either the

[C II]$_{158\mu m}$, [C I](2–1) and IR ratios or the CO SLED suggests that the [C II]$_{158\mu m}$, [C I](2–1), CO, and IR emission can be explained with one PDR component. However, we are not able to rule out a single XDR model or a composite PDR + XDR model, which reproduces the observed CO SLED. In the case of either the XDR or PDR + XDR model, a “moderate density”
<table><caption>**Table 2**

Dust SED and CO SLED Modeling Results</caption><thead><tr><th>Source</th><th>Model</th><th>LS Best Fit</th><th>MCMC</th></tr></thead><tbody><tr><td rowspan="3">Dust SED<br>P215-16</td><td>Model(λ<sub>thick-fixed</sub>, T<sub>dust-fixed</sub>)</td><td></td><td></td></tr><tr><td>Model 1 (100 µm, /)</td><td>T<sub>dust</sub> = 48 K, β = 2.05</td><td>T<sub>dust</sub> = 50<sup>+6</sup><sub>-5</sub> K, β = 1.99<sup>+0.22</sup><sub>-0.22</sub>, L<sub>IR</sub> = 1.80<sup>+0.20</sup><sub>-0.20</sub> × 10<sup>13</sup> L<sub>☉</sub></td></tr><tr><td>Model 2 (264 µm, /)</td><td>T<sub>dust</sub> = 111 K, β = 2.62</td><td>T<sub>dust</sub> = 116<sup>+18</sup><sub>-11</sub> K, β = 2.53<sup>+0.37</sup><sub>-0.34</sub>, L<sub>IR</sub> = 4.61<sup>+0.25</sup><sub>-0.25</sub> × 10<sup>13</sup> L<sub>☉</sub></td></tr><tr><td></td><td>Model 3 ( / , /)</td><td>T<sub>dust</sub> = 80 K, β = 2.33,<br>λ<sub>thick</sub> = 199 µm</td><td>T<sub>dust</sub> = 108<sup>+47</sup><sub>-33</sub> K, β = 2.49<sup>+0.57</sup><sub>-0.44</sub>, λ<sub>thick</sub> = 252<sup>+57</sup><sub>-65</sub> µm, L<sub>IR</sub> =<br>2.70<sup>+0.40</sup><sub>-0.40</sub> × 10<sup>13</sup> L<sub>☉</sub></td></tr><tr><td rowspan="2">J1429+5447</td><td>Model 1 (100 µm, 50 K)</td><td>α = -0.76, β = 2.95</td><td>β = 3.02<sup>+0.50</sup><sub>-0.53</sub>, α = -0.76<sup>+0.04</sup><sub>-0.05</sub>, L<sub>IR</sub> = 6.89<sup>+0.18</sup><sub>-0.18</sub> × 10<sup>12</sup> L<sub>☉</sub></td></tr><tr><td>Model 2 (100 µm, 100 K)</td><td>α = -0.76, β = 2.66</td><td>β = 2.74<sup>+0.58</sup><sub>-0.55</sub>, α = -0.76<sup>+0.04</sup><sub>-0.05</sub>, L<sub>IR</sub> = 3.68<sup>+0.42</sup><sub>-0.42</sub> × 10<sup>13</sup> L<sub>☉</sub></td></tr><tr><td rowspan="5">CO SLED<br>P215-16</td><td>PDR</td><td>n(H<sub>2</sub>) = 10<sup>6.0</sup>, G<sub>0</sub> = 10<sup>6.0</sup></td><td>...</td></tr><tr><td>Extreme PDR</td><td>n(H<sub>2</sub>) = 10<sup>7.1</sup>, G<sub>0</sub> = 10<sup>7.8</sup></td><td>...</td></tr><tr><td>XDR</td><td>n(H<sub>2</sub>) = 10<sup>4.3</sup>, F<sub>x</sub> = 10<sup>1.1</sup></td><td>n(H<sub>2</sub>) = 10<sup>4.7+0.8-0.6</sup>, F<sub>x</sub> = 10<sup>0.87+0.6-0.4</sup></td></tr><tr><td>PDR</td><td>n(H<sub>2</sub>) = 10<sup>6.0</sup>, G<sub>0</sub> = 10<sup>5.4</sup></td><td>n(H<sub>2</sub>) = 10<sup>5.8+0.7-0.3</sup>, G<sub>0</sub> = 10<sup>5.6+0.4-0.7</sup></td></tr><tr><td>XDR</td><td>n(H<sub>2</sub>) = 10<sup>3.7</sup>, F<sub>x</sub> = 10<sup>1.5</sup></td><td>n(H<sub>2</sub>) = 10<sup>3.8+1.1-0.3</sup>, F<sub>x</sub> = 10<sup>0.33+1.4-0.7</sup></td></tr><tr><td rowspan="2">J1429+5447</td><td>PDR + XDR</td><td>PDR: n(H<sub>2</sub>) = 10<sup>5.1</sup>, G<sub>0</sub> = 10<sup>4.2</sup><br>XDR: n(H<sub>2</sub>) = 10<sup>4.3</sup>, F<sub>x</sub> = 10<sup>2.0</sup></td><td>n(H<sub>2</sub>) = 10<sup>5.1+1.1</sup>, G<sub>0</sub> = 10<sup>4.2+0.4-0.5</sup><br>n(H<sub>2</sub>) = 10<sup>4.4+0.9-0.7</sup>, F<sub>x</sub> = 10<sup>1.3+1.4-1.3</sup></td></tr></tbody></table>

**Note.** “LS Best Fit” represents the best-fit result of the least-squares method. “MCMC” represents the fitting result of the MCMC method. Model(λ<sub>thick-fixed</sub>, T<sub>dust-fixed</sub>) represents the model where λ<sub>thick</sub> or T<sub>dust</sub> is fixed to λ<sub>thick-fixed</sub> or T<sub>dust-fixed</sub>. For the CO SLED modeling, n(H<sub>2</sub>), G<sub>0</sub>, and F<sub>x</sub> are in units of cm<sup>-3</sup>, Habing field, and erg s<sup>-1</sup> cm<sup>-2</sup>.

XDR component illuminated by an intense radiation field is suggested.

Results for the gas properties derived from CO SLED modeling are shown in Table 2. In short, both the extreme PDR and the XDR models could explain the observed CO SLEDs (up to J = 10) of P215-16. However, we also notice that the extreme PDR model for P215-16 requires extremely high gas density and FUV radiation field intensity, which are rarely seen in starburst systems and may suggest an AGN contribution. For J1429+5447, the CO SLED could be fitted either with a PDR, an XDR, or a PDR + XDR model. More discussions about the excitation mechanism of the molecular gas can be found in Section 5.

In addition, we robustly detect the H<sub>2</sub>O(3<sub>21</sub>-3<sub>12</sub>) line in the quasar P215-16, and tentatively detect the H<sub>2</sub>O(3<sub>12</sub>-2<sub>21</sub>) in P215-16 and H<sub>2</sub>O(3<sub>12</sub>-2<sub>21</sub>) and H<sub>2</sub>O(3<sub>21</sub>-3<sub>12</sub>) in J1429+5447. The H<sub>2</sub>O molecule, as the third most abundant species after H<sub>2</sub> and CO in the molecular ISM phase, traces obscured star formation in dense molecular clouds. Linear relations were reported between the mid/high-J H<sub>2</sub>O and the IR luminosities in local and high-z infrared-bright galaxies and AGNs (e.g., Yang et al. 2013, 2016; Pensabene et al. 2022). The detections and upper limits of H<sub>2</sub>O(3<sub>12</sub>-2<sub>21</sub>) and H<sub>2</sub>O(3<sub>21</sub>-3<sub>12</sub>) for P215-16 and J1429+5447 follow the H<sub>2</sub>O-IR trend of other z ~ 6 quasars (Yang et al. 2019; Pensabene et al. 2022; D. A. Riechers et al. 2023, in preparation).

## 5. Discussion: The Heating Mechanisms of the Molecular Gas

With the NOEMA observations of the multiple CO transitions, we analyze the CO SLEDs of two z > 5.7 quasars. Both of them are among the most infrared-luminous objects found so far in the early Universe. The IR luminosities are of the order of a few 10<sup>12</sup> to 10<sup>13</sup> L<sub>☉</sub>, revealing massive star formation in a dusty ISM. J1429+5447 is a radio-loud quasar, which enables us to study the CO excitation close to a young and radio-loud AGN for the first time.

To further explore the possible origins of the different CO excitations observed in P215-16 and J1429+5447, we compare our observations with the CO SLEDs of local and high-z starburst galaxies and AGNs. The results are shown in Figure 6. The local and high-z galaxy samples are local (ultra) luminous infrared galaxies (“local (U)LIRGs”; Rosenberg et al. 2015), local star-forming and starburst galaxies (“local SF + SB”; Liu et al. 2015), lensed (sub)millimeter galaxies at z = 2-4 (“z = 2-4 SMGs”; Yang et al. 2017), and (sub) millimeter galaxies at z = 1.2-4.1 (“z = 1.2-4.1 SMGs”; Bothwell et al. 2013; Figure 6(a)). The representative local galaxies comprise three AGNs—NGC 1068 (Spinoglio et al. 2012), Mrk 231 (van der Werf et al. 2010), and NGC 6240 (Rosenberg et al. 2015)—and the starburst galaxy M82 (Weiß et al. 2005; Panuzzo et al. 2010; Figure 6(b)). The z ~ 6 comparison sources include the z ~ 6 quasars that have previously published CO SLEDs with detections of at least four transitions: J2310+1855 (Li et al. 2020b), J1148+5251 (Bertoldi et al. 2003; Walter et al. 2003; Beelen et al. 2006; Riechers et al. 2009; Gallerani et al. 2014), J0100+2802 (Wang et al. 2019), J0439+1634 (Yang et al. 2019), and PJ231-20 (Pensabene et al. 2021; Figure 6(c)). As the CO(6-5) line is not observed in PJ231-20, we normalize its CO SLED using the predicted CO(6-5) flux from the best-fit radiative transfer model of Pensabene et al. (2021). In addition, we collect starburst galaxies at z ~ 6, including SPT 0311 W, SPT 0311 E (Jarugula et al. 2021), and HFLS3 (Riechers et al. 2013), as well the lensed quasar APM 08279 at z = 3.9, which reveals the highest excited CO SLED ever published across redshifts (Weiß et al. 2007; Bradford et al. 2011; Figure 6(d)).

The CO SLEDs of local and high-z starburst galaxy samples typically peak at J ~ 6 and drop rapidly at high J, resulting in low line fluxes at J > 10 (Figure 6(a)). This is consistent with the gas heated by the FUV photons from young massive stars, indicating a PDR origin. On the other hand, the local AGNs and high-z quasars tend to have higher CO excitation with the CO SLEDs peaking at J ~ 6. The bright high-J (J ~ 9) CO transitions are usually associated with other heating mechanisms, e.g., X-ray heating from AGNs, cosmic-ray heating, and
<figure><img src="image_9.png" /><figcaption>**Figure 4.** CO SLED fitting results of the PDR (left), extreme PDR (middle), and XDR models (right) for P215–16 (top) and the PDR (left), XDR (middle), and PDR + XDR (right) models for J1429+5447 (bottom). The fitting results are shown as solid blue (PDR) and red (XDR) lines, and the data are shown as black squares and diamonds for P215–16 and J1429+5447, respectively. The sum of the PDR + XDR model is shown as a solid black line.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 5.** CLOUDY PDR model prediction for the [C II]158µm/[C I](2–1) line ratio (gray contours) with different gas densities (n(H₂)) and radiation fields (G₀). We also show the value of the [C II]158µm/LIR ratio (pink region), the lower limit of the [C II]158µm/[C I](2–1) ratio (magenta solid line), and the upper limit of the [C I](2–1)/LIR ratio (light blue region) for J1429+5447.</figcaption></figure>

mechanical heating by shocks in addition to the FUV heating from young massive stars (Figures 6(b) and (c); e.g., van der Werf et al. 2010; Spinoglio et al. 2012). As is suggested by Figure 6 and predictions from radiative transfer models, the CO

SLEDs of the starburst galaxies and AGNs/quasars are distinguishable through observations of the high-J (e.g., $J \gtrsim 9$) CO transitions (e.g., Vallini et al. 2019).

P215–16 exhibits a highly excited CO SLED, indicating an excitation level that is much higher than that of all local/high-z starburst galaxies (Figures 6(a) and (b)) and local AGNs (Figure 6(b)), with no turnover up to $J = 10$. It has the highest CO(10–9)/CO(6–5) ratio among all published quasars at $z \sim 6$ (Figure 6(c)) and the second highest excited CO ever across redshifts (less excited than APM 08279). We have demonstrated in Section 4.2 that either an XDR component with high gas density and intense radiation ($n(H_2) = 10^{4.3}$ cm⁻³, $F_x = 10^{1.1}$ erg s⁻¹ cm⁻²) or an extreme PDR of $n(H_2) = 10^{7.1}$ cm⁻³ and $G_0 = 10^{7.8}$ could explain such a highly excited CO SLED (Figure 4). Mechanical heating by shocks can also result in highly excited molecular CO, which was reported in some local galaxies (Rosenberg et al. 2015). Shocks are relatively more efficient in heating the gas than the dust, and a typical signature of shock-heated gas is a high CO/IR luminosity ratio of a few times $10^{-4}$ (e.g., Pellegrini et al. 2013). The CO/IR ratio (even for the brightest CO(10–9) transition) of P215–16 is $\sim 10^{-5}$, which disfavors the shock heating scenario. Enhanced cosmic-ray heating could also be responsible for high-J CO line excitation, which is sometimes indistinguishable from XDR. Vallini et al. (2019) investigated several CO excitation mechanisms in high-z galaxies and reported that the most extreme cosmic-ray-heated CO SLED peaked at $J = 9$. The high CO excitation of P215–16, which peaks at $J \ge 10$, is best explained by an XDR or extreme PDR component. In either case, AGNs could be the source powering
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 6.** CO SLEDs of P215−16 and J1429+5447 in comparison with (a) the mean value of local and high-z starburst galaxy samples, (b) local representative AGNs and starburst galaxies, (c) other quasars at $z \gtrsim 6$, and (d) starburst galaxies at $z \gtrsim 6$ and the lensed quasar APM 08279 at $z = 3.9$. The dashed line represents the thermalized CO excitation with constant brightness temperature on the Rayleigh–Jeans scale, i.e., $L_{CO(J-J-1)} = L'_{CO(1-0)}$.</figcaption></figure>

the high CO excitation, by contributing to the X-ray emission or intense FUV emission as suggested by the results of the XDR and extreme PDR models.

In contrast, J1429+5447 exhibits a relatively lower CO excitation, which is comparable to that measured for local and high-z starburst galaxies (Figures 6(a) and (b)) and $z \gtrsim 6$ starburst galaxies, e.g., SPT 0311 W, SPT 0311 E, and HFLS3 (Figure 6(d)). It has the lowest $CO(J-J-1)_{[J \ge 7]}$/CO(6–5) line ratios among the $z \sim 6$ quasars with published CO SLEDs (Figure 6(d)). The CO SLED could be explained by a single PDR model with extreme gas densities of $n(\text{H}_2) = 10^{6.0} \text{ cm}^{-3}$ and radiation of $G_0 = 10^{5.4}$.

We also explore whether the CO SLED of J1429+5447 can be fitted with an XDR model included, given the luminous X-ray detection in this quasar. The CO SLED of J1429+5447 fitted by an XDR or a PDR + XDR model prefers an XDR component with moderate density ($n(\text{H}_2) = 10^{3.7-4.3} \text{ cm}^{-3}$) and illuminated by an intense X-ray radiation field ($F_x = 10^{1.5-2.0} \text{ erg s}^{-1} \text{ cm}^{-2}$). This is at the lower edge of the range of molecular gas density found in $z \sim 6$ quasars (e.g., Li et al. 2020b; Pensabene et al. 2021). As J1429+5447 is a radio-loud object, it is possible that the radio jet expels the gas from the host galaxy, resulting in lower gas density. This is consistent with the indication of a possible outflow found in previous

[C II]$_{158\mu m}$ observations (Khusanova et al. 2022). The outflowing gas mass constrained from the broad [C II]$_{158\mu m}$ component is $1.2 \times 10^{10} M_\odot$, assuming $\alpha_{[CII]} = 30 M_\odot L_\odot^{-1}$ (Zanella et al. 2018). The $OH^+(1_1-0_1)$ absorption is a sensitive tracer of outflows in galaxies and quasars (e.g., Shao et al. 2022). Unfortunately, the low S/N of our $OH^+(1_1-0_1)$ spectrum prevents us from deriving useful constraints on the outflowing molecular gas mass, and deeper observations will be needed to obtain an independent estimate. In this work, we only include a simple plane-parallel geometry and include no dust torus attenuation. A low CO excitation can be produced in an XDR model of sophisticated geometry and high dust attenuation even when the X-ray flux and gas density are high (e.g., Vallini et al. 2019). Future CO SLED observations of more radio-loud quasars will possibly enable us to study whether the low excitation is rare/common in the radio-loud population.

## 6. Summary

We present NOEMA observations and report detections of multiple CO transitions (from CO(5–4) to CO(10–9)), H$_2$O(3$_{21}-3_{12}$), and the underlying continuum emission from the host galaxies of the quasars P215−16 at $z = 5.78$ and J1429 +5447 at $z = 6.18$.
Adopting a general opacity dust SED model, we derive dust temperatures in the range of ~50–110 K for the two quasar host galaxies. The derived infrared luminosities of P215–16 and J1429+5447 differ by factors of ~2.5 and 5.3, respectively, depending on the best-fit models used.

The quasar P215–16 exhibits highly excited molecular CO with a CO SLED peaking at $J \ge 10$. It has the highest CO($J-J=1$)$_{J \ge 7}$/CO(6–5) ratio among the $z \sim 6$ quasars with published CO SLEDs known to date. Analysis of the CO SLED with models suggests that a “dense” ($n(\text{H}_2) = 10^{4.3} \text{ cm}^{-3}$) XDR component illuminated by an intense X-ray radiation field ($F_\chi = 10^{1.1}$) or a “dense” PDR component with $n(\text{H}_2) = 10^{7.1} \text{ cm}^{-3}$ and $G_0 = 10^{7.8}$ reproduces the observed CO SLED. In either case, the central AGN is likely to contribute to the radiation field to excite the high-$J$ CO transitions.

The CO SLED of the radio-loud quasar J1429+5447 reveals the lowest CO excitation among all published quasar CO SLEDs at $z \sim 6$. This is the first reported CO SLED of a radio-loud quasar at this redshift and is comparable to those of local AGNs and starburst galaxies at $z \gtrsim 6$. The CO SLED could be explained by a PDR component. In addition, we also propose two possibilities where XDR is considered to explain the low CO excitation: (1) a luminous X-ray radiation field illuminating a moderate-density gas, which may be attributed to a result of the gas being expelled by an outflow, and (2) an XDR model of sophisticated geometry taking into account the X-ray attenuation by a dust torus. In this case, low CO excitation can be reproduced in a model with high dust attenuation even when the X-ray fluxes and gas densities are high.

In addition, we detect bright $H_2O(3_{21}-3_{12})$ emission in P215–16. This places P215–16 exactly on the linear relation between the $H_2O(3_{21}-3_{12})$ and IR luminosities recently found for other quasars at $z \sim 6$ and local/high-$z$ infrared-bright galaxies.

Future observations of the CO SLEDs toward more radio-loud quasars at $z \gtrsim 6$, possibly up to ($J \gtrsim 10$), will be critical to systematically reveal the physical conditions and heating mechanisms of the molecular gas in the host galaxies of the radio-loud quasar population in the early Universe. Comparisons between the radio-loud and radio-quiet populations will shed light on the possible different/similar impacts of the radio-quiet (loud) AGNs on the molecular gas properties.

## Acknowledgments

We acknowledge support from the National Natural Science Foundation of China (NSFC) with grant Nos. 12173002, 11991052. This work is based on observations carried out under project Nos. W0C3, W18EE, and subsequently S20CW with the IRAM NOEMA Interferometer. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). We gratefully acknowledge the support from the Shuimu Tsinghua Scholar Program of Tsinghua University. A.P. acknowledges support from Fondazione Cariplo grant No. 2020-0902. This work has been supported by the National Key R & D Program of China (grant No. 2018YFA0404503), the National Science Foundation of China (grant No. 12073014), and the science research grants from the China Manned Space Project with Nos. CMS-CSST2021-A05, CMS-CSST-2021-A06. Y.S. acknowledges supports from the NSFC with grant No. 12303001 and the Beihang University grant with No. ZG216S2305.

## ORCID iDs

Jianan Li https://orcid.org/0000-0002-1815-4839
Ran Wang https://orcid.org/0000-0003-4956-5742
Antonio Pensabene https://orcid.org/0000-0001-9815-4953
Fabian Walter https://orcid.org/0000-0003-4793-7880
Bram P. Venemans https://orcid.org/0000-0001-9024-8322
Roberto Decarli https://orcid.org/0000-0002-2662-8803
Eduardo Bañados https://orcid.org/0000-0002-2931-7824
Pierre Cox https://orcid.org/0000-0003-2027-8221
Roberto Neri https://orcid.org/0000-0002-7176-4046
Alain Omont https://orcid.org/0000-0002-4721-3922
Zheng Cai https://orcid.org/0000-0001-8467-6478
Yana Khusanova https://orcid.org/0000-0002-7220-397X
Fuxiang Xu https://orcid.org/0000-0003-0754-9795
Dominik Riechers https://orcid.org/0000-0001-9585-1462
Jeff Wagg https://orcid.org/0000-0002-1815-4839
Yali Shao https://orcid.org/0000-0002-1478-2598
Yuanqi Liu https://orcid.org/0000-0001-9321-6000
Karl M. Menten https://orcid.org/0000-0001-6459-0669
Qiong Li https://orcid.org/0000-0002-3119-9003
Xiaohui Fan https://orcid.org/0000-0003-3310-0131

## References

Becker, R. H., White, R. L., & Helfand, D. J. 1995, ApJ, 450, 559
Beelen, A., Cox, P., Benford, D. J., et al. 2006, ApJ, 642, 694
Bertoldi, F., Carilli, C. L., Cox, P., et al. 2003, A&A, 406, L55
Bothwell, M. S., Smail, I., Chapman, S. C., et al. 2013, MNRAS, 429, 3047
Bradford, C. M., Bolatto, A. D., Maloney, P. R., et al. 2011, ApJL, 741, L37
Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105
Carniani, S., Gallerani, S., Vallini, L., et al. 2019, MNRAS, 489, 3939
da Cunha, E., Hodge, J. A., Casey, C. M., et al. 2021, ApJ, 919, 30
De Rosa, G., Venemans, B. P., Decarli, R., et al. 2014, ApJ, 790, 145
Decarli, R., Pensabene, A., Diaz-Santos, T., et al. 2023, A&A, 673, A157
Decarli, R., Walter, F., Venemans, B. P., et al. 2017, Natur, 545, 457
Decarli, R., Walter, F., Venemans, B. P., et al. 2018, ApJ, 854, 97
Di Mascia, F., Gallerani, S., Behrens, C., et al. 2021, MNRAS, 503, 2349
Esposito, F., Vallini, L., Pozzi, F., et al. 2022, MNRAS, 512, 686
Foreman-Mackey, D., Hogg, D. W., Lang, D., et al. 2013, PASP, 125, 306
Gallerani, S., Ferrara, A., Neri, R., & Maiolino, R. 2014, MNRAS, 445, 2848
Gilli, R., Norman, C., Calura, F., et al. 2022, A&A, 666, A17
Greve, T. R., Leonidaki, I., Xilouris, E. M., et al. 2014, ApJ, 794, 142
Guilloteau, S., & Lucas, R. 2000, in ASP Conf. Ser. 217, Imaging at Radio through Submillimeter Wavelengths, ed. J. G. Mangum & S. J. E. Radford (San Francisco, CA: ASP), 299
Habing, H. J. 1968, BAN, 19, 421
Jarugula, S., Vieira, J. D., Weiss, A., et al. 2021, ApJ, 921, 97
Kamenetzky, J., Rangwala, N., Glenn, J., et al. 2016, ApJ, 829, 93
Khorunzhev, G. A., Meshcheryakov, A. V., Medvedev, P. S., et al. 2021, AstL, 47, 123
Khusanova, Y., Bañados, E., Mazzucchelli, C., et al. 2022, A&A, 664, A39
Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511
Li, J., Wang, R., Cox, P., et al. 2020a, ApJ, 900, 131
Li, J., Wang, R., Riechers, D., et al. 2020b, ApJ, 889, 162
Li, Q., Wang, R., Fan, X., et al. 2020, ApJ, 900, 12
Liu, D., Gao, Y., Isaak, K., et al. 2015, ApJL, 810, L14
Lu, N., Zhao, Y., Xu, C. K., et al. 2014, ApJL, 787, L23
McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565
Medvedev, P., Sazonov, S., Gilfanov, M., et al. 2020, MNRAS, 497, 1842
Meijerink, R., Kristensen, L. E., Weiß, A., et al. 2013, ApJL, 762, L16
Meijerink, R., & Spaans, M. 2005, A&A, 436, 397
Meijerink, R., Spaans, M., & Israel, F. P. 2007, A&A, 461, 793
Meyer, R. A., Walter, F., Cicone, C., et al. 2022, ApJ, 927, 152
Migliori, G., Siemiginowska, A., Sobolewska, M., et al. 2023, MNRAS, 524, 1087
Mingozzi, M., Vallini, L., Pozzi, F., et al. 2018, MNRAS, 474, 3640
Montoya Arroyave, I., Cicone, C., Makroleivaditi, E., et al. 2023, A&A, 673, A13
Morganson, E., De Rosa, G., Decarli, R., et al. 2012, AJ, 143, 142
Panuzzo, P., Rangwala, N., Rykala, A., et al. 2010, A&A, 518, L37
Pellegrini, E. W., Smith, J. D., Wolfire, M. G., et al. 2013, ApJL, 779, L19
Pensabene, A., Carniani, S., Perna, M., et al. 2020, A&A, 637, A84

Pensabene, A., Decarli, R., Bañados, E., et al. 2021, A&A, 652, A66

Pensabene, A., van der Werf, P., Decarli, R., et al. 2022, A&A, 667, A9

Pereira-Santaella, M., Spinoglio, L., van der Werf, P. P., et al. 2014, A&A, 566, A49

Pozzi, F., Vallini, L., Vignali, C., et al. 2017, MNRAS, 470, L64

Ramos Almeida, C., Bischetti, M., García-Burillo, S., et al. 2022, A&A, 658, A155

Riechers, D. A. 2011, ApJ, 730, 108

Riechers, D. A., Bradford, C. M., Clements, D. L., et al. 2013, Natur, 496, 329

Riechers, D. A., Carilli, C. L., Maddalena, R. J., et al. 2011, ApJL, 739, L32

Riechers, D. A., Walter, F., Carilli, C. L., & Lewis, G. F. 2009, ApJ, 690, 463

Riechers, D. A., Walter, F., Carilli, C. L., et al. 2006, ApJ, 650, 604

Rosenberg, M. J. F., van der Werf, P. P., Aalto, S., et al. 2015, ApJ, 801, 72

Shao, Y., Wang, R., Carilli, C. L., et al. 2019, ApJ, 876, 99

Shao, Y., Wang, R., Weiss, A., et al. 2022, A&A, 668, A121

Spaans, M. 2008, EAS, 31, 47

Spinoglio, L., Pereira-Santaella, M., Busquet, G., et al. 2012, ApJ, 758, 108

Uzgil, B. D., Bradford, C. M., Hailey-Dunsheath, S., et al. 2016, ApJ, 832, 209

Valentino, F., Daddi, E., Puglisi, A., et al. 2021, A&A, 654, A165

Vallini, L., Tielens, A. G. G. M., Pallottini, A., et al. 2019, MNRAS, 490, 4502

van der Werf, P. P., Isaak, K. G., Meijerink, R., et al. 2010, A&A, 518, L42

Venemans, B. P., Decarli, R., Walter, F., et al. 2018, ApJ, 866, 159

Venemans, B. P., Walter, F., Neeleman, M., et al. 2020, ApJ, 904, 130

Venemans, B. P., Walter, F., Zschaechner, L., et al. 2016, ApJ, 816, 37

Villarreal Hernández, A. C., & Andernach, H. 2018, arXiv:1808.07178

Walter, F., Bertoldi, F., Carilli, C., et al. 2003, Natur, 424, 406

Walter, F., Neeleman, M., Decarli, R., et al. 2022, ApJ, 927, 21

Wang, F., Wang, R., Fan, X., et al. 2019, ApJ, 880, 2

Wang, R., Carilli, C. L., Wagg, J., et al. 2008, ApJ, 687, 848

Wang, R., Wagg, J., Carilli, C. L., et al. 2011, ApJL, 739, L34

Wang, R., Wagg, J., Carilli, C. L., et al. 2013, ApJ, 773, 44

Wang, R., Wu, X. B., Neri, R., et al. 2016, ApJ, 830, 53

Weiß, A., Downes, D., Neri, R., et al. 2007, A&A, 467, 955

Weiß, A., Walter, F., & Scoville, N. Z. 2005, A&A, 438, 533

Willott, C. J., Bergeron, J., & Omont, A. 2015, ApJ, 801, 123

Willott, C. J., Delorme, P., Reylé, C., et al. 2010, AJ, 139, 906

Yang, C., Gao, Y., Omont, A., et al. 2013, ApJL, 771, L24

Yang, C., Omont, A., Beelen, A., et al. 2016, A&A, 595, A80

Yang, C., Omont, A., Beelen, A., et al. 2017, A&A, 608, A144

Yang, J., Venemans, B., Wang, F., et al. 2019, ApJ, 880, 153

Zanella, A., Daddi, E., Magdis, G., et al. 2018, MNRAS, 481, 1976