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

# Constraining the Excitation of Molecular Gas in Two Quasar-starburst Systems at $z \sim 6$

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

<sup>1</sup>Department of Astronomy, School of Physics, Peking University, Beijing 100871, 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

<sup>3</sup>Department of Astronomy, Tsinghua University, Beijing 100084, People's Republic of China

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

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

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

<sup>7</sup>School of Space and Earth Sciences, Beihang University, Beijing, People's Republic of China

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

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

<sup>10</sup>Argelander-Institut für Astronomie, University at Bonn, Auf dem Hügel 71, 53121 Bonn, Germany

<sup>11</sup>INAF-Osservatorio Astronomico di Trieste, Via G. Tiepolo 11, I-34143 Trieste, Italy

<sup>12</sup>IFPU-Institut for Fundamental Physics of the Universe, via Beirut 2, I-34151 Trieste, Italy

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

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

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

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

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

<sup>18</sup>World Food Programme, No. 2 Jawatte Ave., Colombo 00500, Sri Lanka

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

<sup>20</sup>Steward Observatory, University of Arizona, 933 North Cherry Ave., Tucson, AZ 85721, USA

Received 2024 August 13; revised 2024 October 27; accepted 2024 October 30; published 2024 December 12

## Abstract

We present Northern Extended Millimeter Array observations of CO(8–7), (9–8), and (10–9) lines, as well as the underlying continuum for two far-infrared luminous quasars: SDSS J2054-0005 at $z = 6.0389$ and SDSS J0129-0035 at $z = 5.7788$. Both quasars were previously detected in CO (2–1) and (6–5) transitions, making them candidates for studying the CO spectral line energy distribution (SLED) of quasars at $z \sim 6$. Utilizing the radiative transfer code CLOUDY, we fit the CO SLED with two heating mechanisms, including the photodissociation region (PDR) and X-ray-dominated region (XDR) for both objects. The CO SLEDs can be fitted by either a dense PDR component with an extremely strong far-ultraviolet radiation field (gas density $n_H \sim 10^6$ cm$^{-3}$ and field strength $G_0 \gtrsim 10^6$) or a two-component model including a PDR and an XDR. However, the line ratios, including $L_{\text{TIR}}$ and previous [C II]$_{158 \, \mu m}$ and [C I]$_{369 \, \mu m}$ measurements, argue against a very high PDR radiation field strength. Thus, the results prefer a PDR+XDR origin for the CO SLED. The excitation of the high-$J$ CO lines in both objects is likely dominated by the central active galactic nucleus (AGN). We then check the CO (9–8)-to-(6–5) line luminosity ratio $r_{96}$ for all $z \sim 6$ quasars with available CO SLEDs (seven in total) and find that there are no clear correlations between $r_{96}$ and both $L_{\text{FIR}}$ and the AGN UV luminosities. This further demonstrates the complexity of the CO excitation powered by both the AGN and nuclear star formation in these young quasar host galaxies.

Unified Astronomy Thesaurus concepts: AGN host galaxies (2017); Galaxies (573); Quasars (1319); Molecular gas (1073); Submillimeter astronomy (1647)

## 1. Introduction

Over the past two decades, more than 500 quasars at $z \gtrsim 5.3$ have been discovered through optical/near-infrared (NIR) wide-field multiband surveys (X. Fan et al. 2000; L. Jiang et al. 2015; F. Wang et al. 2019; J. Yang et al. 2019; Y. Matsuoka et al. 2022; X. Fan et al. 2023). These quasars represent the actively accreting supermassive black holes at the earliest epoch (e.g., L. Jiang et al. 2007, 2016; Y. Shen et al. 2019). A large percent of the high-$z$ quasars were detected in far-infrared (FIR) dust continuum, while most high-$z$ FIR-luminous quasars were detected in various CO emission lines at (sub)millimeter wavelengths, revealing dynamical masses of $\sim 10^{10}-10^{11} M_\odot$,

star formation rates (SFRs) of $\sim 100-1000 M_\odot$ yr$^{-1}$, and molecular gas reservoirs of $\sim 10^{10}-10^{11} M_\odot$ within a few kiloparsec scales in their host galaxies (e.g., F. Bertoldi et al. 2003; D. A. Riechers et al. 2009; R. Wang et al. 2010, 2011; C. Carilli & F. Walter 2013; R. Wang et al. 2013, 2016; B. P. Venemans et al. 2017, 2018, 2019; A. Pensabene et al. 2020; B. P. Venemans et al. 2020; M. Neeleman et al. 2021; F. Walter et al. 2022; Y. Shao et al. 2022; M. Neeleman et al. 2023).

The interstellar medium (ISM) in galaxies serves as a tool for diagnosing various heating and excitation mechanisms based on its abundant spectral lines and continuum emission from the thermal dust. Recently, (sub)millimeter and radio facilities, such as the Atacama Large Millimeter/submillimeter Array (ALMA) and the Northern Extended Millimeter Array (NOEMA), have enabled us to measure the multiphase ISM emission from high-$z$ quasar host galaxies with unprecedented

<figure><img src="image_2.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.
sensitivity and resolution. An essential tracer of the ISM in these galaxies is the [C II]$_{158}$ µm emission line, one of the primary coolants of the ISM, originating from the ionized and neutral gas surrounding regions of intense star formation. Notably, observations have revealed significant [C II] emission in $z \sim 6$ quasar host galaxies, providing compelling evidence for widespread photodissociation region (PDR) heating and ongoing star formation activity (B. P. Venemans et al. 2012; R. Wang et al. 2013; B. P. Venemans et al. 2016; R. Wang et al. 2016; R. Decarli et al. 2018; B. P. Venemans et al. 2020; Y. Shao et al. 2022).

The CO emission lines have been used to study the physical conditions in a variety of local and high-$z$ systems (e.g., A. Weiß et al. 2005; D. A. Riechers et al. 2006; A. Weiß et al. 2007; C. M. Bradford et al. 2009, 2011; L. Spinoglio et al. 2012; S. Gallerani et al. 2014; M. J. F. Rosenberg et al. 2015; F. Wang et al. 2019; J. Yang et al. 2019; J. Li et al. 2020, 2024). Low-$J$ (e.g., $J \lesssim 2$) CO emission lines trace the bulk of the molecular gas content. The relations between the line luminosities of different CO transitions and the FIR luminosity are widely discussed, which reflects the tight connection between the dense molecular gas and star formation (e.g., T. R. Greve et al. 2014; D. Liu et al. 2015; J. Kamenetzky et al. 2016). The low and mid-$J$ CO transitions (e.g., $3 \lesssim J \lesssim 8$) were observed in samples of high-$z$ submillimeter galaxies and starburst quasar hosts which trace the PDR region powered by the ultraviolet (UV) photons from newly formed high-mass stars (e.g., M. S. Bothwell et al. 2013; C. Yang et al. 2017; R. Decarli et al. 2022). Bright high-$J$ (e.g., $J \gtrsim 9$) CO transitions were detected in quasar host galaxies from local to the highest redshift, suggesting highly excited molecular gas associated with the powerful active galactic nucleus (AGN) activity (C. M. Bradford et al. 2003; L. Spinoglio et al. 2012; R. Meijerink et al. 2013; S. Gallerani et al. 2014; J. Li et al. 2020; A. Pensabene et al. 2021).

By integrating observational data with modeling techniques, we can effectively determine representative values for crucial ISM physical parameters, including the gas kinetic temperature ($T_{\text{kin}}$), gas density ($n_H$), and interstellar radiation field intensity (UV intensity $G_0$ from PDRs or X-rays intensity $F_X$ from AGN). A common approach employs the large velocity gradient model for single/multiple-component fitting of CO spectral line energy distributions (SLEDs; A. Weiß et al. 2005, 2007; D. A. Riechers et al. 2009; L. Spinoglio et al. 2012; F. Wang et al. 2019; J. Yang et al. 2019; J. Li et al. 2020). This method has unveiled the properties of molecular gas through a combination of cold ($T_{\text{kin}} \sim 20–100$ K) low-density ($n_H \sim 10^3–10^4$ cm$^{-3}$) and hot ($T_{\text{kin}} \sim 100–300$ K) high-density ($n_H \sim 10^4–10^5$ cm$^{-3}$) components. Radiation transfer models for PDRs and X-ray-dominated regions (XDRs) are explored to further constrain the excitation power and mechanisms. For example, studies of the low- and mid-$J$ CO transitions in the starburst systems reveal PDR conditions with a typical interstellar radiation field strength of $G_0 \sim 10^2–10^3$ and $n_H \sim 10^2–10^4$ cm$^{-3}$ (P. P. van der Werf et al. 2010; R. Meijerink 2011; L. Spinoglio et al. 2012; T. R. Greve et al. 2014; D. Liu et al. 2015; M. J. F. Rosenberg et al. 2015). In comparison, high-$J$ CO emission in high-$z$ quasars requires an extra, highly dense ($n_H \sim 10^4–10^5$ cm$^{-3}$) XDR component or other possible scenarios such as mechanical heating, cosmic rays, and shocks (R. Meijerink et al. 2007; L. Spinoglio et al. 2012; S. Gallerani et al. 2014; A. Pensabene et al. 2021; R. Decarli et al. 2022; J. Li et al. 2024). This highlights the

potential contribution of AGN to ISM heating, particularly in the early Universe.

The studies of CO SLEDs for $z \sim 6$ quasars have been mainly focused on the most optically luminous ($M_{1450} \lesssim -26$) objects so far. To better address the impact of AGN on the nuclear ISM, it is essential to expand the study to less UV luminous and FIR bright systems. In this work, we analyzed the CO SLED in the host galaxies of two FIR-luminous quasars, SDSS J2054-0005 at $z=6.0379$ and SDSS J0129-0035 at $z=5.7787$. These two objects are fainter in the optical ($M_{1450}=-26.1$ and $-24.4$ mag; L. Jiang et al. 2016) compared to quasars at $z \sim 6$ from previous CO SLED studies (S. Gallerani et al. 2014; F. Wang et al. 2019; J. Yang et al. 2019; J. Li et al. 2020; A. Pensabene et al. 2021; J. Li et al. 2024). In particular, J0129-0035 is about two magnitudes fainter than the previously faintest one, J1429+5447 (J. Li et al. 2024). The two quasars were detected in the dust continuum, as well as in the [C II]$_{158}$ µm, CO(6–5), and CO(2–1) line (R. Wang et al. 2013, 2019; Y. Shao et al. 2019). The derived FIR luminosities (integrated over the range $42.5–122.5$ µm) are about $8 \times 10^{12}$ L$_\odot$ for J2054-0005 and $5 \times 10^{12}$ L$_\odot$ for J0129-0035, placing them among the most FIR-luminous objects at $z \sim 6$. The ALMA continuum and [C II]$_{158}$ µm images reveal compact nuclear star formation on scales of a few kiloparsecs (R. Wang et al. 2013, 2019). Thus, they provide ideal targets for studying ISM excitation in environments with intense star formation and less luminous AGNs. In this paper, we present new observations of CO(8–7), CO(9–8), and CO(10–9) lines, as well as the dust continuum obtained with NOEMA toward J2054-0005 and J0129-0035. The goal is to investigate and further constrain the molecular gas excitation mechanisms in their nuclear regions and explore the dust heating and star formation activity.

The paper is organized as follows. In Section 2, we describe the sample, the observations, and the data reduction. In Section 3, we present our results. In Section 4, we first analyze the dust temperature, mass, and SFR for J2054-0005 and J0129-0035 through SED fitting of the dust continuum. We then use the CLOUDY radiative transfer code (version c.17.01, G. J. Ferland et al. 2017) to investigate the physical conditions of the molecular gas through an analysis of the CO SLED. Finally, we compare the CO excitation of J2054-0005 and J0129-0035 with local and high-$z$ AGN and galaxy samples. In Section 5, we summarize our results.

Throughout this work, we adopt a ΛCDM cosmology with $H_0 = 70$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_m = 0.3$, and $\Omega_\Lambda = 0.7$.

## 2. Observations

We carried out observations of the CO(8–7) ($\nu_{\text{rest}} = 921.80$ GHz), CO(9–8) ($\nu_{\text{rest}} = 1036.9$ GHz), CO(10–9) ($\nu_{\text{rest}} = 1152.0$ GHz) emission lines, together with the underlying continuum from SDSS J2054-0005 and SDSS J0129-0035 using NOEMA (Project ID: S21DD). The observations were carried out from 2021 September to November, using 9 to 10 antennas in C/D configurations. We used the Band 2 receiver and the PolyFix correlator covering the CO lines in two setups. In one setup, we observe the CO (8–7) line in the lower sideband (LSB) and the CO (9–8) line in the upper sideband (USB), and in the other setup, we put the CO (10–9) line in the USB and CO (9–8) in the LSB. Each setup covers a total bandwidth of 7.74 GHz and a channel width of 2 MHz. MWC 349 was used as the flux calibrator for both J2054-0005 and J0129-0035. The phase/amplitude calibrators were J2059+034 for J2054-0005 and
<table><caption>**Table 1**
Summary of the Observations</caption><thead><tr><th>Line</th><th>t<sub>on_source</sub><br>(hr)</th><th>n<sub>ant</sub></th><th>Beam<br>(arcsec<sup>2</sup>)</th><th>rms (Channel Width)<br>mJy Beam<sup>-1</sup> (km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td colspan="5"><strong>J2054-0005</strong></td></tr><tr><td>CO(8-7)</td><td>7.89</td><td>9-10</td><td>3.78 &times; 3.01</td><td>0.31 (64)</td></tr><tr><td>CO(9-8)</td><td>15.89</td><td>9-10</td><td>2.41 &times; 1.06</td><td>0.26 (65)</td></tr><tr><td>CO(10-9)</td><td>8.00</td><td>9-10</td><td>2.91 &times; 1.57</td><td>0.36 (77)</td></tr><tr><td colspan="5"><strong>J0129-0035</strong></td></tr><tr><td>CO(8-7)</td><td>6.81</td><td>9-10</td><td>4.73 &times; 2.46</td><td>0.31 (66)</td></tr><tr><td>CO(9-8)</td><td>11.41</td><td>9-10</td><td>3.43 &times; 2.54</td><td>0.27 (63)</td></tr><tr><td>CO(10-9)</td><td>4.60</td><td>9-10</td><td>3.23 &times; 2.52</td><td>0.59 (74)</td></tr></tbody></table>

Notes. Column (1): line's name. Column (2): on-source integration time. Column (3): number of antennas. Column (4): beam size in FWHM. Column (5): rms and corresponding channel width.

J0106+013 for J0129-0035. The typical calibration uncertainty is <10% in Band 2.

Data reduction was performed using the Grenoble Image and Line Data Analysis System software (S. Guilloteau & R. Lucas 2000) packages CLIC and MAPPING. We extracted the continuum from all channels free of spectral lines in the *uv* plane with UV\_AVERAGE. The *uv* table for spectral lines was produced via UV\_SUBTRACT with the subtraction of the underlying continuum. The *uv* table of the continuum and spectral lines underwent cleaning with the HOGBOM algorithm, while NATURAL weighting was applied to achieve the maximum signal-to-noise ratio (S/N). The CO(9–8) line data from two tunings were combined to improve the S/N. The processed data were binned to spectral resolutions of 64–77 km s⁻¹ to optimize the sensitivity. The details are listed in Table 1.

### 3. Results

With the new NOEMA observations, we detected the CO (8–7), CO (9–8), and CO (10–9) lines from J2054-0005 in the velocity-integrated maps up to 6σ, 8σ, and 5σ, respectively (Figure 1). We also detected the CO (8–7) and CO (9–8) lines from J0129-0035 up to 4σ and 5σ, respectively (Figure 1). The sources are spatially unresolved in the velocity-integrated intensity maps. The unresolved continuum emission was also detected at high S/N by taking the peak surface brightness value on the continuum map averaging over line-free channels (Figure 1). We summarize the measurements of the CO lines in Table 2 and the continuum in Table 3.

#### 3.1. J2054-0005

This quasar was initially detected in Sloan Digital Sky Survey (SDSS) stripe 82 with a magnitude of $m_{1450\AA} = 20.60$ mag (L. Jiang et al. 2008). Previous Very Large Array (VLA) observations (Y. Shao et al. 2019) detected the CO (2–1) line with a flux of $0.06 \pm 0.01$ Jy km s⁻¹ and an FWHM line width of $270 \pm 47$ km s⁻¹. The CO (6–5) line detection in previous Plateau de Bure Interferometer (PdBI) observations (R. Wang et al. 2010) reported a line flux of $0.34 \pm 0.07$ Jy km s⁻¹ and an FWHM line width of $356 \pm 106$ km s⁻¹. However, new ALMA observations (R. Tripodi et al. 2024) show a line flux of $0.29 \pm 0.04$ Jy km s⁻¹ and an FWHM line width of $229 \pm 20$ km s⁻¹. Considering that the ALMA data have significantly improved sensitivity and more accurate measurement of the line

profile, we adopt the ALMA CO (6–5) line measurements for our analysis throughout this paper. R. Decarli et al. (2022) reported a much broader line width for the CO (7–6) and [C I]$_{369\mu m}$ lines. This inconsistency may be due to the low S/N in their line spectra and the strong noise close to the atmosphere features. We, therefore, exclude the CO (7–6) emission line from the CO SLED analysis in Section 4.2. The intensity maps of CO(8–7), CO(9–8), and CO(10–9), created by integrating over $1.5 \times \text{FWHM}_{\text{line}}$ around the redshifted CO line, are presented in the top-row panels of Figure 1 from left to right. The corresponding continuum images are displayed in the panels of the second row of Figure 1. We also fitted a single Gaussian profile from the peak position to the line spectra and obtained line fluxes of $0.36 \pm 0.07$ Jy km s⁻¹, $0.27 \pm 0.04$ Jy km s⁻¹, and $0.24 \pm 0.07$ Jy km s⁻¹ for CO(8–7), CO(9–8), and CO(10–9), respectively (Figure 2). The spectra of CO (8–7) and (10–9) have lower S/N compared to that of the other two transitions and the line widths determined by the fitting suffer large uncertainties. In particular, the line width derived directly from the least-square fitting of the CO (8–7) line spectrum yields a value of $331 \pm 62$ km s⁻¹, which is much larger than that determined with the CO (6–5) and (9–8) emission lines (Table 2). Therefore, for these two lines, we employ a Bayesian Markov Chain Monte Carlo (MCMC) method and adopt the CO (6–5) line width of $229$ km s⁻¹ as a prior to fit the line width. To ensure consistency, we also applied MCMC to the CO(9–8) line and found that the result is consistent with the least-square fitting within the uncertainties.

In Figure 1, the peak positions of the emission are marked with green–yellow crosses, all of which coincide with the optical position ($20^h54^m06.4966 \pm 0^s0.347, -00^\circ05'14.5'500 \pm 0^\prime521$, Y. Shao et al. 2019). The CO(8–7), CO(9–8), and CO(10–9) emission lines show consistent line widths with those reported for CO(6–5) (R. Tripodi et al. 2024) and [C II]$_{158\mu m}$ (B. P. Venemans et al. 2020).

#### 3.2. J0129-0035

This quasar was initially identified in SDSS stripe 82 with a magnitude of $m_{1450\AA} = 22.28$ mag (L. Jiang et al. 2009). In previous VLA observations (Y. Shao et al. 2019), the CO (2–1) line was detected with a flux of $0.036 \pm 0.005$ Jy km s⁻¹ and an FWHM line width of $195 \pm 41$ km s⁻¹. The CO (6–5) line detection in previous PdBI observations (R. Wang et al. 2011) reported a line flux of $0.37 \pm 0.07$ Jy km s⁻¹ and an FWHM line width of $283 \pm 87$ km s⁻¹. However, new ALMA observations (R. Tripodi 2024) detected the CO (6–5) line with a line flux of $0.18 \pm 0.03$ Jy km s⁻¹ and an FWHM line width of $183 \pm 16$ km s⁻¹. Similar to J2054-0005, we adopt the ALMA measurements for the CO (6–5) line in our analysis throughout this paper. The intensity maps of CO(8–7), CO(9–8), and CO(10–9), which are also created by integrating over a velocity range of $1.5 \times \text{FWHM}_{\text{line}}$ around the redshifted CO line, are showcased in the third-row panels of Figure 1. The corresponding continuum images are displayed in the bottom-row panels. For the detections, we extracted the line spectra from the peak position and fitted a single Gaussian profile to the line spectra (Figure 2). The measured line fluxes are $0.21 \pm 0.09$ Jy km s⁻¹ for the CO (8–7) line and $0.21 \pm 0.05$ Jy km s⁻¹ for CO (9–8). Note that the S/N of the CO (8–7) line spectra is insufficient for a robust fitting. Thus, we follow a similar formula as we did for J2054-0005, using the ALMA CO (6–5) FWHM line width of
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Line intensity maps of CO(8–7), (9–8), (10–9) (from left to right) and corresponding continuum of two objects observed by NOEMA. The top two row panels are for J2054-0005, and the bottom two rows are for J0129-0035. The green–yellow cross on the line intensity map and continuum represents the optical position. The shaded white ellipse in the bottom-left corner of each panel shows the FWHM of the beam. The contours denote the line map contours to [-2, 2, 3, 4, 6] ×σ and the continuum contours to [-2, 2, 4, 8, 16] ×σ.</figcaption></figure>
<table><caption>**Table 2**

Spectral Line Observations</caption><thead><tr><th>Line</th><th>z<sub>line</sub></th><th>FWHM (km s<sup>-1</sup>)</th><th>Flux (Jy km<sup>-1</sup>)</th><th>Beam Size (arcsec)</th><th>Luminosity (10<sup>8</sup>L<sub>☉</sub>)</th><th>Facilities</th><th>References</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></tr></thead><tbody><tr><td colspan="8"><strong>J2054-0005</strong></td></tr><tr><td>[C II]</td><td>6.0389 ± 0.0001</td><td>236 ± 12</td><td>3.23 ± 0.14</td><td>0.15 × 0.11</td><td>30.8 ± 1.4</td><td>ALMA</td><td>V20</td></tr><tr><td>[C I]</td><td>6.0397 ± 0.0003</td><td>487 ± 115</td><td>&lt; 0.261</td><td>4.0 × 2.6</td><td>&lt;1.1</td><td>NOEMA</td><td>D22</td></tr><tr><td>CO(2-1)</td><td>6.0394 ± 0.0004</td><td>270 ± 47</td><td>0.06 ± 0.01</td><td>2.42 × 2.08</td><td>0.069 ± 0.012</td><td>VLA</td><td>S19</td></tr><tr><td>CO(6-5)</td><td>6.0379 ± 0.0022</td><td>360 ± 110</td><td>0.34 ± 0.07</td><td>5.91 × 5.12</td><td>1.17 ± 0.24</td><td>PdBI</td><td>W10</td></tr><tr><td>CO(6-5)</td><td>6.0391 ± 0.0002</td><td>229 ± 20</td><td>0.29 ± 0.05</td><td>0.39 × 0.30</td><td>1.00 ± 0.20</td><td>ALMA</td><td>T24a</td></tr><tr><td>CO(7-6)</td><td>6.0397 ± 0.0003</td><td>487 ± 115</td><td>0.235 ± 0.08</td><td>4.0 × 2.6</td><td>0.95 ± 0.38</td><td>NOEMA</td><td>D22</td></tr><tr><td>CO(8-7)</td><td>6.0396 ± 0.0005</td><td>246 ± 29</td><td>0.36 ± 0.07</td><td>3.78 × 3.01</td><td>1.66 ± 0.32</td><td>NOEMA</td><td>this work</td></tr><tr><td>CO(9-8)</td><td>6.0380 ± 0.0003</td><td>202 ± 32</td><td>0.27 ± 0.04</td><td>2.91 × 1.57</td><td>1.40 ± 0.21</td><td>NOEMA</td><td>this work</td></tr><tr><td>CO(10-9)</td><td>6.0391 ± 0.0007</td><td>233 ± 34</td><td>0.24 ± 0.07</td><td>2.41 × 1.06</td><td>1.38 ± 0.40</td><td>NOEMA</td><td>this work</td></tr><tr><td colspan="8"><strong>J0129-0035</strong></td></tr><tr><td>[C II]</td><td>5.7788 ± 0.0001</td><td>206 ± 9</td><td>2.15 ± 0.08</td><td>0.22 × 0.16</td><td>19.2 ± 0.7</td><td>ALMA</td><td>V20</td></tr><tr><td>CO(2-1)</td><td>5.7783 ± 0.0004</td><td>195 ± 41</td><td>0.036 ± 0.005</td><td>0.69 × 0.64</td><td>0.039 ± 0.005</td><td>VLA</td><td>S19</td></tr><tr><td>CO(6-5)</td><td>5.7794 ± 0.0008</td><td>283 ± 87</td><td>0.37 ± 0.07</td><td>~5 × 5</td><td>1.20 ± 0.23</td><td>PdBI</td><td>W11</td></tr><tr><td>CO(6-5)</td><td>5.7785 ± 0.0002</td><td>183 ± 16</td><td>0.18 ± 0.03</td><td>0.45 × 0.34</td><td>0.59 ± 0.09</td><td>ALMA</td><td>T24b</td></tr><tr><td>CO(8-7)</td><td>5.7800 ± 0.0014</td><td>283 ± 187</td><td>0.21 ± 0.09</td><td>4.73 × 2.46</td><td>0.91 ± 0.39</td><td>NOEMA</td><td>this work</td></tr><tr><td>CO(9-8)</td><td>5.7780 ± 0.0005</td><td>205 ± 53</td><td>0.21 ± 0.05</td><td>3.43 × 2.54</td><td>1.02 ± 0.24</td><td>NOEMA</td><td>this work</td></tr><tr><td>CO(10-9)</td><td>...</td><td>...</td><td>&lt;0.26</td><td>3.23 × 2.52</td><td>...</td><td>NOEMA</td><td>this work</td></tr></tbody></table>

**Notes.** Column (1): line ID. Columns (2)–(4): redshift, line width in FWHM, and line flux. Note that the line flux is calculated through a single Gaussian fit to the line profile. Column (5): beam size in FWHM. Column (6): line luminosity and calibration uncertainties are not included in the error bars. Column (7): facilities. Column (8): references.

**References:** R. Wang et al. (2010, W10); R. Wang et al. (2011, W11); Y. Shao et al. (2019, S19); B. P. Venemans et al. (2020, V20); R. Decarli et al. (2022, D22); R. Tripodi et al. (2024, T24a); R. Tripodi (2024, T24b).

183 km s⁻¹ as a prior constraint and employing the MCMC method in the fitting. The 3σ upper limit for the CO (10–9) line is calculated as 3σ<sub>channel</sub>(Δv<sub>channel</sub>Δv<sub>line</sub>)¹/²(in Jy km s⁻¹), where Δv<sub>line</sub> is adopted as the velocity range 275 km s⁻¹ from the ALMA CO (6–5) measurement, Δv<sub>channel</sub> is the binned channel width of 74 km s⁻¹, and σ<sub>channel</sub> is the corresponding point-source rms noise value in Jy, respectively (E. R. Seaquist et al. 1995; R. Wang et al. 2010). This estimates the upper limit to be <0.26 Jy km s⁻¹. For completeness, we also performed an MCMC fit to the CO(9–8) line and confirmed that the result agrees with the least-square fit, falling within the margin of error.

In Figure 1, the peak positions of the emission lines are marked with black crosses, all of which match the quasar optical position (01h29m58.5150 ± 0.0347, −00°35′39″810 ± 0.0′521, Y. Shao et al. 2019). The CO(8–7) and CO(9–8) emission lines show consistent line widths with those reported for CO(6–5) (R. Tripodi 2024) and [C II]158 µm (B. P. Venemans et al. 2020).

# 4. Analysis and Discussion

## 4.1. Continuum Emission

To estimate the total infrared luminosity, $L_{\text{TIR}}$, and the underlying SFR of the host galaxies, we present our SED fitting results of the two objects in this section. Here, $L_{\text{TIR}}$ refers to the luminosity integrated from 8 to 1000 µm, while $L_{\text{FIR}}$ is integrated over the range from 42.5 to 122.5 µm (G. Helou et al. 1985).

Recent studies suggest that the optically thin assumption may not be valid at FIR wavelength in compact and dusty starburst systems (J. M. Simpson et al. 2017; C. M. Casey et al. 2019). Therefore, we adopt the general dust opacity modified blackbody model (Gτ-MBB) in the continuum SED fitting to determine the dust temperature and dust mass (E. da Cunha et al. 2021).

The dust heating from the cosmic microwave background (CMB) at the quasar redshift with a temperature of $T_{\text{CMB}}(z)$ is also taken into account (E. da Cunha et al. 2013). The SED of the thermal dust emission is modeled as:

$$S_{\nu \text{obs}}^{\text{obs}} = \frac{\Omega}{(1+z)^3} [B_{\nu}(T_{\text{dust}}(z)) - B_{\nu}(T_{\text{CMB}}(z))](1 - e^{-\tau_{\nu}}), \quad (1)$$

where $\Omega = (1+z)^4AD_L^{-2}$ is the solid angle covered by the source with A, and $D_L$ is the surface area and luminosity distance of the source, respectively. $B_\nu(T)$ is the Planck function, $T_{\text{dust}}$ is the dust temperature, $T_{\text{CMB}}(z)$ is the CMB temperature at redshift z, and $\tau_\nu$ is the dust optical depth.

The dust optical depth $\tau_\nu$ can be expressed as

$$\tau_{\nu} = \kappa_{\nu} \Sigma_{\text{d}} = \frac{M_{\text{dust}}}{A} \kappa_{0} \left(\frac{\nu}{\nu_{0}}\right)^{\beta}, \quad (2)$$

where β is the emissivity index, $M_{\text{dust}}$ is the dust mass, and $\kappa_\nu = \kappa_0 \left(\frac{\nu}{\nu_0}\right)^\beta$ is the frequency-dependent dust opacity. In this paper, we adopt $\kappa_0 = 0.77 \text{ cm}^2 \text{ g}^{-1}$ at $\nu_0 = 353 \text{ GHz}$ (i.e., $\lambda_0 = 850 \text{ µm}$), which is consistent with E. da Cunha et al. (2008, 2015, 2021).

While star formation in the quasar host galaxy dominates the FIR dust emission, the central AGN may also contribute to the dust heating to a large scale (~1 kpc). For J2054-0005, the Herschel data measures the continuum emission in the rest-frame MIR-to-FIR bands. To subtract the contribution from the AGN, we utilize a simple SED fitting method to fit the NIR/MIR and FIR multiwavelength flux density measurements using the Python lmfit package. This method adopted a clump torus model (CAT3D) model from S. F. Hönig &
<table><caption>**Table 3**
Dust Continuum Observations</caption><thead><tr><th>Wavelength<br>(1)</th><th>Telescope<br>(2)</th><th>J2054-0005<br>(3)</th><th>J0129-0035<br>(4)</th></tr></thead><tbody><tr><td><i>i</i> (mag)</td><td>SDSS</td><td>23.30 ± 0.22 (a)</td><td>24.52 ± 0.25 (b)</td></tr><tr><td><i>z</i> (mag)</td><td>SDSS</td><td>20.72 ± 0.09 (a)</td><td>22.16 ± 0.11 (b)</td></tr><tr><td><i>J</i> (mag)</td><td>SDSS</td><td>19.18 ± 0.06 (a)</td><td>21.78 ± 0.15 (b)</td></tr><tr><td>F<sub>100 µm</sub> (mJy)</td><td>Herschel/PACS</td><td>&lt;2.7 (c)</td><td>...</td></tr><tr><td>F<sub>160 µm</sub> (mJy)</td><td>Herschel/PACS</td><td>9.8 ± 1.3 (c)</td><td>...</td></tr><tr><td>F<sub>250 µm</sub> (mJy)</td><td>Herschel/<br>SPIRE</td><td>15.2 ± 5.4 (c)</td><td>&lt;12.2 (d)</td></tr><tr><td>F<sub>350 µm</sub> (mJy)</td><td>Herschel/<br>SPIRE</td><td>12.0 ± 4.9 (c)</td><td>&lt;11.4 (d)</td></tr><tr><td>F<sub>444 µm</sub> (mJy)</td><td>ALMA</td><td>9.87 ± 0.94 (l)</td><td>...</td></tr><tr><td>F<sub>500 µm</sub> (mJy)</td><td>Herschel/<br>SPIRE</td><td>&lt;19.5 (c)</td><td>&lt;15.2 (d)</td></tr><tr><td>F<sub>488 GHz</sub> (mJy)</td><td>ALMA</td><td>10.35 ± 0.15 (h)</td><td>...</td></tr><tr><td>F<sub>346 GHz</sub> (mJy)</td><td>ALMA</td><td>5.723 ± 0.009 (j)</td><td>...</td></tr><tr><td>F<sub>280 GHz</sub> (mJy)</td><td>ALMA</td><td>...</td><td>2.61 ± 0.06 (i)</td></tr><tr><td>F<sub>270 GHz</sub> (mJy)</td><td>ALMA</td><td>3.15 ± 0.10 (i)</td><td>...</td></tr><tr><td>F<sub>262 GHz</sub> (mJy)</td><td>ALMA</td><td>2.98 ± 0.05 (e)</td><td>...</td></tr><tr><td>F<sub>287 GHz</sub> (mJy)</td><td>ALMA</td><td>...</td><td>2.57 ± 0.06 (e)</td></tr><tr><td>F<sub>250 GHz</sub> (mJy)</td><td>IRAM</td><td>2.38 ± 0.53 (f)</td><td>2.37 ± 0.49 (g)</td></tr><tr><td>F<sub>166.5 GHz</sub> (µJy)</td><td>NOEMA</td><td>...</td><td>602 ± 35</td></tr><tr><td>F<sub>161.0 GHz</sub> (µJy)</td><td>NOEMA</td><td>542 ± 24</td><td>...</td></tr><tr><td>F<sub>153.0 GHz</sub> (µJy)</td><td>NOEMA</td><td>...</td><td>494 ± 37</td></tr><tr><td>F<sub>147.7 GHz</sub> (µJy)</td><td>NOEMA</td><td>438 ± 23</td><td>...</td></tr><tr><td>F<sub>134.7 GHz</sub> (µJy)</td><td>NOEMA</td><td>...</td><td>285 ± 21</td></tr><tr><td>F<sub>133.3 GHz</sub> (µJy)</td><td>NOEMA</td><td>270 ± 19</td><td>...</td></tr><tr><td>F<sub>102 GHz</sub> (µJy)</td><td>IRAM</td><td>...</td><td>140 ± 40 (g)</td></tr><tr><td>F<sub>92.3 GHz</sub> (µJy)</td><td>ALMA</td><td>82 ± 9 (k)</td><td>...</td></tr></tbody></table>

**References.** (a) L. Jiang et al. (2008); (b) L. Jiang et al. (2009); (c) C. Leipski et al. (2014); (d) Y. Shao et al. (2019); (e) R. Wang et al. (2013); (f) R. Wang et al. (2008); (g) R. Wang et al. (2011); (h) T. Hashimoto et al. (2019); (i) B. P. Venemans et al. (2020); (j) D. Salak et al. (2024); (k) R. Tripodi et al. (2024); (l) R. Tripodi (2024).

M. Kishimoto (2010), which models the NIR and MIR contributions from the AGN dust torus.

For J2054-0005, we approached the SED fitting in two steps. First, we separately fit the SED by combining the G7-MBB model with each torus template from the CAT3D model. The CAT3D model includes a range of dust torus SED templates generated using various parameters, such as the inclination angle and opening angle. In this initial step, we freely fit the G7-MBB parameters $T_{\text{dust}}$ and $\beta$ to achieve the best fit for the SED, using the dust continuum size ($0.9 \times 0.7 \text{ kpc}^2$) measured in B. P. Venemans et al. (2020) as a nonfitting parameter. We identified the optimal torus template by choosing the one with the smallest $\chi^2$ value. In the second step, we refined the SED fitting by incorporating the chosen torus template alongside the G7-MBB model to extract the best-fitting parameters. This method, despite its inherent degeneracy and uncertainty, was aimed primarily at quantifying the AGN's relative contribution without needing to fit the torus' specific parameters precisely. Therefore, our free parameters included the torus component scale factor, $k_t$, the torus SED shape determined by the selected template, along with the dust temperature $T_{\text{dust}}$, dust mass $M_{\text{dust}}$, and emissivity index $\beta$.

Figure 3 displays the best-fit SED for J2054-0005, where the G7-MBB model yielded parameters of $\beta = 2.17 \pm 0.09$, a dust temperature $T_{\text{dust}} = 52 \pm 1$ K, and a dust mass $M_{\text{dust}} =$

$(1.3 \pm 0.2) \times 10^8 M_{\odot}$. T. Hashimoto et al. (2019) also fitted the SED of J2054-0005, using only four flux density measurements and assuming thin optical depth. They reported a dust temperature $T_{\text{dust}} = 50 \pm 2$ K and a dust emissivity index $\beta = 1.8 \pm 0.1$. While their $T_{\text{dust}}$ is similar to ours, they found a lower $\beta$, which can be attributed to their assumption of thin optical depth and the limited number of data points used in their fit. Our calculated total infrared luminosity is $L_{\text{TIR}} = 6.0 \pm 1.1 \times 10^{12} L_{\odot}$, which does not include the torus component heated by AGN. Adopting a Kroupa initial mass function and the scaling relation $\text{SFR}_{\text{IR}}/(M_{\odot} \text{yr}^{-1}) = 3.88 \times 10^{-44} L_{\text{TIR}}/(\text{erg s}^{-1})$ from E. J. Murphy et al. (2011), we estimated the SFR to be $896 \pm 163 M_{\odot} \text{yr}^{-1}$, a value approximately half of that reported by T. Hashimoto et al. (2019). This discrepancy primarily arises from our inclusion of the AGN torus in the SED fitting, which T. Hashimoto et al. (2019) omitted. Additionally, as F. Duras et al. (2017) suggested, about 50% of the total IR luminosity in AGNs with $L_{\text{bol}} > 10^{47} \text{ erg s}^{-1}$ is attributed to dust heated by QSOs, a finding consistent with our results.

For J0129-0035, due to the lack of Herschel data at the dust emission peak, $T_{\text{dust}}$ has significant uncertainty. We have to assume a SED shape to compute the total infrared (TIR) luminosity. We adopted the graybody with $\beta = 1.6$ and $T_{\text{dust}} = 47$ K that was often assumed in the previous optically thin model (A. Beelen et al. 2006; B. P. Venemans et al. 2020). We use the size of the dust continuum ($1.2 \times 1.0 \text{ kpc}^2$) measured in B. P. Venemans et al. (2020) as the galaxy size. Therefore, during our fitting procedure, the free parameter to fit is the dust mass $M_{\text{dust}}$. With this assumption, we measured the TIR luminosity to be $L_{\text{TIR}} = (4.9 \pm 0.4) \times 10^{12} L_{\odot}$, and the dust mass to be $M_{\text{dust}} = (3.9 \pm 0.1) \times 10^8 M_{\odot}$. The SFR estimated from $L_{\text{TIR}}$ is $733 \pm 59 M_{\odot} \text{yr}^{-1}$. However, due to the unknown AGN contribution, we considered this an upper limit for the real SFR.

To address the uncertainty in TIR luminosity and $M_{\text{dust}}$ introduced by different assumptions of parameters, we also show the SED models adopting different values of $\beta$ and $T_{\text{dust}}$ in Figure 4. Among these two parameters, the different assumptions of $T_{\text{dust}}$ introduce the largest discrepancies in $M_{\text{dust}}$ and $L_{\text{TIR}}$. In Figure 4(a), with $T_{\text{dust}}$ fixed to 40 K, we obtained a dust mass of $M_{\text{dust}} = (6.1 \pm 0.2) \times 10^8 M_{\odot}$, while with $T_{\text{dust}} = 70$ K, the dust mass is reduced by a factor of 3.7, and $L_{\text{TIR}}$ increase by a factor of 7.6. In case we fixed $T_{\text{dust}}$ to 40 K and changed $\beta$ from 1.5 to 2.2 (Figure 4(b)), $M_{\text{dust}}$ will decrease by a factor of 3.0 while $L_{\text{TIR}}$ moderately increases by a factor of 1.3.

## 4.2. Radiative Transfer Analysis of the CO SLED

We further investigate the physical conditions and excitation mechanism of the molecular gas in the quasar host galaxies based on the detections of the CO transitions. We here assume that radiative processes dominate the ISM excitation over nonradiative mechanisms and compare our observed CO SLEDs of J2054-0005 and J0129-0035 with the predictions from PDR and XDR models. Our analysis employs the CLOUDY spectral synthesis code (version c.17.03, G. J. Ferland et al. 2017). Following the formula described in A. Pensabene et al. (2021), we have developed an expanded parameter grid comprising 418 × 2 models that encompass both PDR and XDR scenarios. These models span a range of the total hydrogen density, $\log(n_H [\text{cm}^{-3}])$, from 2 to 7 (in 19 models with steps of approximately 0.29 dex), and incident radiation field strength ($\log G_0$) from 1 to 7 for PDRs and
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Spectra of CO(8–7), (9–8), (10–9) (from top to bottom) of two objects (left: J2054-0005; right: J0129-0035) and the fitting results by one Gaussian mode. The dotted histograms mark the ±1σ range in the rms noise. The red lines show the best-fit model of the spectra. The yellow bins represent the velocity range within 1.5 × FWHM of the detected line.</figcaption></figure>

log $F_X$[erg s⁻¹ cm⁻²] from −2.4 to 2.4 for XDRs. We focus on a specific case with a total hydrogen column density of log($N_H$/cm⁻²) = 23 to fully sample the molecular component typically located at $N_H > 10^{22}$ cm⁻², as observed in giant molecular clouds (C. F. McKee & J. P. Ostriker 1977; A. Pensabene et al. 2021). We assume that all line emission originates from the neutral or molecular gas phase, disregarding

the contribution of H II regions to the [C II]158 µm emission (A. Pensabene et al. 2021; R. Decarli et al. 2022). From the CLOUDY simulations, we obtained intensities of the lines of interest and the thermal dust continuum, including CO $J=1-0$, [C II]158 µm, [C I]369 µm, and, $L_{\text{TIR}}$.

We then fit the grids described above to the observed CO SLEDs of the two objects by minimizing the χ² statistic over
<figure><img src="image_5.png" /><figcaption>**Figure 3.** SED fitting results of J2054-0005. Open circles and errors are broadband photometric measurements from Herschel observation (black, C. Leipski et al. 2014), ALMA observation (orange, T. Hashimoto et al. 2019; B. P. Venemans et al. 2020; D. Salak et al. 2024; R. Tripodi et al. 2024), and NOEMA/PdBI observation (blue, R. Wang et al. 2008, 2013). The red-orange dashed line represents the best-fitted CAT3D model and the cyan dashed line shows the best-fitted host galaxy cold dust component.</figcaption></figure>

the grids of log G₀ and log (nH/cm⁻³):

$$ \chi^2 = \sum_{i=1}^{n-m} z_i + \sum_{j=1}^{m} \left( -2 \log \left( \frac{1 + \operatorname{erf}(z_j)}{2} \right) \right). \quad (3) $$

Here, n represents the measurements of CO lines, including both detections and upper limits, and m denotes the number of upper limits. The quantity zᵢ is defined as zᵢ = [yᵢ - f(xᵢ; nH, G₀)]² / σᵢ² (M. Sawicki 2012; M. Boquien et al. 2019). The quantities yᵢ and σᵢ are the flux and the corresponding uncertainty, respectively. The function f(xᵢ; nH, G₀) represents the model prediction for the line fluxes at each observed transition, and erf(zᵢ) is the error function. The modified χ² formulation, which includes the error function erf(zᵢ), allows us to incorporate upper-limit information into the analysis. In particular, for detected lines, the standard χ² term is used. For upper limits, the contribution to χ² is weighted by the error function, which effectively penalizes the model less for exceeding upper limits than for deviating from actual detections. This approach helps constrain the model parameters by using all available data, including nondetections.

### 4.2.1. Single-PDR Model

We first utilize a single PDR component to fit the CO SLEDs for both targets. The results are shown in Figures 5(a) and (d). In subsequent fittings of the CO SLED, we normalize the intensities to CO(6–5). Therefore, the fitting parameters include the total hydrogen density of the PDR, denoted as log (nH,PDR), and the radiation field strength log G₀ (where G₀ is in Habing flux units, i.e., 1.6 × 10⁻³ erg s⁻¹ cm⁻²; see H. J. Habing 1968). We found that the CO SLEDs of both

objects can be fitted with a dense gas component and intense FUV radiation field. The best-fitting parameters are log (nH/cm⁻³) = 6, log G₀ = 6.5 for J2054-0005 (Figure 5(a)), and log (nH/cm⁻³) = 6.0, log G₀ = 6.0 (Figure 5(d)) for J0129-0035.

### 4.2.2. Double-PDR Model

As the single-PDR model’s extreme physical conditions may oversimplify the actual conditions in PDRs, we now explore a scenario with double PDR components for the CO SLEDs. One PDR component is closer to the typical conditions of the extended molecular ISM and contributes mainly to low-J CO lines. The other PDR component could be denser and give rise to the bright high-J CO lines. As demonstrated by the results of the single PDR, very dense and highly irradiated PDR clumps could likely play a significant role, indistinguishable from XDR regions. We fit the SLEDs of both targets as shown in Figures 5(b) and (e) using a combined low-density PDR (PDR I) + dense PDR (PDR II) model. We normalize the CO SLED to CO(6–5), defining the contributions of two components to CO(6–5) luminosity as kPDR I and kPDR II, with kPDR I + kPDR II = 1. The fitting parameters include the total hydrogen density (log (nH,PDR I)) and the radiation field strength (log G₀,PDR I) of the PDR I, the total hydrogen density of the PDR II (log (nH,PDR II)) and its radiation field strength (log G₀,PDR II), as well as the normalization parameter for the PDR I component (kPDR I). For PDR I, we fix the gas density log (nH/cm⁻³) = 3.4 and the incident FUV field log G₀ = 1.9 in the grids as closely described in the local AGNs (log (nH/cm⁻³) ~ 3.5, log G₀ ~ 2; R. Meijerink et al. 2007; R. Meijerink 2011; L. Spinoglio et al. 2012). These specific values are used instead of the exact literature values due to the discrete nature of our grids. They are sufficiently representative of environments from the extended region encompassing the star-forming ring. By adopting this, we assume there is a similar common star-forming ISM component in these quasar host galaxies. The more dense PDR II component’s parameters are fitted without constraints. For J2054-0005, the results in Figure 5(b) reproduce the observation of CO SLED with log (nH,PDR II/cm⁻³) = 6.5 and log G₀,PDR II = 7.0. For J0129-0035, the fitting results in Figure 5(e) match the observed CO SLED with log (nH,PDR II/cm⁻³) = 6.5 and log G₀,PDR II = 7.0. The fitting results found that PDR I contributes 56% and 53% of the luminosity of the CO (2–1) line for J2054-0005 and J0129-0035, respectively. This indicates that PDR I dominates the CO (2–1) line flux and molecular gas mass, assuming the same conversion factor αCO ≈ 0.8 M☉ (K km s⁻¹ pc²)⁻¹ (C. Carilli & F. Walter 2013). The fitting results of PDR II for both quasars suggest high gas densities with very strong FUV radiation fields, which are quite close to those from the single-PDR models for both objects. This indicates that the best model points to extreme conditions in PDRs, making the contribution of the PDR I negligible for the mid- and high-J CO transitions.

### 4.2.3. PDR+XDR Model

We also investigate the possibility of whether an XDR component is involved. We define the contributions of two components to CO(6–5) luminosity as kPDR and kXDR, with kPDR + kXDR = 1. The fitting parameters include the total hydrogen density of the PDR (log (nH,PDR)) and the radiation
<figure><img src="image_6.png" /><figcaption>**Figure 4.** SED fitting results of J0129-0035. The left panel displays the variation of the fitting results with respect to different $T_{\text{dust}}$ and fixed $\beta = 1.6$. The right panel shows the variation of the fitting results for different $\beta$ and fixed $T_{\text{dust}} = 47$ K.</figcaption></figure>

field strength (log $G_0$), the total hydrogen density of the XDR (log ($n_{\text{H,XDR}}$)) and its radiation field strength (log ($F_X$/erg s⁻¹ cm⁻²)), as well as the normalization parameter for the PDR component ($k_{\text{PDR}}$). However, the observational data are insufficient to constrain these five parameters of the two components simultaneously. In the double-PDR model (Section 4.2.2), since both components contribute to the observed FIR dust continuum and fine-structure lines, determining their contributions is challenging. As a result, we fixed PDR I based on literature values. For the PDR+XDR model, we assume that the FIR dust continuum and fine-structure lines primarily originate from the star-forming PDR component. This allows us to use the observational line and FIR luminosities to constrain the PDR’s physical conditions. Given the limited data and model grid constraints, we fixed the PDR parameters derived from these constraints, instead of fitting a narrower range, to avoid boundary issues.

In Figure 6 (a), we utilize the predicted ratio of the [C II]158 µm, [C I]369 µm, and TIR luminosity to estimate the physical conditions of PDR in J2054-0005 (A. Pensabene et al. 2021). Our PDR simulations using CLOUDY also calculate the theoretical ratios between line and continuum intensities. We adopt the measurements of [C II]158 µm from B. P. Venemans et al. (2020) and $L_{\text{TIR}}$ from Section 4.1, as well as the 3σ upper limit of [C I]369 µm from R. Decarli et al. (2022). Although these ratios offer poor constraints to the hydrogen density, they effectively constrain the intensity of the radiation field, indicating log $G_0 \lesssim 4.5$ for [C II]158 µm/$L_{\text{TIR}}$ and a more stringent limit of $2.5 \lesssim \log G_0 \lesssim 4$ for [C II]369 µm/$L_{\text{TIR}}$ and [C II]158 µm/[C I]369 µm. The lower bound of [C II]158 µm/[C I]369 µm becomes more sensitive to the density at lower

$n_{\text{H}}$ values and increasingly sensitive to log $G_0$ at log ($n_{\text{H}}$/cm⁻³) $\gtrsim 4$, suggesting the host galaxy resides in a dense and intensely irradiated environment. Despite large uncertainties and its general insensitivity to the gas density, our diagnostic suggests a moderate radiation field strength, explicitly excluding scenarios with log $G_0 \gtrsim 4.5$. Consequently, our model adopts PDR conditions of log ($n_{\text{H}}$/cm⁻³) = 3.4 and log $G_0$ = 3.9. This assumption is based on the intersection of different constraints from Figure 6(a), offering a reasonable approximation of the PDR conditions in J2054-0005.

For J0129-0035, we constrain the PDR conditions using the luminosity ratios of [C II]158 µm and the $L_{\text{TIR}}$. The TIR luminosity in this object has very large uncertainties, as discussed in the previous section. [C II]158 µm has been estimated to arise predominantly from the neutral medium in $z \sim 6$ QSO hosts (M. Novak et al. 2019; A. Pensabene et al. 2021). We here calculate the upper and lower limits of [C II]158 µm to TIR ratio adopting the $L_{\text{TIR}}$ values at $T_{\text{dust}} = 40$ K and 70 K, which is shown in the gray area in Figure 6(b). These constrain the PDR conditions in J0129-0035 to be log $G_0 \lesssim 5.1$. With a similar dust temperature value described in Section 4.1, the [C II]158 µm/$L_{\text{TIR}}$ ratio of J0129-0035 is close to that of J2054-0005. Considering this, we also fix the PDR condition of J0129-0035 to log ($n_{\text{H}}$/cm⁻³) = 3.4 and log $G_0$ = 3.9 in the PDR+XDR model fitting.

With the constraints mentioned above, we fit the CO SLED of both targets as shown in Figures 5(c) and (f) using a combined PDR+XDR model. For J2054-0005, the results in Figure 5(c) reproduce the observation of CO SLED with log ($n_{\text{H,PDR}}$/cm⁻³) = 3.4, log $G_0$ = 3.9, log ($n_{\text{H,XDR}}$/cm⁻³) = 4.9, and log ($F_X$/erg s⁻¹ cm⁻²) = 0.1. It exhibits a pattern where
<figure><img src="image_7.png" /><figcaption>**Figure 5.** CO SLED fitted with PDR and/or XDR models. The black points with error bars represent the CO fluxes and the corresponding uncertainties of J2054-0005 and J0129-0035. The data and/or upper limits for CO(6–5), (8–7), (9–8), and (10–9) are from this work, with CO(6–5) obtained from ALMA observations (R. Tripodi et al. 2024), and (8–7), (9–8), and (10–9) from NOEMA observations. CO(2–1) data are sourced from Y. Shao et al. (2019). We determined the optimal parameter combination for the SLED of two sources by fitting the minimum $\chi^2$ of the residuals. For J0129-0035, we considered the impact of upper limits on the fitting by incorporating an error function into the residuals. In the upper panel, from left to right, the fits include a single-PDR model of J2054-0005, low-density PDR+dense PDR models, and PDR+XDR modeling with the line constraints ($\log(n_{\text{H}}/\text{cm}^{-3}) = 3.4$, $\log G_0 = 3.9$). The lower panel presents similar cases for J0129-0035, with the line constraints of $\log(n_{\text{H}}/\text{cm}^{-3}) = 3.4$, $\log G_0 = 3.9$. Blue lines represent the PDR component, red lines indicate the XDR component and gray lines depict the total contribution from both components.</figcaption></figure>

the excitation at low-J CO transitions is dominated by the PDR component, while the mid to high-J excitations are predominantly driven by the XDR. To estimate the uncertainty of the best-fitting parameters, we employed a Bayesian MCMC fitting algorithm, incorporating Gaussian distribution priors of ($\log(n_{\text{H,PDR}}/\text{cm}^{-3}) = 3.1$, $\log G_0 = 3.9$, $k_{\text{PDR}} = 0.32$, $\log(n_{\text{H,XDR}}/\text{cm}^{-3}) = 4.9$ and $\log(F_X/\text{erg s}^{-1}\text{cm}^{-2}) = 0.10$, with $\sigma = 0.1$ for $\log(n_{\text{H,PDR}}/\text{cm}^{-3})$ and $\log G_0$, while $\sigma = 0.5$ for the others), as predicted by the constrained minimum $\chi^2$ fitting results using line ratio constraints. The best-fitting results are displayed in Figure 7(a), yielding a derived gas density of $\log(n_{\text{H}}/\text{cm}^{-3}) = 3.39_{-0.31}^{+0.31}$ and radiation field strength $\log G_0 = 3.88_{-0.27}^{+0.34}$, as well as an XDR model with a gas density of $\log(n_{\text{H}}/\text{cm}^{-3}) = 4.93_{-0.61}^{+0.67}$ and radiation field strength $\log(F_X/\text{erg s}^{-1}\text{cm}^{-2}) = 0.11_{-0.13}^{+0.72}$.

For J0129-0035, as shown in Figure 5(f), the model also successfully reproduces the observed data with $\log(n_{\text{H,PDR}}/\text{cm}^{-3}) = 3.4$, $\log G_0 = 3.9$, $\log(n_{\text{H,XDR}}/\text{cm}^{-3}) = 4.9$ and $\log(F_X/\text{erg s}^{-1}\text{cm}^{-2}) = 0.1$. We employed a similar Bayesian MCMC fitting algorithm, taking into account the upper limit of the CO (10–9) line, to fit the combined PDR+XDR model with a Gaussian distribution prior of ($\log(n_{\text{H,PDR}}/\text{cm}^{-3}) = 3.4$, $\log G_0 = 3.9$, $k_{\text{PDR}} = 0.36$, $\log(n_{\text{H,XDR}}/\text{cm}^{-3}) = 4.9$ and $\log(F_X/\text{erg s}^{-1}\text{cm}^{-2}) = 0.1$, with $\sigma = 0.1$ for $\log(n_{\text{H,PDR}}/\text{cm}^{-3})$ and $\log G_0$, while $\sigma = 0.5$ for the others), as predicted by the constrained minimum $\chi^2$ fitting results using line ratio constraints. The optimal fitting results are presented in Figure 7(b), yielding a derived gas density of $\log(n_{\text{H}}/\text{cm}^{-3}) = 3.38_{-0.29}^{+0.32}$ and radiation field strength $\log G_0 = 3.86_{-0.27}^{+0.38}$, along with an XDR model featuring a gas density of $\log(n_{\text{H}}/\text{cm}^{-3}) = 4.91_{-0.60}^{+0.58}$ and

radiation field strength $\log(F_X/\text{erg s}^{-1}\text{cm}^{-2}) = 0.16_{-0.26}^{+0.72}$. These results align with those obtained from the constrained minimum $\chi^2$ fitting methods.

The best-fitting CO SLED models of the two sources allow us to estimate molecular gas masses and the relative impact of PDR versus XDR. Assuming a conversion factor $\alpha_{\text{CO}} \approx 0.8 M_{\odot}(\text{K km s}^{-1}\text{pc}^2)^{-1}$ and the molecular gas mass $M_{\text{H}_2} = \alpha_{\text{CO}}L_{\text{CO}(1-0)}$ (C. Carilli & F. Walter 2013), we estimate the total CO(1–0) luminosities to be $5.7 \pm 1.1 \times 10^6 L_{\odot}$ and $3.5 \pm 0.5 \times 10^6 L_{\odot}$, and the gas masses to be $9.2 \pm 1.8 \times 10^9 M_{\odot}$ and $5.7 \pm 0.9 \times 10^9 M_{\odot}$ for J2054-0005 and J0129-0035, respectively. We estimate the molecular gas mass of the XDR component to be about $3.9 \pm 0.8 \times 10^9 M_{\odot}$ and $2.2 \pm 0.3 \times 10^9 M_{\odot}$, contributing approximately 42% and 39% of the total molecular gas mass for J2054-0005 and J0129-0035, respectively.

For both sources, we found that the CO SLEDs can be modeled with a PDR component with a very high gas density and an intense FUV radiation field ($n_{\text{H}} \sim 10^6 \text{cm}^{-3}$, $G_0 \gtrsim 10^6$). However, the diagnostics from the fine-structure lines prefer PDR components with lower $n_{\text{H}}$ and $G_0$ values, which motivates us to consider PDR+XDR models. When we use these line ratios to constrain the PDR conditions and employ a PDR+XDR model to fit the CO SLED, it can be well-modeled by a PDR ($n_{\text{H}} \sim 10^3 \text{cm}^{-3}$, $G_0 \sim 10^4$) and an XDR ($n_{\text{H}} \sim 10^5 \text{cm}^{-3}$, $F_X \sim 1 \text{ erg s}^{-1}\text{cm}^{-2}$). The PDR gas densities are consistent with the typical conditions of the star-forming component in quasar host galaxies from local to high-z. We note that N. Ishii et al. (2024) recently reported a similar PDR component for J2054-0005, though the $n_{\text{H}}$ and $G_0$ values from their diagnostic are slightly different due to the additional constraints from the $[O\text{I}]_{63\text{\mu m}}$ line and a high FIR luminosity
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Grids of [C II]158 µm/[C I]369 µm intensity line ratio as a function of log G₀ and log (nH/cm⁻³) in the PDRs. Gray contours indicate model values. The different color lines indicate the line ratio measured in the J2054-0005. The red solid line and gray region represent the [C II]158 µm/TIR luminosity ratio reported in B. P. Venemans et al. (2020). The blue dashed line means the upper limit of [C I]369 µm/TIR luminosity ratio where 3σ upper limit of [C I]369 µm in the J2054-0005 was reported in R. Decarli et al. (2022). The green dotted line indicates the lower limit of the [C II]158 µm/[C I]369 µm luminosity ratio. (b) Grids of [C II]158 µm/TIR intensity line ratio as a function of log G₀ and log (nH/cm⁻³) in the PDRs. Gray contours indicate model values. The red solid line and the line width represent the [C II]158 µm/TIR luminosity ratio and the uncertainty reported in B. P. Venemans et al. (2020).</figcaption></figure>

from a different SED fitting method (D. Salak et al. 2024). The XDR components have a higher gas density compared to that of the PDR component, which is reasonable as the X-ray could reach the denser interior of the molecular cloud (R. Meijerink et al. 2007; L. Spinoglio et al. 2012; S. Gallerani et al. 2014; A. Pensabene et al. 2021; J. Li et al. 2024). However, we cannot fully rule out the scenario that the high-J CO lines are from a dense PDR component. Observations of higher-J CO transitions (J > 10) are required to distinguish between XDR and dense PDR models.

## 4.3. Complex CO Excitation in the Starburst Quasar Hosts at z ~ 6

The CO SLED is an essential tool for comprehending the physical conditions and molecular gas in galaxies, which is critical for studying the star formation process and AGN power feedback. To investigate the connection between CO excitation and the power source, we compared the CO SLEDs of the two quasars with those of local and high-z starburst galaxies and AGNs. Figure 8(a) displays a comparison with local and high-z galaxy samples, including local ULIRGs samples (M. J. F. Rosenberg et al. 2015), local normal+starburst samples (D. Liu et al. 2015), (sub)millimeter galaxies at z = 1.2–4.1 (C. Yang et al. 2017), and lensed (sub)millimeter galaxies at z = 2–4 (M. S. Bothwell et al. 2013). These starburst galaxies show CO excitation, peaking at J ~ 6 and dropping rapidly at higher J, indicating that the excitation is dominated by PDR. J2054-0005 and J0129-0035 exhibit higher CO excitation with a peak at J = 8, which requires PDR with extremely high radiation field strength and density or additional heating mechanisms as described in Section 4.2. Figure 8(b) presents a comparison to other z ~ 6 quasars that have

published CO SLEDs with at least four detected transitions (see references in the figure caption). These high-z quasars usually exhibit high CO excitation, with CO SLEDs peaked at J ~ 6 and bright high-J (J ~ 9) CO transitions. Moreover, we did not find a common shape for the CO SLED in these quasar-starburst systems at z ~ 6. For example, among the objects shown in Figure 8(b), the intensity ratios between the high-J and mid-J CO transitions vary in a wide range.

We further investigate the correlations between the CO line ratio r96 and the AGN 1450 Å absolute magnitude M1450, as well as the thermal dust emission powered by host star formation, for a sample of 10 z ~ 6 quasars. rJ2J1 is the brightness temperature ratio between two transitions and is often used to measure the CO excitation, which is expressed as

$$r_{J_2, J_1} = \frac{L'_{CO J_2 \rightarrow J_2-1}}{L'_{CO J_1 \rightarrow J_1-1}} = \frac{S_{CO J_2 \rightarrow J_2-1}}{S_{CO J_1 \rightarrow J_1-1}} \left( \frac{J_1}{J_2} \right)^2. \quad (4)$$

For some sources in the sample, the absence of mid-infrared observations prevents us from constraining the AGN's contribution to L_TIR using multicomponent SED fitting. However, since the contribution of the AGN torus is primarily concentrated in the 1–40 µm range in the rest frame, L_FIR is less affected by the AGN contribution compared to L_TIR. Therefore, in this section, we will adopt the L_FIR obtained from a single graybody model fitting to examine the correlation.

As shown in Figure 9, we employed weighted least squares to fit the correlation between r96 and both M1450 and L_FIR. There is no significant correlation between r96 and M1450 with a p-value much greater than 0.05. On the other hand, the relation between r96 and L_FIR is also very marginal with p = 0.071. Moreover, if we exclude the lensed object APM 08279 at z = 3.9113, which has the highest luminosity and line ratio, the
<figure><img src="image_9.png" /><figcaption>**Figure 7.** CO SLED fitting results for J2054-0005 ((a), left) and J0129-0035 ((b), right), employing priors obtained from best-fitting parameters through the least squares method with constrained PDR conditions and predicted line ratios. Blue (PDR) and red (XDR) lines represent the fitting results of the mode of the parameter posterior distribution (cyan). Black dots represent the observed data, and the gray line illustrates the combined PDR+XDR model.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** CO SLED normalized to CO(6–5) in normalized units of Jy km s⁻¹. Left column: J2054-0005 (blue squares) and J0129-0035 (orange squares) in comparison with the mean of four galaxy samples: local (U)LIRGs (M. J. F. Rosenberg et al. 2015, green); local normal + starburst galaxies (D. Liu et al. 2015, brown); z ~ 1.2–4.1 SMGs (M. S. Bothwell et al. 2013, red); and strongly lensed SMGs at z ~ 2–4 (C. Yang et al. 2017, purple). Right column: J2054-0005 (blue squares) and J0129-0035 (orange squares) compared with high-z quasars. The plotted quasars are J2310+1855 (J. Li et al. 2020, green), PSOJ215.15 (J. Li et al. 2024, red), 1429+5447 (J. Li et al. 2024, purple), J1148+5251 (F. Bertoldi et al. 2003; F. Walter et al. 2003; A. Beelen et al. 2006; D. A. Riechers et al. 2009; S. Gallerani et al. 2014, brown), J0100+2802 (F. Wang et al. 2019, magenta), J0439+1634 (J. Yang et al. 2019, gray), and APM 08279+5255 (A. Weiß et al. 2007; D. A. Riechers et al. 2009; C. M. Bradford et al. 2011, yellow).</figcaption></figure>

correlation disappears (*p* = 0.25). This suggests that the correlation is significantly influenced by the extreme case of APM 08279. The lack of correlations between the CO line ratio and both the AGN and star formation luminosities highlights the complexity of the CO excitation in the starburst AGN hosts. The high-*J* CO lines from different objects may be dominated by different excitation mechanisms (e.g., PDR versus XDR),

and even in a single object, the CO lines could be dominated by complex gas components with different densities and/or power resources. That is, the CO (9–8) transition could be contributed by both PDR and XDR, which confounds the connection between *r*<sub>96</sub> and the AGN power. Additionally, the X-ray luminosity is more suitable than *M*<sub>1450</sub> for investigating XDR. However, only four targets in the comparison sample have
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Relation between the CO brightness temperature ratio $r_{96}$ and the AGN absolute magnitude $M_{1450}$ (a) and the FIR luminosity $L_{FIR}$ (b). The colored points indicate the value for high-$z$ quasars. The hollow inverted triangle is the bold estimate of $r_{96}$ derived from $r_{106}$ for J0100+2802.</figcaption></figure>

detected X-ray luminosity: APM 08279, J0100+2802, J2310+1855, and J1429+5447 (F. Vito et al. 2019; E. Bertola et al. 2022; G. Migliori et al. 2023). A larger sample of high-$z$ quasars covered by both CO and X-ray measurements is required with further observations to improve the statistics here.

Future observations of higher-$J$ CO ($J \gtrsim 10$) transitions are needed, as the significant role of XDR can be more easily distinguished from PDR in these transitions. Additionally, diagnostics of other FIR lines such as $[C\,II]_{158\,\mu m} / [C\,I]_{369\,\mu m}$ can help differentiate between the two excitation mechanisms, with the maximum possible ratio for XDR being around 15 (e.g., B. P. Venemans et al. 2017; M. Novak et al. 2019; A. Pensabene et al. 2021; R. A. Meyer et al. 2022). High-resolution CO observations can map the spatial distribution of these lines, distinguishing between different regions such as star-forming areas and AGN-dominated zones, and thus provide insights into the physical conditions and processes driving the excitation (e.g., L. Spinoglio et al. 2012).

## 5. Summary

We have reported NOEMA observations of the CO(8–7), (9–8), and (10–9) lines and the underlying dust continuum in two $z \sim 6$ quasars. These two objects are about 1–2 mag fainter in rest-frame UV ($M_{1450} \sim -24$) compared to the previous $z \sim 6$ sample with CO SLED measurements, which allow us to explore the impact of AGN luminosity on CO excitation in the less luminous quasar population at the epoch of cosmic reionization. We also analyze the physical conditions of the dust and molecular gas through SED fitting and CO SLED modeling, respectively, and compare the CO SLED of the two quasars with local/high-$z$ starburst systems and other $z \sim 6$ quasars. We discuss the relation between CO excitation with AGN and FIR luminosity. The main results are summarized as follows.

1. We present a detailed analysis of the dust continuum of J2054-0005 at $z = 6.0389$ based on Herschel and ALMA data and derive the dust temperature $T_{\text{dust}} \sim 52$ K, the dust mass $M_{\text{dust}} \sim 1.3 \times 10^8 M_{\odot}$, and the total infrared luminosity $L_{\text{FIR}} = 6.0 \times 10^{12} L_{\odot}$. We also found that the contribution from AGN is indispensable for the SED

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

fitting in the infrared band. For J0129-0035 at $z = 5.7788$, we fit the SED using the G7-MBB model and derive $M_{\text{dust}} \sim 3.9 \times 10^8 M_{\odot}$ and $L_{\text{TIR}} = 4.9 \times 10^{12} L_{\odot}$ with a fixed $T_{\text{dust}} = 47$ K and $\beta = 1.6$.

2. We analyze the CO SLEDs of J2054-0005 and J1029-0035 with available PDR and XDR models. The best-fitting results of both single-PDR and double-PDR models require a gas component with a high gas density ($n_{\text{H}} \sim 10^6 \text{ cm}^{-3}$) and a powerful radiation field ($G_0 \gtrsim 10^6$) to explain the bright high-$J$ CO transitions. However, this is inconsistent with the physical condition derived from the diagnostics of the ratios of the fine-structure lines. Adopting the prior of PDR conditions from the line ratios, we found that the best-fitting results show a PDR component with $n_{\text{H}} \sim 10^3 \text{ cm}^{-3}$ and $G_0 \sim 10^4$ for the two targets and an XDR component with $n_{\text{H}} \sim 10^5 \text{ cm}^{-3}$ and $F_{\text{X}} \sim 1 \text{ erg s}^{-1} \text{ cm}^{-2}$. It reveals that the AGN may play a significant role in the CO SLED at $J \gtrsim 6$ even for the optically faint quasar with $M_{1450} \sim -24$. However, we still cannot completely rule out the possibility of an additional extremely dense PDR region.

3. The CO SLEDs of the two targets exhibit higher excitation compared to local/high-$z$ starburst systems and a similar excitation to most other $z \sim 6$ quasars. This result is consistent with the CO SLED fitting results and underlines the AGN’s critical role in the CO excitation.

4. By analyzing the relation between $r_{96}$ with $M_{1450}$ and $F_{\text{FIR}}$, we found there is no significant correlation between the CO line ratio $r_{96}$ neither with $M_{1450}$ nor with $L_{\text{TIR}}$, highlighting the complexity of the CO excitation. Large samples, other more sensitive gas tracers, and higher angular resolution observations are required to further explore the excitation mechanisms of the molecular gas.

The study of CO SLEDs in these two sources expands our understanding of quasars at $z \sim 6$ with $M_{1450} \gtrsim -26$. We can further explore the different heating mechanisms with future detections of CO SLEDs, especially high-$J$ CO transitions, in more $z \sim 6$ quasars. Once there are several QSO hosts with consistently modeled CO SLEDs using PDR/XDR, it will be possible to compare gas density and radiation field strength
values across different objects to identify any trends with QSO properties derived from other fine-structure lines.

## Acknowledgments

We acknowledge support from the National Natural Science Foundation of China (NSFC) with grant Nos. 12173002, 11991052, 11721303, and the National Key R & D Program of China grant No. 2016YFA0400703. This work is based on observations carried out under project number S21DD with the IRAM NOEMA Interferometer. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). This paper makes use of the following ALMA data: ADS/JAO. ALMA#2018.1.01289.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSTC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. A.P. acknowledges support from Fondazione Cariplo grant no. 2020-0902. Y.S. acknowledges support from the National Natural Science Foundation of China (NSFC) with grant No. 12303001, the Fundamental Research Funds for the Central Universities, and the Beihang University grant No. ZG216S2305.

## ORCID iDs

Fuxiang Xu https://orcid.org/0000-0003-0754-9795

Ran Wang https://orcid.org/0000-0003-4956-5742

Jianan Li https://orcid.org/0000-0002-1815-4839

Roberto Neri https://orcid.org/0000-0002-7176-4046

Antonio Pensabene https://orcid.org/0000-0001-9815-4953

Roberto Decarli https://orcid.org/0000-0002-2662-8803

Yali Shao https://orcid.org/0000-0002-1478-2598

Eduardo Bañados https://orcid.org/0000-0002-2931-7824

Pierre Cox https://orcid.org/0000-0003-2027-8221

Frank Bertoldi https://orcid.org/0000-0002-1707-1775

Chiara Feruglio https://orcid.org/0000-0002-4227-6035

Fabian Walter https://orcid.org/0000-0003-4793-7880

Bram P. Venemans https://orcid.org/0000-0001-9024-8322

Alain Omont https://orcid.org/0000-0002-4721-3922

Dominik Riechers https://orcid.org/0000-0001-9585-1462

Jeff Wagg https://orcid.org/0000-0002-1815-4839

Karl M. Menten https://orcid.org/0000-0001-6459-0669

Xiaohui Fan https://orcid.org/0000-0003-3310-0131

## References

Beelen, A., Cox, P., Benford, D. J., et al. 2006, ApJ, 642, 694

Bertola, E., Vignali, C., Lanzuisi, G., et al. 2022, A&A, 662, A98

Bertoldi, F., Cox, P., Neri, R., et al. 2003, A&A, 409, L47

Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103

Bothwell, M. S., Smail, I., Chapman, S. C., et al. 2013, MNRAS, 429, 3047

Bradford, C. M., Aguirre, J. E., Aikin, R., et al. 2009, ApJ, 705, 112

Bradford, C. M., Bolatto, A. D., Maloney, P. R., et al. 2011, ApJL, 741, L37

Bradford, C. M., Nikola, T., Stacey, G. J., et al. 2003, ApJ, 586, 891

Carilli, C., & Walter, F. 2013, ARA&A, 51, 105

Casey, C. M., Zavala, J. A., Aravena, M., et al. 2019, ApJ, 887, 55

da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595

da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13

da Cunha, E., Hodge, J. A., Casey, C. M., et al. 2021, ApJ, 919, 30

da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110

Decarli, R., Pensabene, A., Venemans, B., et al. 2022, A&A, 662, A60

Decarli, R., Walter, F., Venemans, B. P., et al. 2018, ApJ, 854, 97

Duras, F., Bongiorno, A., Piconcelli, E., et al. 2017, A&A, 604, A67

Fan, X., Bañados, E., & Simcoe, R. A. 2023, ARA&A, 61, 373

Fan, X., White, R. L., Davis, M., et al. 2000, AJ, 120, 1167

Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, RMxAA, 53, 385

Gallerani, S., Ferrara, A., Neri, R., & Maiolino, R. 2014, MNRAS, 445, 2848

Greve, T. R., Leonidaki, I., Xilouris, E. M., et al. 2014, ApJ, 794, 142

Guilloteau, S., & Lucas, R. 2000, in Imaging at Radio Through Submillimeter Wavelengths, ed. J. G. Mangum & S. J. E. Radford, Vol. 217 (San Francisco, CA: ASP), 299

Habing, H. J. 1968, BAN, 19, 421

Hashimoto, T., Inoue, A. K., Tamura, Y., et al. 2019, PASJ, 71, 109

Helou, G., Soifer, B. T., & Rowan-Robinson, M. 1985, ApJL, 298, L7

Hönig, S. F., & Kishimoto, M. 2010, A&A, 523, A27

Ishii, N., Hashimoto, T., Ferkinhoff, C., et al. 2024, arXiv:2408.09944

Jiang, L., Fan, X., Annis, J., et al. 2008, ApJ, 135, 1057

Jiang, L., Fan, X., Bian, F., et al. 2009, ApJ, 138, 305

Jiang, L., Fan, X., Vestergaard, M., et al. 2007, AJ, 134, 1150

Jiang, L., McGreer, I. D., Fan, X., et al. 2015, AJ, 149, 188

Jiang, L., McGreer, I. D., Fan, X., et al. 2016, ApJ, 833, 222

Kamenetzky, J., Rangwala, N., Glenn, J., Maloney, P. R., & Conley, A. 2016, ApJ, 829, 93

Leipski, C., Meisenheimer, K., Walter, F., et al. 2014, ApJ, 785, 154

Li, J., Wang, R., Pensabene, A., et al. 2024, ApJ, 962, 119

Li, J., Wang, R., Riechers, D., et al. 2020, ApJ, 889, 162

Liu, D., Gao, Y., Isaak, K., et al. 2015, ApJ, 810, L14

Matsuoka, Y., Iwasawa, K., Onoue, M., et al. 2022, ApJS, 259, 18

McKee, C. F., & Ostriker, J. P. 1977, ApJ, 218, 148

Meijerink, R. 2011, in IAU Symp. 280, The Molecular Universe (Cambridge: Cambridge Univ. Press), 177

Meijerink, R., Kristensen, L. E., Weiß, A., et al. 2013, ApJL, 762, L16

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

Murphy, E. J., Condon, J. J., Schinnerer, E., et al. 2011, ApJ, 737, 67

Neeleman, M., Novak, M., Venemans, B. P., et al. 2021, ApJ, 911, 141

Neeleman, M., Walter, F., Decarli, R., et al. 2023, ApJ, 958, 132

Novak, M., Bañados, E., Decarli, R., et al. 2019, ApJ, 881, 63

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

Riechers, D. A., Walter, F., Bertoldi, F., et al. 2009, ApJ, 703, 1338

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

Salak, D., Hashimoto, T., Inoue, A. K., et al. 2024, ApJ, 962, 1

Sawicki, M. 2012, PASP, 124, 1208

Seaquist, E. R., Ivison, R. J., & Hall, P. J. 1995, MNRAS, 276, 867

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

Shen, Y., Wu, J., Jiang, L., et al. 2019, ApJ, 873, 35

Simpson, J. M., Smail, I., Swinbank, A. M., et al. 2017, ApJ, 839, 58

Spinoglio, L., Pereira-Santaella, M., Busquet, G., et al. 2012, ApJ, 758, 108
Tripodi, R. 2024, PhD thesis, Universita Degli Studi di Trieste https://hdl.handle.net/11368/3069363

Tripodi, R., Feruglio, C., Fiore, F., et al. 2024, A&A, 689, A220

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., McMahon, R. G., Walter, F., et al. 2012, ApJL, 751, L25

Venemans, B. P., Neeleman, M., Walter, F., et al. 2019, ApJL, 874, L30

Venemans, B. P., Walter, F., Decarli, R., et al. 2017, ApJ, 845, 154

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

Vito, F., Brandt, W. N., Bauer, F. E., et al. 2019, A&A, 630, A118

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., Neri, R., et al. 2010, ApJ, 714, 699

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

Wang, R., Shao, Y., Carilli, C. L., et al. 2019, ApJ, 887, 40

Wang, R., Wagg, J., Carilli, C. L., et al. 2011, AJ, 142, 101

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., et al. 2005, A&A, 438, 533

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