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

# Probing the Full CO Spectral Line Energy Distribution (SLED) in the Nuclear Region of a Quasar-starburst System at $z = 6.003$

Jianan Li<sup>1,2</sup>①, Ran Wang<sup>2,17</sup>①, Dominik Riechers<sup>3,4</sup>①, Fabian Walter<sup>4</sup>①, Roberto Decarli<sup>5</sup>①, Bram P. Venamans<sup>4</sup>, Roberto Neri<sup>6</sup>①, Yali Shao<sup>1,2</sup>①, Xiaohui Fan<sup>7</sup>①, Yu Gao<sup>8,15</sup>, Chris L. Carilli<sup>9,10</sup>①, Alain Omont<sup>11</sup>①, Pierre Cox<sup>11</sup>, Karl M. Menten<sup>12</sup>①, Jeff Wagg<sup>13</sup>, Frank Bertoldi<sup>14</sup>①, and Desika Narayanan<sup>16</sup>①

<sup>1</sup>Department of Astronomy, School of Physics, Peking University, Beijing 100871, People's Republic of China; jiananl@pku.edu.cn

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

<sup>3</sup>Department of Astronomy, Cornell University, Space Sciences Building, Ithaca, NY 14853, USA

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

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

<sup>6</sup>Institute de Radioastronomie Millimetrique, St. Martin d'Heres, F-38406, France

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

<sup>8</sup>Purple Mountain Observatory & Key Laboratory for Radio Astronomy, Chinese Academy of Sciences, 10 Yuanhua Road, Nanjing 210033, People's Republic of China

<sup>9</sup>Cavendish Laboratory, 19 J. J. Thomson Avenue, Cambridge CB3 0HE, UK

<sup>10</sup>National Radio Astronomy Observatory, Socorro, NM 87801-0387, USA

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

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

<sup>13</sup>SKA Organization, Lower Withington Macclesfield, Cheshire SK11 9DL, UK

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

<sup>15</sup>Department of Astronomy, Xiamen University, Xiamen, Fujian 361005, China

<sup>16</sup>Department of Astronomy, University of Florida, 211 Bryant Space Science Center, Gainesville, FL 32611, USA

Received 2019 November 13; revised 2019 December 18; accepted 2019 December 27; published 2020 February 4

## Abstract

We report Atacama Large Millimeter/submillimeter Array observations of CO(8–7), (9–8), H₂O(2₀,₂–1₁,₁), and OH⁺(l₁–O₁) and NOrthern Extended Millimeter Array observations of CO(5–4), (6–5), (12–11), and (13–12) toward the $z = 6.003$ quasar SDSS J231038.88+185519.7, aiming to probe the physical conditions of the molecular gas content of this source. We present the best sampled CO spectral line energy distribution (SLED) at $z = 6.003$, and analyzed it with the radiative transfer code MOLPOP-CEP. Fitting the CO SLED to a one-component model indicates a kinetic temperature $T_{\text{kin}} = 228$ K, molecular gas density log(n(H₂)/cm⁻³) = 4.75, and CO column density log(N(CO)/cm⁻²) = 17.5; although, a two-component model better fits the data. In either case, the CO SLED is dominated by a “warm” and “dense” component. Compared to samples of local (Ultra) Luminous Infrared Galaxies, starburst galaxies, and high-redshift submillimeter galaxies, J2310+1855 exhibits higher CO excitation at (J ≥ 8), like other high-redshift quasars. The high CO excitation, together with the enhanced L_H₂O/L_IR, L_H₂O/L_CO, and L_OH⁺/L_H₂O ratios, suggests that besides the UV radiation from young massive stars, other mechanisms such as shocks, cosmic-rays, and X-rays might also be responsible for the heating and ionization of the molecular gas. In the nuclear region probed by the molecular emissions lines, any of these mechanisms might be present due to the powerful quasar and the starburst activity.

Unified Astronomy Thesaurus concepts: Submillimeter astronomy (1647); Quasars (1319); Starburst galaxies (1570); Infrared galaxies (790); AGN host galaxies (2017); High-redshift galaxies (734); Interstellar medium (847); Molecular gas (1073); Interstellar phases (850); Early universe (435)

## 1. Introduction

The quasars discovered at $z \gtrsim 6$ represent the first generation of super massive black holes (SMBHs) and host galaxies. Many among these earliest systems host SMBHs of ~$10^9 M_{\odot}$ (e.g., Jiang et al. 2007, 2016), and the strong dust continuum and [C II] detections reveal dynamical masses of ~$10^{10}-10^{11} M_{\odot}$ and star formation rate (SFR) of ~$10^2-10^3 M_{\odot}$ yr⁻¹ in the host galaxies (e.g., Maiolino et al. 2005; Wang et al. 2008, 2016, 2019b; Walter et al. 2009; Carilli & Walter 2013; Venemans et al. 2017b, 2019; Decarli et al. 2018; Neeleman et al. 2019). These suggest that the SMBH and galaxy coevolution is already in place in these $z \sim 6$ quasar-starburst systems. In the meantime, bright molecular CO emission lines are widely detected in starburst quasar hosts, which reveal a molecular gas content of ~$10^9-10^{10} M_{\odot}$ within a few kpc scale (e.g., Bertoldi et al. 2003;

Walter et al. 2003; Riechers et al. 2009; Wang et al. 2010, 2011a, 2013, 2016; Carilli & Walter 2013; Venemans et al. 2017a). In particular, $z \gtrsim 6$ quasars are detected in very high (rotational quantum number) J (e.g., $J \ge 9$) CO transitions, indicating high CO excitation comparable to that found in local extreme (Ultra)luminous Infrared Galaxies ((U)LIRGs) and active galactic nuclei (AGNs; e.g., Gallerani et al. 2014; Wang et al. 2019a; Yang et al. 2019). Taking advantage of the most powerful sub/millimeter and radio facilities, such as Atacama Large Millimeter/submillimeter Array (ALMA), NOrthern Extended Millimeter Array (NOEMA), and the Karl G. Jansky Very Large Array (VLA), extensive observations at sub-mm/ mm wavelengths have recently been carried out to search for the emission lines from the ionized, atomic, and molecular interstellar medium (ISM; e.g., Decarli et al. 2018; Hashimoto et al. 2018; Walter et al. 2018; Novak et al. 2019). These observations are crucial for our understanding of the physical,

<footnote n="17"><sup>17</sup> Corresponding author.</footnote>
chemical conditions and kinematics of the multiphase ISM in
these young quasar hosts at the earliest epoch and allow us to
study the coevolution of SMBHs and their host galaxies at the
earliest evolutionary phase.

CO emission lines have long been the workhorse probing
molecular gas at rest-frame sub-mm band in the local and high-
redshift universe. Low *J* (e.g., *J* ≲ 3) CO emission lines are
easy to excite in typical molecular cloud conditions (i.e., the
lowest CO transition requires only ~5 K above ground and
densities of ~100 cm⁻³); thus, it traces the bulk of the
molecular gas content. Mid *J* (e.g., 4 ≲ *J* ≲ 8) CO transitions
are found to be linearly correlated with the far-infrared
luminosity and trace the SFR (e.g., Greve et al. 2014; Liu
et al. 2015). UV photons from newly formed high-mass stars
are responsible for the molecular gas heating at this regime.
The excitation of high *J* (e.g., *J* ≳ 9) CO transitions require
both high temperature and high density, which are usually
related to processes such as shocks, X-rays from AGNs, and
cosmic-rays (Bradford et al. 2003; Spinoglio et al. 2012;
Meijerink et al. 2013; Gallerani et al. 2014).

CO spectral line energy distribution (SLED)—the CO flux as a
function of rotational quantum number *J* is a probe of the
molecular gas physical conditions (e.g., temperature, density and
illuminating radiation field strength). It has been used to study the
physical conditions in a variety of local and high-redshift systems
(e.g., Weiß et al. 2005, 2007; Riechers et al. 2006; Bradford et al.
2009, 2011; Spinoglio et al. 2012; Gallerani et al. 2014; Wang
et al. 2019a; Yang et al. 2019). Normal star-forming galaxies have
CO SLEDs that peak at relatively low *J* CO transitions. The CO
SLED of the Milky Way's inner disk peaks at around *J* = 3–4
(Fixsen et al. 1999), and similar CO SLEDs are found in other
local star-forming galaxies (Daddi et al. 2015). Galaxies that
experience higher star formation activity than normal galaxies, e.g.,
starburst galaxies, (U)LIRGs, and submillimeter galaxies (SMGs),
have moderately excited CO SLEDs that peak at higher *J*
compared to the star-forming galaxies. One of the closest examples
is the local starburst galaxy M82, whose CO SLED peaks at *J* = 5
in the central region and shows little line intensity at *J* > 9 (Weiß
et al. 2005; Panuzzo et al. 2010). The molecular CO in luminous
AGNs is usually highly excited. Very high *J* CO transitions
(*J* ≳ 9) are detected in well-known AGNs such as NGC 1068
(Spinoglio et al. 2012), Mrk 231 (van der Werf et al. 2010), the
Cloverleaf quasar at *z* = 2.56 (Bradford et al. 2009; Uzgil et al.
2016), and the *z* = 3.9 lensed quasar APM 08279+5255 (Weiß
et al. 2007; Riechers et al. 2009; Bradford et al. 2011). And known
to date, bright CO emission lines at *J* ≳ 10 were detected in the
host galaxies of three *z* ≳ 6 quasars, i.e., SDSS J114816.64
+525150.3 (hereafter J1148+5251) at *z* = 6.4 (e.g., Riechers et al.
2009; Gallerani et al. 2014), SDSS J010013.02+280225.8
(hereafter J0100+2802) at *z* = 6.3 (Wang et al. 2019a), and
UHS J043947.08+163415.7 (hereafter J0439+1634) at *z* = 6.5
(Yang et al. 2019). The CO lines at *J* ≳ 10 likely arose from warm
gas with a kinetic temperature of *T*<sub>kin</sub> ≳ 100 K. Kinetic
temperature is proportional to kinetic energy through *T*<sub>kin</sub> = (2/*E*<sub>n</sub>)/3 *k*<sub>p</sub>,
where *k*<sub>B</sub> is Boltzmann constant and *E*<sub>n</sub> is the kinetic energy of the
molecule. In the high-redshift quasars, X-rays from AGNs are
frequently proposed to explain CO excitation at high *J* transitions
(Bradford et al. 2009; Gallerani et al. 2014; Uzgil et al. 2016).

The H₂O and OH⁺ lines provide additional diagnostics of
heating and ionization sources of the molecular gas (e.g.,
Cosmic-rays, UV radiation, X-rays, and shocks) in addition to
CO. The H₂O molecule traces the warm and dense molecular

regions. It is found to be bright in infrared luminous galaxies
and can even reach luminosities comparable to CO in these
galaxies (van der Werf et al. 2011; Omont et al. 2013; Yang
et al. 2013, 2016; Jarugula et al. 2019). Recent studies found a
nearly linear relation between the water luminosity and the
infrared luminosity in local and high-redshift systems over
three orders of magnitude (González-Alfonso et al. 2010, 2014;
Omont et al. 2013; Riechers et al. 2013; Yang et al. 2013,
2017; Jarugula et al. 2019). The brightest water lines are
detected in the presence of shocks or X-rays (González-
Alfonso et al. 2010; Pellegrini et al. 2013). Accordingly, the
H₂O emission may even act as a tracer of the powering source
of molecular gas, e.g., if the molecular gas is heated by the UV
radiation or other mechanisms like shocks and X-rays. In
addition, a variety of chemical processes are enrolled in the
formation of H₂O and OH⁺. The gas-phase H₂O molecule is
formed by either solid-phase or gas-phase chemical reactions
(Pellegrini et al. 2013; Yang et al. 2016). Neutral-neutral and
ion-neutral reactions are two mechanisms to form H₂O
molecules in the gas phase. The former is usually related to
the shocks, while the latter is associated with PDRs, cosmic-
ray-dominated regions, and X-ray-dominated regions (Yang
et al. 2016). The molecular ions (e.g., H₂O⁺, OH⁺) as
intermediates of the ion-neutral reactions, play important roles
in distinguishing between shocks and PDR/XDR/cosmic-ray-
dominated region. The bright H₂O⁺ line detection in a sample
of high-redshift lensed SMGs, for example, is probably
initiated by cosmic-rays (Yang et al. 2016).

The OH⁺ line, although not as strong as the turbulent gas tracer CH⁺ that has been detected in high-redshift systems (e.g., Falgarone et al. 2017), traces the turbulent gas components as well, e.g., inflow or outflows. The formation of OH⁺ requires both atomic and molecular hydrogen, at which column density the cosmic-rays or X-rays are more capable of penetrating and ionizing the neutral and molecular gas (e.g., van der Werf et al. 2010; Meijerink et al. 2011; González-Alfonso et al. 2018). The OH⁺ line has been detected both in absorption (that probes the cold turbulent gas) and in emission, where the chemical structure is dominated by cosmic-ray ionization or X-ray radiation from AGNs (e.g., van der Werf et al. 2010; González-Alfonso et al. 2018). Limited by the weak strength and the P Cygni line profile (presence of both absorption and emission in the profile of the same spectral line), there is only one reported OH⁺ line detection at *z* ≳ 6 in absorption in the starburst galaxy HFLS3 (Riechers et al. 2013).

In order to understand the physical conditions and the
heating mechanisms of the ISM in the complex environment
with both AGN and nuclear starburst activities in these young
quasar hosts at *z* ~ 6, we here present a study of the CO SLED
in one of the most far-infrared and CO luminous quasars at
*z* ~ 6, SDSS J231038.88+185519.7 (hereafter J2310+1855)
at *z* = 6.003. J2310+1855 hosts a SMBH of ≈4 × 10⁹ M☉
(Jiang et al. 2016). It was detected in bright dust continuum,
CO(2–1) and (6–5), [C II] 158 µm, and [O III] 88 µm lines
(Wang et al. 2013; Feruglio et al. 2018; Hashimoto et al. 2018;
Odorico et al. 2018; Carniani et al. 2019; Shao et al. 2019).
This quasar is also detected in bright CO(10–9) emission with a
line flux of 1.04 ± 0.17 Jy km s⁻¹ (D. A. Riechers et al. 2020,
in preparation). The far-infrared dust emissions suggest the host
galaxy is actively forming stars with an SFR of ≈2400 M☉ yr⁻¹
and is abundant in dust with a dust mass of ≈1.7 × 10⁹ M☉
(Wang et al. 2013; Hashimoto et al. 2018; Carniani et al. 2019;
Shao et al. 2019). Wang et al. (2013), for the first time, spatially resolved the large amount of gas ($\approx 9.6 \times 10^{10} M_{\odot}$) residing in a $0.75 \times 0.40$ (~3 kpc) disk based on the [C II] 158 µm observations. The CO(2–1) emission line has a size comparable to the [C II] 158 µm emission with an associated molecular gas mass of $\approx 4.3 \times 10^{10} M_{\odot}$ (Shao et al. 2017). The [O III] 88 µm emission indicates an $L[\text{O III}]/L_{\text{IR}}$ ratio comparable to local systems at similar $L_{\text{IR}}$. The quasar is luminous in X-ray as well and has a derived $L_{2-10\,\text{kev}} = 6.93 \times 10^{44} \, \text{erg s}^{-1}$ (Vito et al. 2019). In addition, there is evidence of companions close to the quasar although further confirmation is needed (Feruglio et al. 2018; Odorico et al. 2018). All together, J2310+1855 is an extraordinary quasar-starburst sample enabling us to study in detail an SMBH and galaxy coevolution at $z \sim 6$.

In this paper, we present our new ALMA observations of CO(8–7), (9–8), $H_2O(2_{0,2}-1_{1,1})$ and $OH^+(1_1-0_1)$ and NOEMA observation of CO(5–4), (6–5), (12–11), and (13–12) toward J2310+1855, aiming to investigate the molecular gas excitation mechanisms in environments of both intense star formation activity and luminous AGNs. The paper is organized as follows: In Sections 2 and 3, we present the observations and results. In Section 4, we analyze the CO SLED with a radiative transfer code to probe the physical conditions of molecular gas. In Section 5, we compare the CO excitation in J2310+1855 with local and high-redshift galaxy samples and AGNs and discuss the heating mechanisms of molecular gas as well. Finally, we summarize the results in Section 6. We adopt a flat $\Lambda$CDM cosmology with $H_0 = 70 \, \text{km s}^{-1} \text{Mpc}^{-1}$ and $\Omega_m = 0.3$, where $1''$ corresponds to 5.7 kpc at the J2310+1855 redshift ($z = 6.0031$), and the luminosity distance to J2310+1855 is 57763 Mpc.

# 2. Observations

## 2.1. ALMA

We observed the CO(8–7) ($\nu_{\text{rest}} = 921.7997$ GHz), CO(9–8) ($\nu_{\text{rest}} = 1036.9124$ GHz), $H_2O(2_{0,2}-1_{1,1})$ ($\nu_{\text{rest}} = 987.9268$ GHz, hereafter $H_2O$) and $OH^+(1_1-0_1)$ ($\nu_{\text{rest}} = 1033$ GHz, hereafter $OH^+$) emission lines, as well as the underlying continuum toward $z = 6.003$ quasar J2310+1855 with ALMA (Cycle 3, ID 2015.1.01265.S). All of the observations were executed between April and November in 2016 with beam sizes between $0.6''$ and $0.8''$. The CO(8–7), (9–8), $H_2O$, and $OH^+$ lines were observed in ALMA Band 4 with two separate executions, where 36 to 44 12 m diameter antennas were used during observations. For each observation, we used four spectral windows, each with a width of 1.875 GHz consisting of 128 channels, with two of the windows in the lower sideband (LSB) and the other two in the upper sideband (USB). The CO(8–7) and $H_2O$ lines were observed in one spectral setup, with one spectral window centered on the CO(8–7) observed frequency of 131.6274 GHz, one window covering the $H_2O$ emission, and the other two covering line-free dust continuum. In the other turning, we observed the CO(9–8) line centered at the frequency of 148.0648 GHz with the $OH^+$ line also covered in the same spectral window, while the other three windows measured the dust continuum. The fluxes were calibrated using the standard flux calibrator *Pallas*, while SDSS J2253+1608 was used as both the phase calibrator and the bandpass calibrator. The typical calibration uncertainty is <5% in ALMA band 4; here, we use 15% uncertainty that also includes the uncertainties of the old *Pallas* flux model in early casa versions (Stanley et al. 2019).

The observational data were calibrated and reduced with the Common Astronomy Software Applications (CASA) software

package version 4.7.0 (McMullin et al. 2007), using the standard ALMA pipeline. The maps were generated using the CLEAN task in CASA, and we apply the robust weighting algorithm with a Briggs parameter of 2 equivalent to natural weighting. This results in an FWHM synthesized beam size of $0.79 \times 0.75$ and $0.77 \times 0.63$ at CO(8–7)($H_2O$) and CO(9–8)($OH^+$) observing frequency, and $0.75 \times 0.72$ and $0.80 \times 0.65$ at 136.6 and 141.1 GHz for the continuum. The continuum level was determined using a first-order polynomial, and the emission lines were imaged from the continuum subtracted data cube with all of the line-emitting channels included. With a total on-source time of 34.9 and 30.3 minutes for CO(8–7)($H_2O$) and CO(9–8)($OH^+$), we finally binned the data of CO(8–7)($H_2O$) and CO(9–8)($OH^+$) to 36 and 32 km s$^{-1}$, and the corresponding rms sensitivities were 0.17 and 0.19 mJy beam$^{-1}$, respectively. The sensitivity of the underlying continuum was $15 \, \mu\text{Jy beam}^{-1}$.

## 2.2. NOEMA

We observed CO(5–4) ($\nu_{\text{rest}} = 576.2679$ GHz), CO(6–5) ($\nu_{\text{rest}} = 691.4731$ GHz), CO(12–11) ($\nu_{\text{rest}} = 1381.9951$ GHz), CO(13–12) ($\nu_{\text{rest}} = 1496.9229$ GHz), and the underlying continuum of this quasar with NOEMA (Project W18EE). The CO(5–4) and (6–5) lines were observed in one tuning with the PolyFix correlator in Band 1 (3 mm), with the CO(5–4) line in the LSB and CO(6–5) in the USB, each with 7.744 GHz bandwidth. The observations were executed in the A configuration on 2019 January 18 with a total observing time of 2 hr, with 1.18 hr on source, while the rest of the time was expended for calibrations. A total of eight or nine antennas was used. The CO(12–11) and (13–12) lines were observed in the C/D configuration in Band 3 (1 mm) with one frequency setup, with the CO(12–11) line in the LSB and CO(13–12) in the USB, each with a 7.744 GHz bandwidth. The observations started on 2019 April 17 and ended on 2019 May 1. The total observing time was 8 hr with 6.2 hr on source, and a total number of 8–10 antennas were used in the observations. 3C454.3 was used as a phase calibrator throughout the entire CO observation. The typical calibration uncertainty is <10% in the 3 mm band and <20% in the 1 mm band.

The data were reduced with the Grenoble Image and Line Data Analysis System software (GILDAS; Guilloteau & Lucas 2000) packages CLIC and MAPPING. We extracted the continuum from all line free channels in the uv plane with UV_AVERAGE. The uv table of spectral lines was generated through UV_SUBTRACT with the underlying continuum subtracted. Both the uv table of the continuum and spectral lines were cleaned with the HOGBOM algorithm, and NATURAL weighting was used to ensure the maximum signal-to-noise ratio (S/N). This results in an FWHM synthesized beam size of $1.67 \times 1.37/1.42 \times 1.19$ for CO(5–4)/(6–5) and $2.08 \times 1.62/1.91 \times 1.53$ for CO(12–11)/(13–12). The calibrated data of CO(5–4) and (6–5) were smoothed by a factor of eight in frequency, resulting in a spectral resolution of 16 MHz (~60 km s$^{-1}$), and the calibrated CO(12–11) and (13–12) data were binned to 40 MHz (~60 km s$^{-1}$) resolution. The sensitivity was $0.32$ mJy beam$^{-1}$ for CO(5–4) and (6–5) and $0.54$ mJy beam$^{-1}$ for CO(12–11) and (13–12) per binned channels.

The observational details are listed in Table 1.
<table><caption>**Table 1**
Observational Details</caption><thead><tr><th>Line ID<br>(1)</th><th>ν<sub>obs</sub><br>(GHz)<br>(2)</th><th>δν(δν)<br>MHz (km s<sup>-1</sup>)<br>(3)</th><th>Band<br>(4)</th><th>T<sub>source</sub><br>(min)<br>(5)</th><th>T<sub>tot</sub><br>(min)<br>(6)</th><th>rms<br>(mJy beam<sup>-1</sup>)<br>(7)</th></tr></thead><tbody><tr><td>CO(5–4)</td><td>82.2875</td><td>16(58)</td><td>NOEMA 3 mm</td><td>70.8</td><td>120</td><td>0.32</td></tr><tr><td>CO(6–5)</td><td>98.7381</td><td>16(49)</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CO(8–7)</td><td>131.6274</td><td>15.625(36)</td><td>ALMA Band 4</td><td>34.9</td><td>53.3</td><td>0.17</td></tr><tr><td>H<sub>2</sub>O(2<sub>0,2</sub>−1<sub>1,1</sub>)</td><td>141.0699</td><td>15.625(36)</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CO(9–8)</td><td>148.0648</td><td>15.625(32)</td><td>ALMA Band 4</td><td>30.3</td><td>48.5</td><td>0.19</td></tr><tr><td>OH<sup>+</sup>(1<sub>1</sub>−0<sub>1</sub>)</td><td>147.5061</td><td>15.625(32)</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CO(12–11)</td><td>197.3405</td><td>40(61)</td><td>NOEMA 1 mm</td><td>372</td><td>600</td><td>0.54</td></tr><tr><td>CO(13–12)</td><td>213.7515</td><td>40(56)</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

Note. Column 1: Line ID; Column 2: Line center frequency in the observer frame; Column 3: Binned spectral resolution in frequency (velocity); Columns 4–7: Observing band, on-source time, total observing time, and achieved sensitivity per binned channel. The lines without on-source time are observed in the same frequency setup as the upper ones with on-source time listed in the table.

### 3. Results

With ALMA, we detect the CO(8–7), (9–8), and H₂O emission lines at 25σ, 22σ, and 15σ, respectively. The CO emission lines and the H₂O line are marginally resolved. We also obtained a tentative signal (4σ) for the OH⁺ line close to the CO(9–8) line. The line intensity maps integrated over the line-emitting channels are presented in Figure 1. The right panel of Figure 1 shows the line spectra integrated within the 2σ contour in the intensity map. We calculate the line widths, redshift, and fluxes of CO(8–7), (9–8), and H₂O by fitting a Gaussian profile to the spectra. We fit a 2D Gaussian component to the intensity maps of CO and H₂O lines, and the source sizes are derived by deconvolving the fitted component with the beam. The spectral profiles of CO(9–8), (8–7), and H₂O are similar (Figure 7), suggesting that the high J CO and H₂O lines are probing similar regions. The source sizes measured from the CO(9–8), (8–7), and H₂O lines of ≈(0''4 ± 0''1) × (0''3 ± 0''1) are slightly smaller than values found from the previous CO(2–1) observation of (0''6 ± 0''2) × (0''4 ± 0''2) and [C II] 158 µm observation of (0''6 ± 0''1) × (0''4 ± 0''1) at similar spatial resolution (Wang et al. 2013; Feruglio et al. 2018; Shao et al. 2019). This may imply that the CO(9–8), (8–7), and H₂O lines trace similar dense molecular regions that are closer to the central SMBH compared to CO(2–1) and [C II] 158 µm. The redshift measured with the CO(9–8), (8–7), and H₂O lines are within the uncertainties consistent with that from previous [C II], CO(2–1), and (6–5) observations (Wang et al. 2013; Feruglio et al. 2018; Shao et al. 2019). The OH⁺ line is not as strong as the other detections, so we fix the center frequency to the [C II] 158 µm redshift and fit a Gaussian profile to the spectra extracted from the peak pixel. As for the line widths, all of the ALMA detections show line widths of ~400 km s⁻¹ consistent with previous CO and [C II] observations. From the OH⁺ spectra, we find that there is an absorption-like feature in the line center frequency, but the current S/N is insufficient to confirm this feature. The continuum detections were published in Shao et al. (2019). The derived continuum source sizes of ≈(0''30 ± 0''04) × (0''22 ± 0''06) are comparable to that measured with the CO and H₂O lines.

We detected CO(5–4), (6–5), (12–11), and (13–12) with NOEMA. All four CO lines are unresolved. For the line widths, redshift, and fluxes calculation, we fit a Gaussian profile to the spectra extracted from the peak intensity pixel. The redshift measured with the CO(5–4), (6–5), (12–11), and (13–12) are

consistent with our ALMA detections as well as previous CO and [C II] detections. The line widths detected in ALMA and NOEMA observations are consistent. The CO(6–5) line has been previously detected with ALMA and the (pre-NOEMA) IRAM Plateau de Bure Interferometer (PdBI) with different spatial resolutions and the spectral line flux measured with ALMA is found to be only 70% of that found with the PdBI (Wang et al. 2013; Feruglio et al. 2018). In addition to the calibration uncertainties in different observations, it is also possible that the low resolution PdBI data include more flux from the extended region. The new NOEMA observation yield a CO(6–5) flux of 1.05 ± 0.07 Jy km s⁻¹, consistent with the results obtained with ALMA (Feruglio et al. 2018). We also detected the underlying continuum at high S/N. The measurements of redshift, line widths, line fluxes, and deconvolved source sizes of our ALMA and NOEMA observations as well as previous detections are summarized in Table 2, and the continuum measurements are listed in Table 3. The continuum, line intensity maps, and spectra of CO(5–4), (6–5), (8–7), (9–8), (12–11), (13–12), H₂O, and OH⁺ are shown in Figures 2 (NOEMA) and 1 (ALMA).

Figure 3 shows the velocity and velocity dispersion maps of [C II] 158 µm (Wang et al. 2013), CO(8–7), (9–8), and H₂O. The velocity fields of CO(8–7), (9–8), and H₂O overall follow the velocity gradient observed in [C II] 158 µm from north to south, which indicates that the emission might trace a rotating molecular gas disk. The velocity field of CO(8–7) shows a high-velocity dispersion part in the western part, which is not observed in CO(9–8) and H₂O. Such irregular velocity structure in CO(8–7) is likely to be a result of the low S/N. As for the velocity dispersion, CO(8–7), (9–8), and H₂O show velocity dispersion of <100 km s⁻¹ in the outskirts (that is not likely influenced much by the beam smearing effect). Higher S/N observations, possibly at even higher angular resolution, are required in order to constrain the kinematic structures of the dense molecular gas.

### 4. Radiative Transfer Analysis of the CO SLED

In the CO SLED analysis, we also include a new detection of the CO (10–9) line from NOEMA at high S/N. More details of the observation will be described in Riechers et al. (2020, in preparation). Our ALMA and NOEMA data, together with previous detections of the CO(2–1) (Shao et al. 2019) and the CO(10–9) emission line, enable us to probe the CO SLED
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Continuum, line intensity maps, and spectrum (from left to right) of CO(8–7), (9–8), H₂O(2₀,₂–l₁,₁), and OH⁺(l₁–0₁) (from top to bottom) observed by ALMA. The white cross represents the Gaia position of the quasar (Shao et al. 2019). The size of the white cross demonstrates the astrometric uncertainty of the quasar position. The filled white ellipse on the lower left shows the FWHM of the beam. First column: continuum maps. The white contours denote [−2, 2, 4, 8, 16, 32, 64, 74] × σ (1σ = 15 µJy beam⁻¹) at 136.6 GHz and [−2, 2, 4, 8, 16, 32, 64] × σ (1σ = 15 µJy beam⁻¹) at 141.1 GHz. Second column: spectra line intensity maps. The white contours denote [−2, 2, 4, 8, 16, 25] × σ (1σ = 0.042 mJy beam⁻¹ km s⁻¹), [−2, 2, 4, 8, 16, 22] × σ (1σ = 0.041 mJy beam⁻¹ km s⁻¹), [−2, 2, 4, 8, 15] × σ (1σ = 0.038 mJy beam⁻¹ km s⁻¹) and [−2, 2, 3, 4] × σ (1σ = 0.029 mJy beam⁻¹ km s⁻¹) for CO(8–7), (9–8), H₂O(2₀,₂–l₁,₁) and OH⁺(l₁–0₁) respectively. Third column: the yellow histogram represents the spectra extracted from the 2σ contour on the intensity map for spatially resolved CO(8–7), (9–8), H₂O(2₀,₂–l₁,₁) lines, and from the peak pixel for spatially unresolved OH⁺(l₁–0₁) emission. The red solid line is a single Gaussian profile fit to the spectral line. The Gaussian fit to the OH⁺(l₁–0₁) spectral line should be taken with caution since there are possibly both emission and absorption features, and the current sensitivity is not enough to confirm these features. The spatial offsets between the continuum and spectra line emissions are within the uncertainty of the quasar position.</figcaption></figure>

of J2310+1855 from *J* = 2 to *J* = 13, making it the most
complete CO SLED ever obtained for a *z* ≳ 6 quasar. We here
use the radiative transfer model MOLPOP-CEP to investigate

the physical conditions of the molecular gas, including the
kinetic temperature *T*<sub>kin</sub>, molecular hydrogen density *n*(*H*₂),
and CO column density *N*(*CO*).
<table><caption>**Table 2**

Spectral Line Observations</caption><thead><tr><th>Line</th><th>z<sub>line</sub></th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>S<sub>δν</sub><br>(Jy km s<sup>-1</sup>)</th><th>Beam Size<br>(arcsec)</th><th>Source Size<br>(arcsec)</th><th>Luminosity<br>(10<sup>9</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><th>(9)</th></tr></thead><tbody><tr><td>CO(2–1)</td><td>6.0029 ± 0.0005</td><td>484 ± 48</td><td>0.18 ± 0.02</td><td>0.61 × 0.59</td><td>(0.60 ± 0.18) × (0.40 ± 0.21)</td><td>0.021 ± 0.002</td><td>VLA</td><td>S19</td></tr><tr><td>CO(5–4)</td><td>6.0030 ± 0.0004</td><td>409 ± 44</td><td>0.89 ± 0.09</td><td>1.67 × 1.37</td><td>...</td><td>0.254 ± 0.026</td><td>NOEMA</td><td>L19</td></tr><tr><td>CO (6–5)</td><td>6.0025 ± 0.0007</td><td>456 ± 64</td><td>1.52 ± 0.13</td><td>5.4 × 3.9</td><td>...</td><td>0.520 ± 0.045</td><td>PdBI</td><td>W13</td></tr><tr><td>CO(6–5)</td><td>6.0028 ± 0.0003</td><td>361 ± 9</td><td>1.12 ± 0.06</td><td>0.6 × 0.4</td><td>(0.33 ± 0.06) × (0.20 ± 0.04)</td><td>0.383 ± 0.021</td><td>ALMA</td><td>F18</td></tr><tr><td>CO(6–5)</td><td>6.0030 ± 0.0003</td><td>422 ± 30</td><td>1.05 ± 0.07</td><td>1.42 × 1.19</td><td>(0.74 ± 0.34) × (0.46 ± 0.28)</td><td>0.359 ± 0.024</td><td>NOEMA</td><td>L19</td></tr><tr><td>CO(8–7)</td><td>6.0028 ± 0.0001</td><td>390 ± 15</td><td>1.53 ± 0.05</td><td>0.79 × 0.75</td><td>(0.46 ± 0.09) × (0.21 ± 0.10)</td><td>0.699 ± 0.023</td><td>ALMA</td><td>L19</td></tr><tr><td>CO(9–8)</td><td>6.0031 ± 0.0002</td><td>376 ± 18</td><td>1.31 ± 0.06</td><td>0.77 × 0.63</td><td>(0.41 ± 0.10) × (0.32 ± 0.11)</td><td>0.673 ± 0.030</td><td>ALMA</td><td>L19</td></tr><tr><td>CO(10–9)</td><td>...</td><td>...</td><td>1.04 ± 0.17</td><td>...</td><td>...</td><td>0.594 ± 0.097</td><td>...</td><td>Rpr</td></tr><tr><td>CO(12–11)</td><td>6.0030 ± 0.0008</td><td>451 ± 81</td><td>0.78 ± 0.13</td><td>2.08 × 1.62</td><td>...</td><td>0.534 ± 0.089</td><td>NOEMA</td><td>L19</td></tr><tr><td>CO(13–12)</td><td>6.0028 ± 0.0007</td><td>324 ± 75</td><td>0.49 ± 0.11</td><td>1.91 × 1.53</td><td>...</td><td>0.363 ± 0.082</td><td>NOEMA</td><td>L19</td></tr><tr><td>H<sub>2</sub>O</td><td>6.0028 ± 0.0003</td><td>398 ± 28</td><td>0.70 ± 0.05</td><td>0.72 × 0.68</td><td>(0.39 ± 0.14) × (0.25 ± 0.21)</td><td>0.343 ± 0.024</td><td>ALMA</td><td>L19</td></tr><tr><td>OH<sup>+</sup></td><td>...</td><td>320 ± 313</td><td>0.13 ± 0.10</td><td>0.77 × 0.63</td><td>...</td><td>0.067 ± 0.051</td><td>ALMA</td><td>L19</td></tr><tr><td>[CII]<sub>1588nm</sub></td><td>6.0031 ± 0.0002</td><td>393 ± 21</td><td>8.83 ± 0.44</td><td>0.72 × 0.51</td><td>(0.55 ± 0.05) × (0.40 ± 0.07)</td><td>8.310 ± 0.414</td><td>ALMA</td><td>W13</td></tr></tbody></table>

**Note.** Column 1: Line ID; Columns 2–4: Redshift, linewidth 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: Source size deconvolved from the beam in FWHM; Column 7: Line luminosity, and calibration uncertainties are not included in the error bars; Column 8: Facilities; Column 9: References: This paper (L19); Shao et al. (2019) (S18); Wang et al. (2013) (W13); Feruglio et al. (2018) (F18), Riechers et al. (2020, in preparation) (Rpr).

<table><caption>**Table 3**

Continuum Properties</caption><thead><tr><th>Frequency<br>(GHz)<br>(1)</th><th>S<sub>δν</sub><br>(mJy)<br>(2)</th><th>Rms<br>(μJy beam<sup>-1</sup>)<br>(3)</th><th>Beam Size<br>(arcsec)<br>(4)</th><th>Source Size<br>(arcsec)<br>(5)</th></tr></thead><tbody><tr><td>80.6</td><td>0.22 ± 0.04</td><td>16</td><td>1.68 × 1.37</td><td>(1.59 ± 0.52) × (0.21 ± 0.45)</td></tr><tr><td>96.0</td><td>0.29 ± 0.03</td><td>15</td><td>1.42 × 1.19</td><td>(0.83 ± 0.25) × (0.31 ± 0.39)</td></tr><tr><td>136.6</td><td>1.28 ± 0.03</td><td>15</td><td>0.75 × 0.72</td><td>(0.34 ± 0.04) × (0.22 ± 0.06)</td></tr><tr><td>141.1</td><td>1.42 ± 0.03</td><td>15</td><td>0.80 × 0.65</td><td>(0.27 ± 0.03) × (0.22 ± 0.06)</td></tr><tr><td>200.9</td><td>3.88 ± 0.04</td><td>45</td><td>2.06 × 1.62</td><td>...</td></tr><tr><td>215.9</td><td>4.46 ± 0.05</td><td>45</td><td>1.92 × 1.54</td><td>...</td></tr></tbody></table>

**Note.** Column 1: continuum frequency in observed frame; Columns 2–3: continuum flux density and rms; Columns 4–5: beam size and source size deconvolved from beam in FWHM.

## 4.1. Method

MOLPOP-CEP is a universal code that enables exact solutions of multi-level line emissions radiative transfer problems for all of the atoms/molecules that have atomic/molecular data in the Leiden Atomic and Molecular Database database (Elitzur & Asensio Ramos 2006; Asensio Ramos & Elitzur 2018). This code assumes a slab geometry with the emitting region divided into several zones and treats the radiative transfer problem with a coupled escape probability (CEP) method that aims to solve the coupled level population equations of different zones under consideration. Dividing the geometry into several zones makes it possible to solve the level populations as a function of depth into the line-emitting region and also leads to more accurate solutions compared to previous large velocity gradient models.

The physical parameters as inputs control the number of zones into which the geometry is divided and the physical parameters of the individual zone in the slab geometry. Even for a uniform physical parameter setup in the whole geometry, the division of geometry into zones is necessary to increase the accuracy of the results. This is because, for optically thick lines, the strength of radiative reactions changes with distances to the surface, and the transition level population distributions depend on positions in

the geometry (Asensio Ramos & Elitzur 2018). In the slab, each zone in principle can have different physical parameter setups including: (1) the zone width $\Delta L$, (2) the gas density within the zone $n(H_2)$, (3) kinetic temperature $T_{kin}$, (4) molecular abundance, and (5) local linewidth (which corresponds to the line absorption/emission profile in each point in the geometry). In addition to these, MOLPOP-CEP allows the inclusion of an external radiation field. Given these physical parameters, the code will then solve the coupled level population problem between zones and finally predict the emergent intensities of the emission lines that can be directly compared to the observations. We assume uniform parameters for each slab. The accuracy of MOLPOP-CEP solutions increases with the number of zones. Here, we divide the geometry into 10 zones for the model calculation, as is suggested by Asensio Ramos & Elitzur (2018).

We generate a grid of slab models through varying the physical parameters of greatest interest—the gas density $n(H_2)$, temperature $T_{kin}$, and zone width $\Delta L$. For the other two zone parameters, we fix the molecular abundance of CO to $X_{CO} = 10^{-4}$ (Milky Way, Blake et al. 1987), and local linewidth to 1 km s$^{-1}$ for all of the model calculations. In addition, we include the CMB at the quasar redshift of 19.12 K, because the hot CMB at high-redshift will (1) act as an extra heating source of the CO emission, and (2) serve as a continuum background (Da Cunha et al. 2013;
<figure><img src="image_3.png" /></figure>

**Figure 2.** Continuum, line intensity maps, and spectrum (from left to right) of CO(5–4), (6–5), (12–11), and (13–12) (from top to bottom) observed by NOEMA. The white cross on continuum and line intensity map represents the *Gaia* position of the quasar (Shao et al. 2019). The size of the white cross demonstrates the quasar location astrometric uncertainty. The filled white ellipse on the lower left shows the FWHM of the beam. For the continuum maps (first column), contours denote [−2, 2, 4, 8] × σ (1σ = 16 µJy beam⁻¹) at 80.6 GHz, [−2, 2, 4, 8, 14] × σ (1σ = 15 µJy beam⁻¹) at 96.0 GHz, [−2, 2, 4, 8, 16, 32, 64] × σ (1σ = 45 µJy beam⁻¹) at 200.9 GHz and [−2, 2, 4, 8, 16, 32, 64] × σ (1σ = 45 µJy beam⁻¹) at 215.9 GHz. For the spectra line intensity maps (second column), contours denote [−2, 2, 4, 8, 10] × σ (1σ = 0.08 mJy beam⁻¹·km s⁻¹), [−2, 2, 4, 8, 11] × σ (1σ = 0.08 mJy beam⁻¹·km s⁻¹), [−2, 2, 4, 7] × σ (1σ = 0.10 mJy beam⁻¹·km s⁻¹) and [−2, 2, 4, 5] × σ (1σ = 0.10 mJy beam⁻¹·km s⁻¹) for CO(5–4), (6–5), (12–11), and (13–12), respectively. Third column: yellow histogram represents spectrum extracted from the peak pixel (all spectra lines are spatially unresolved), and the red solid line is a single Gaussian profile fit to the spectral line. The peak positions of the continuum and spectra line emissions are within uncertainty of the quasar position.
<figure><img src="image_4.png" /><figcaption>Figure 3. Velocity (first row) and velocity dispersion (second row) of [C II]158 µm, CO(8–7), (9–8), and H₂O(2₀,₂–1₁,₁) (left to right). The black ellipse on the left shows the FWHM of the beam. We generate the velocity and velocity dispersion maps with the pixels of ≥3.0σ values in the line-emitting channels, and the zero velocity corresponds to the [C II]158 µm redshift of 6.0031 (Wang et al. 2013). The white contours on the velocity maps are [−4, −2, 0, 2, 4, 6] × 20 km s⁻¹ for [C II]158 µm, [−6, −4, −2, 0, 2, 4] × 20 km s⁻¹ for CO(8–7), [−6, −4, −2, 0, 2, 4] × 20 km s⁻¹ for CO(9–8), and [−6, −4, −2, 0, 2, 4] × 20 km s⁻¹ for H₂O(2₀,₂–1₁,₁). We find in the velocity maps that CO(8–7), (9–8), and H₂O(2₀,₂–1₁,₁) generally follow the velocity gradient observed from [C II]158 µm (Wang et al. 2013) from northeast to southwest.</figcaption></figure>

Zhang et al. 2016). In each model, all 10 zones have the same physical parameter setup; thus, each model in the grid can be described by a uniform $n(H_2)$, $T_{kin}$ and $\Delta L$. The CO column density within the geometry, which is the sum of the column densities of 10 zones, is proportional to the zone width $\Delta L$ through:

$$N(\text{CO}) = 10 \times n(H_2) \times \Delta L \times X_{\text{CO}}. \quad (1)$$

In the remainder of the paper, we use $n(H_2)$, $T_{kin}$, and $N(\text{CO})$ to characterize the physical condition of each slab grid. As it is easier to use $N(\text{CO})$ rather than $\Delta L$ to make comparisons with constraints from observations (see Section 4.2). The final grid covers the typical physical conditions of the molecular clouds with a temperature range of 20–800 K, density of $10^3–10^8$ cm⁻³, and CO column density of $10^{14}–10^{21}$ cm⁻² (this corresponds to a H₂ column density of $10^{18}–10^{25}$ cm⁻² for our assumed $X_{CO}$ of $10^{-4}$). More details about the grid are listed in Table 4.

We use the model grid to fit our observed CO SLED of J2310+1855. The fitting procedure is as follows: we first apply the least-squares method to find the best-fitting results. We also use the Bayesian code emcee (Foreman-Mackey et al. 2013) to efficiently explore the parameter space and get the posterior probability distributions of all parameters considered. Emcee is an extensible, pure-Python implementation that is designed for Bayesian parameter estimation using Ensemble samplers with affine invariance (Goodman & Weare 2010).

## 4.2. Parameter Constraints

During the fitting procedure, we set constraints for the three parameters: $T_{kin}$, $n(H_2)$, and $N(\text{CO})$. The fact that $T_{kin}$ is hotter than the background CMB radiation at redshift 6 is a prior, which sets $T_{kin} > 19.12$ K. As for column density, one prior is that the total amount of gas producing the observed CO luminosities should be no more than the total dynamical mass of this system, this leads to:

$$N(\text{CO})\Phi_A < \frac{M_{\text{dyn}}X_{\text{CO}}}{\mu m_{\text{H}_2}A} [\text{cm}^{-2}], \quad (2)$$

where $\Phi_A$ is the filling factor, $M_{dyn}$ is the dynamic mass, $A$ is the source area in cm², $\Delta V$ is the CO linewidth in km s⁻¹, $\mu$ is the mean molecular weight, $m_{H_2}$ is the H₂ molecule mass, and $X_{CO}$ is the CO abundance. The source size and gas dynamical mass are adopted from Shao et al. (2019), with $M_{dyn} \approx 4.3 \times 10^{10} M_\odot$ and $A \approx \frac{\pi}{4} \times 0.60 \times 0.40 \text{ arcsec}^2$. Adopting a CO-to-H₂ abundance ratio $X_{CO} = 10^{-4}$, and $\mu = 1.4$ for mean molecular weight, the final constraint is:

$$N(\text{CO})\Phi_A < 3.2 \times 10^{19} [\text{cm}^{-2}]. \quad (3)$$

The other prior for column density is that it should be less than the source gas volume density integrated along the line of sight for line-emitting regions. This leads to:

$$N(\text{CO}) < n(H_2) \times X_{\text{CO}} \times S [\text{cm}^{-2}], \quad (4)$$
<table><caption>**Table 4**
MOLPOP-CEP Grid Parameter Ranges</caption><thead><tr><th>Input Parameters<br>(1)</th><th>Range<br>(2)</th><th>Grid Step<br>(3)</th><th>Grid Number<br>(4)</th><th>Unit<br>(5)</th></tr></thead><tbody><tr><td>Kinetic Temperature (T<sub>kin</sub>)</td><td>20–800</td><td>ΔT<sub>kin</sub> = 13</td><td>61</td><td>K</td></tr><tr><td>Volume Density (n(H<sub>2</sub>))</td><td>10<sup>3</sup>–10<sup>8</sup></td><td>ΔLog(n<sub>H<sub>2</sub></sub>/cm<sup>-3</sup>) = 0.25</td><td>21</td><td>cm<sup>-3</sup></td></tr><tr><td>Column Density (N(CO))</td><td>10<sup>14</sup>–10<sup>21</sup></td><td>ΔLog(N<sub>CO</sub>/cm<sup>-2</sup>) = 0.5</td><td>15</td><td>cm<sup>-2</sup></td></tr></tbody></table>

Note. Column 1: input parameters to generate the grid; Column 2: parameter ranges; Column 3: steps of parameters in log space; Column 4: the resulting number of grid for a specific parameter; Column 5: units

where *S* is the source size along the line of sight, and *n*(*H*₂) is the *H*₂ volume density. Assuming a largest diameter of CO-emitting region of 0.60 (Shao et al. 2019), we set a constraint on both *N*(*CO*) and *n*(*H*₂) as follows:

$$ \frac{N(\text{CO})}{n(\text{H}_2)} < 1.1 \times 10^{18} [\text{cm}]. \quad (5) $$

### 4.3. Fitting Results

We endeavor to probe the physical conditions of the molecular gas in J2310+1855 by fitting a one-component model to the CO SLED. Calibration uncertainties are included in all of the modeling processes throughout the paper. We first fit the grid models to the observed CO SLED with the least-squares method. The left panel of Figure 4 shows the best-fitting result, suggesting a “warm” and “dense” gas component with kinetic temperature of $T_{\text{kin}} = 228$ K, density of $\log(n(\text{H}_2)/\text{cm}^{-3}) = 4.75$, and column density of $\log(N(\text{CO})/\text{cm}^{-2}) = 17.5$. We then search for all possible physical conditions that fit the observational CO SLED with the emcee code. The posterior probability distributions of the three parameters are shown in Figure 5. This indicates that the data can be fitted with a “warm” and “dense” gas component with parameter range of $T_{\text{kin}} \approx 167_{-0.56}^{+1.53}$ K, $\log(n(\text{H}_2)/\text{cm}^{-3}) \approx 5.11_{-0.58}^{+1.83}$, and $\log(N(\text{CO})/\text{cm}^{-2}) \approx 17.28_{-0.42}^{+0.33}$.18 But the current best-fit model fails to reproduce the very high CO(8–7) line flux detected.

Previous CO SLED modeling from local to high-redshift galaxies/AGNs suggest different gas physical properties in different systems. Exempli gratia in the *z* = 2.56 quasar cloverleaf (Riechers et al. 2011) and the *z* = 6.34 starburst galaxy HLFS3 (Riechers et al. 2013), a single gas component is able to reproduce the observed CO SLED. In addition, more than one gas component is found in the CO SLED analysis of local starburst galaxies, (U)LIRGs, and even the quasars at the highest redshift (e.g., M82, NGC 1068, Mrk 231, APM 08279 +5255, J0100+2802, J0439+1634, and J1148+5251 (Weiß et al. 2005, 2007; Panuzzo et al. 2010; van der Werf et al. 2010; Spinoglio et al. 2012; Gallerani et al. 2014; Wang et al. 2019a; Yang et al. 2019)). The CO SLED of the *z* = 6.3 quasar J0100 +2802 suggests two components of gas, with a “cold” component with $T_{\text{kin}} \approx 24$ K, $\log(n(\text{H}_2)/\text{cm}^{-3}) \approx 4.5$ and a “warm” component with $T_{\text{kin}} \approx 224$ K, $\log(n(\text{H}_2)/\text{cm}^{-3}) \approx 3.6$ (Wang et al. 2019a). The CO SLED of J0439+1634 at *z* = 6.5 indicates a “cold” component with $T_{\text{kin}} \approx 23$ K and $\log(n(\text{H}_2)/\text{cm}^{-3}) \approx 4.1$ in combination with a “warm” component with $T_{\text{kin}} \approx 140$ K and $\log(n(\text{H}_2)/\text{cm}^{-3}) \approx 4.5$ (Yang et al. 2019). The cold component was thought to be associated

with the sub-mm/mm-detected dust component powered by active star formation with temperatures of 40~60 K (assuming optically thin; Beelen et al. 2006; Leipski et al. 2013). The dust continuum SED fitting of J2310+1855 indicates a dust temperature of ~40 K in the optically thin dust assumption of Shao et al. (2019) or 76 K in the optically thick assumption of Carniani et al. (2019). Both are much lower than the one-component fitting result of $T_{\text{kin}} = 228$ K. As the one-component best-fit model fails to explain the observed CO (8–7) flux (left panel of Figure 4), it is possible that there is an additional “cold” gas component physically associated with the sub-mm/mm-detected dust in J2310+1855. Motivated by the above, we examine whether the data can be explained with a two-component model.

The two-component model fitting to the data requires eight parameters, including the physical parameters ($T_{\text{kin}}$, $n(\text{H}_2)$, and $N(\text{CO})$) and the normalization of each component. Our data are insufficient to constrain all eight parameters. Considering that the “cold” component is usually physically connected to and has a similar temperature to that of the cold dust (as is explained in detail in the previous paragraph), and the dust temperature of J2310 +1855 is not well constrained ($T_{\text{dust}}$ ranges between 40 and 80 K depending on the dust model assumed), we will fix the “cold” component to the typical “cold” gas physical conditions observed in *z* ~ 6 quasar in the following analysis. With this assumption, we are fitting only five instead of eight parameters. As a consequence, the model parameters could be better constrained. We adopt a typical “cold” gas component with physical parameters of $T_{\text{kin}} \approx 50$ K and $\log(n(\text{H}_2)/\text{cm}^{-3}) = 4.2$ that is observed in a “typical” *z* ~ 6 quasar J1148+5251 (Riechers et al. 2009). Because column density is not one of the model parameters in Riechers et al. (2009), we use a column density of $\log(N(\text{CO})/\text{cm}^{-2}) = 18.0$ for the “cold” model. We note that the resulting $X_{\text{CO}}/\text{dv}/\text{dr}$ here is different from that in Riechers et al. (2009), because in the MOLPOP-CEP model, the $\text{dv}/\text{dr}$ is derived from the first principle and is different in different places in the whole geometry. The final “cold” model using the set of parameters we adopted here can well represent the observational CO SLED of J1148+5251 presented in Riechers et al. (2009). The right panel of Figure 4 shows the minimum $\chi^2$ fitting result of the two-component model. The resulting “warm” component with the minimum $\chi^2$ has a physical condition of $T_{\text{kin}} = 306$ K, $\log(n(\text{H}_2)/\text{cm}^{-3}) = 5.25$, and $\log(N(\text{CO})/\text{cm}^{-2}) = 15.5$. We find in the fitting result that the “cold” component (the J1148 +5251 model) dominates the low *J* (*J* = 2) part and contributes to 77% of the observed CO(2–1) flux. As we have already mentioned before, the low *J* CO emission lines trace the total molecular gas mass; thus, the “cold” component dominates the total molecular gas mass. In the mid-*J* (*J* = 5, 6) part, the contribution of the “cold” component decreases and only accounts for ~30% of the observed CO fluxes. And in the high *J* (*J* ≥ 8) part, the “cold”

<footnote n="18">18 The resulting parameter ranges are consistent within $1\sigma$ between including and excluding the CO(10–9) line in the fitting procedure.</footnote>
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** CO SLED fitted with one- (left panel) or two-component (right panel) models. The black squares are the CO fluxes of J2310+1855 with the calibration uncertainties included. The CO(5–4), (6–5), (8–7), (9–8), (12–11), and (13–12) data are from this work. The CO(2–1) data is taken from Shao et al. (2019), the CO(10–9) data will be presented in Riechers et al. (2020, in preparation), and the upper limit of the CO(17–16) line is from Carniani et al. (2019). Details are presented in Table 2. Left panel: the black solid line represents the best one-component fitting result with the minimum $\chi^2$ ($T_{\text{kin}}$ = 228 K, $\log(n(\text{H}_2)/\text{cm}^{-3})$ = 4.75, and $\log(N(\text{CO})/\text{cm}^{-2})$ = 17.5). Right panel: least-squares fitting result with two-component model. The blue solid line represents the J1148+5251 model with $T_{\text{kin}}$ = 50 K, $\log(n(\text{H}_2)/\text{cm}^{-3})$ = 4.20 and $\log(N(\text{CO})/\text{cm}^{-2})$ = 18.0 (Riechers et al. 2009). The red solid line is the “warm” component with $T_{\text{kin}}$ = 306 K, $\log(n(\text{H}_2)/\text{cm}^{-3})$ = 5.25 and $\log(N(\text{CO})/\text{cm}^{-2})$ = 15.5.</figcaption></figure>

component contribution is negligible. The “warm” component, that barely contributes to the total molecular gas mass, dominates the CO SLED from the mid- to the high-$J$ ($J \ge 5$) part of the overall CO SLED. The posterior probability distributions of the parameters calculated by the emcee code are shown in Figure 5. It suggest a “warm” and “dense” component with parameter range of $T_{\text{kin}} \approx 306_{-149}^{+263}$ K, $\log(n(\text{H}_2)/\text{cm}^{-3}) \approx 5.22_{-0.49}^{+1.04}$, and $\log(N(\text{CO})/\text{cm}^{-2}) \approx 15.29_{-1.17}^{+1.34}$$^{19}$.

To summarize, the best one-component “warm” and “dense” model reproduces the observed CO SLED, in general, except for an underestimation of the CO(8–7) flux. The two-component fitting result suggests that the CO SLED at $J \ge 5$ is dominated by a “warm” and “dense” gas component, while the “cold” component barely contributes to the mid- to the high-$J$ CO fluxes but dominates the total molecular gas mass. Either the one- or two-component model suggests that the CO SLED detected within the nuclear region (source size of ~2 kpc) of the quasar host is dominated by a “warm” and “dense” gas component at $J \ge 5$.

# 5. Discussion

## 5.1. CO Emission in J2310+1855 Compared with Local Starburst Systems

The CO SLED reveals the physical conditions of molecular gas (e.g., the illuminating radiation field strength, kinetic temperature, volume density, and column density). We first compare the CO SLED of J2310+1855 with local starburst systems. Figure 6(a) shows the CO SLED of J2310+1855 compared with two local starburst samples. They are local (U) LIRGs sample consisting of 29 (U)LIRGs (Rosenberg et al. 2015) and local normal + starburst galaxy sample consisting of 43 star-forming galaxies (SFGs) and 124 (U)LIRGs (Liu et al. 2015). Although there are actually a small number of AGNs

contained in the two (U)LIRGs comparison samples, they are confirmed to impact negligibly on both the CO flux and the infrared luminosity (Rosenberg et al. 2015). For the Rosenberg et al. (2015) sample, we exclude NGC 6240 in a sample mean calculation (because this AGN represents a very extreme CO SLED, see details in Section 5.2). Through comparisons, we find that the peak of the CO SLED is $\lesssim 4$ for all local starburst samples, while the J2310+1855 CO SLED peaks at much higher $J$ transitions at $J = 8$. The CO emission lines of J2310+1855 show higher excitation compared to all local (U)LIRG samples (i.e., the CO flux is higher than the average of all local ULIRG samples and is also well above the range of all comparison samples, especially for $J \ge 8$). In addition, we compare J2310+1855 with a representative example of local starburst galaxy M82 in Figure 6(b). We get a similar result to that of the starburst samples. The CO SLED of M82 peaks at $J = 5$ and decreases dramatically at $J \ge 8$, contrary to J2310+1855, which peaks at $J = 8$ and is luminous even at $J \ge 10$.

The high CO excitation detected in the J2310+1855 nuclear region (source size of ~2 kpc in FWHM) may indicate other heating mechanisms besides the UV heating from massive young stars (e.g., mechanical heating by shocks, X-ray heating from AGNs), or very intense UV radiation field (e.g., large UV photon flux produced by a result of both the quasar and the star formation, see Section 5.5 for further discussions about these heating mechanisms).

## 5.2. CO Emission in J2310+1855 Compared with Local AGNs

We compare the CO emission lines of J2310+1855 with some representative local AGNs (Mrk 231, NGC 1068, and NGC 6240). The AGNs we selected are thoroughly studied local AGNs that represent different CO heating mechanisms. NGC 1068 is one of the closest AGNs, whose high $J$ CO emission lines in the circumnuclear disk (CND) are best explained by an XDR model, and the starburst ring at larger radii that dominates the molecular gas mass is best fitted with a

<footnote n="19">¹⁹ The resulting parameter ranges are consistent within $1\sigma$ between including and excluding the CO(10–9) line in the fitting procedure.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Posterior probability distributions of the three parameters $T_{\text{kin}}$ (K), $\log(n(\text{H}_2)/\text{cm}^{-3})$, and $\log(N(\text{CO})/\text{cm}^{-2})$. The plotted contours show the 95% and 68% confidence intervals. Left panel: posterior probability distribution of the parameters for the one-component model. The resulting MCMC result is $T_{\text{kin}} \approx 167_{-56}^{+153}$ K, $\log(n(\text{H}_2)/\text{cm}^{-3}) \approx 5.11_{-0.58}^{+1.83}$, and $\log(N(\text{CO})/\text{cm}^{-2}) \approx 17.28_{-0.42}^{+0.33}$. We note that the median and uncertainties here are calculated based on the 16th, 50th, and 84th percentiles of the samples in the marginalized distributions. Right panel: posterior probability distribution of the “warm” component in the two-component model fit to the data. The fitting result suggests the “warm” component with $T_{\text{kin}} \approx 306_{-149}^{+263}$ K, $\log(n(\text{H}_2)/\text{cm}^{-3}) \approx 5.22_{-0.49}^{+1.04}$, and $\log(N(\text{CO})/\text{cm}^{-2}) \approx 15.29_{-1.17}^{+1.34}$.</figcaption></figure>

PDR model (Spinoglio et al. 2012). Mrk 231 requires an XDR model to fit the high *J* CO emission lines in the central 160 pc molecular region, although the low *J* CO lines are mainly from the PDR component at larger distances (van der Werf et al. 2010). NGC 6240 is a local AGN with three nuclei, the CO emission lines are unlikely to correlate with the position of either AGN nuclei, and mechanical heating is proposed to interpret the extremely excited CO SLED and also the optical ISM emission lines (Meijerink et al. 2013; Kollatschny et al. 2019).

Figure 6(b) shows the CO SLED of J2310+1855 and the local AGNs. The CO SLED of NGC 1068 peaks at $J \le 4$, and then decreases rapidly with increasing *J*. Mrk 231 shows an extreme CO SLED that peaks at *J* = 5 and displays a high normalized CO flux even at *J* ≈ 10–13. The differences between NGC 1068 and Mrk 231 can be explained by different X-ray energy input to the CO heating, evident from the higher X-ray flux derived in Mrk 231 (28 erg cm⁻² s⁻¹, van der Werf et al. 2010) than NGC 1068 (9 erg cm⁻² s⁻¹, Spinoglio et al. 2012). NGC 6240 shows the most extreme CO SLED among these three AGNs. The quasar J2310+1855 peaks at higher *J* (*J* = 8) compared to these local AGNs, i.e., the gas in the nuclear region of J2310+1855 has higher excitation compared to that of the local AGNs.

### 5.3. CO Emission in J2310+1855 Compared with High-redshift Systems

We also compare the CO SLED of J2310+1855 with high-redshift systems, including high-redshift (lensed) SMGs and quasars. Figure 6(a) shows the CO SLED of J2310+1855 and high-redshift SMGs: a sample of $z \sim 1.2-4.1$ SMGs (Bothwell et al. 2013) and a sample of $15 \le z \sim 2-4$ lensed

SMGs from Carilli & Walter (2013) and Yang et al. (2017). The high-redshift SMG CO SLEDs peak at $J \lesssim 6$, while J2310+1855 peaks at higher *J* (*J* = 8) than the two SMG samples. This is similar to the results when comparing the J2310+1855 CO SLED to local starburst samples.

We also select some well-known high-redshift quasars for comparison. Including two lensed quasars, the Cloverleaf at $z = 2.56$ and APM 08279+5255 at $z = 3.91$. The fit to the CO SLED detected in the very central region of APM 08279+5255 requires an XDR component dominating the high *J* CO emission lines (Bradford et al. 2011). We also include three $z \gtrsim 6$ quasars that are detected in at least four CO transitions: J1148+5251, J0439+1634, and J0100+2802 (Bertoldi et al. 2003; Walter et al. 2003; Beelen et al. 2006; Riechers et al. 2009; Gallerani et al. 2014; Wang et al. 2019a; Yang et al. 2019). Together with J2310+1855, this allows us to perform a systematic study of the CO emission lines in the quasar-starburst systems at the highest redshift. Figure 6(c) shows the CO SLED of J2310+1855 compared to other high-redshift quasars. APM 08279+5255 exhibits the most extreme CO SLED that is detected in the nuclear 550 pc size scale, and it represents a highly excited nuclear CO SLED exposed to the intense X-rays from the quasar. The limited CO detections in J0439+1634 and Cloverleaf suggests that maybe the CO SLED peaks at *J* = 8 ~ 10, similarly to J2310+1855. As for J1148+5251 and J0100+2802, we are not able to determine the CO SLED peak. This is because for J0100+2802, there is no CO observation at *J* = 8, 9 and *J* > 11, while the CO SLED of J1148+5251 is observed at *J* ≤ 7 and *J* = 17. The comparison of CO SLEDs between J2310+1855 and other high-*z* quasars suggest that a highly excited molecular gas component is common in the nuclear region of the quasar hosts. However, as
<figure><img src="image_8.png" /><figcaption>**Figure 6.** CO SLED normalized to CO(6–5) in normalized unit of Jy km s⁻¹. Left column: J2310+1855 (black squares) in comparison with the mean of four galaxy samples: the local (U)LIRGs (Rosenberg et al. 2015) (red); local normal + starburst galaxies (Liu et al. 2015) (blue); *z* ~ 1.2–4.1 SMGs (Bothwell et al. 2013) (magenta); and strongly lensed SMGs at *z* ~ 2–4 (Yang et al. 2017) (green). Middle column: J2310+1855 (black squares) in comparison with local starburst systems and AGNs. M82 (Weiß et al. 2005; Panuzzo et al. 2010) (green) is a representative example of the local starburst galaxy. The local representative AGNs are NGC 1068 (Spinoglio et al. 2012) (magenta), Mrk 231 (van der Werf et al. 2010) (blue), and NGC 6240 (Rosenberg et al. 2015)(red). Right column: J2310+1855 (black squares) in comparison with high-redshift quasars. The plotted quasars are APM 08279+5255 (Weiß et al. 2007; Riechers et al. 2009; Bradford et al. 2011) (gray), Cloverleaf (Bradford et al. 2009; Uzgil et al. 2016) (red), J1148+5251 (Bertoldi et al. 2003; Walter et al. 2003; Beelen et al. 2006; Riechers et al. 2009; Gallerani et al. 2014) (blue), J0439+1634 (Yang et al. 2019) (green), and J0100+2802 (Wang et al. 2019a) (magenta).</figcaption></figure>

described above, the shape of the CO SLED of these systems are different from object to object. Due to the lensing effect, the CO SLED of APM 08279+5255 may represent molecular gas on <1 kpc scale. The CO SLED of J2310+1855 is not as extreme as APM 08279+5255, and more comparable to the Cloverleaf. When compared to the two *z* > 6 quasars, J0100+2802 and J0439+1634, that have available CO data at *J* = 8 ~ 11, J2310+1855 is more single-peaked no flatten or turn over around *J* = 6. It is possible that the cold star-forming component contribute more to the flux of the mid-*J* CO in the cases of J0100+2802 and J0439+1634.

## 5.4. H₂O and OH⁺ Emission

Recent studies found linear relations between the infrared luminosity and the H₂O(2₀,₂−I₁,₁) luminosity in local and high-redshift infrared bright systems, suggesting that the excitation of this water transition is dominated by infrared pumping (van der Werf et al. 2011; Yang et al. 2016; Jarugula et al. 2019). We first compare the H₂O detection of J2310+1855 to the local and high-redshift (U)LIRGs. Figure 7 shows the L<sub>H₂O</sub> and L<sub>IR</sub> for local and high-redshift (U/Hy)LIRGs and AGNs. The black dashed line is the best fit to the local and high-redshift (U/Hy)LIRGs presented in Figure 3 of Yang et al. (2016). We consider the infrared luminosity of J2310+1855 in two cases: (1) we adopt total infrared luminosity from both the quasar and the host galaxy (L<sub>IR(total)</sub>); (2) we use the infrared luminosity only from the host galaxy (L<sub>IR(galaxy)</sub>, Shao et al. 2019). The linear relation is fitted with (U/Hy)LIRGs (local and high *z*) and dusty star-forming galaxies, while all four plotted AGNs including J2310+1855 are well below this relation. This is because the AGNs not only provide the power source of water emission but also contribute significantly to the infrared luminosity (González-Alfonso et al. 2010; van der Werf et al. 2011). As for the host galaxy infrared luminosity case, J2310+1855 reveals a slightly higher water luminosity given its galaxy IR luminosity compared to the linear relation. The velocity dispersion map of the H₂O shows a velocity dispersion of <100 km s⁻¹, suggesting that it may not be the large velocity dispersion that contributes to the luminous water emission. Higher spatial resolution observations are required to confirm this. We find in J2310+1855 a slightly higher L<sub>H₂O</sub>/L<sub>IR</sub>

ratio than that for local and high-redshift (U/Hy)LIRGs. At *z* ~ 6, only a few quasars are detected in water emission (Bañados et al. 2015; Yang et al. 2019). A similar result is found in the *z* = 6.52 quasar J0439+1634, where a higher L<sub>H₂O</sub>(3₂,₁−3₁,₂)/L<sub>IR</sub> ratio is found as compared to the linear relation (Yang et al. 2019).

To further investigate the heating sources of molecular gas, we also study the ratio between H₂O and CO in J2310+1855. Extremely luminous H₂O emission is not expected in typical PDRs. Exempli gratia, in the Orion bar (a representative dense PDR with *n* ~ 10⁵ cm⁻³ illuminated by an intense FUV radiation field of G₀ = 4 × 10⁴), the L<sub>H₂O</sub>/L<sub>CO</sub>(6–5) ratio is 0.20 (Habart et al. 2010; Putaud et al. 2019). Another example is the local starburst galaxy M82, which shows a ratio of 0.06 (Kamenetzky et al. 2012). If the physical/chemical condition is dominated by shocks or X-rays, then the L<sub>H₂O</sub>/L<sub>CO</sub>(6–5) ratio can even be as high as unity. NGC 1266 is an S0 galaxy highly excited in molecular gas, and Pellegrini et al. (2013) found an L<sub>H₂O</sub>/L<sub>CO</sub>(6–5) ratio of 0.96 that can be only explained by shocks. Mrk 231 is representative of the molecular gas heated by X-rays, and it has an L<sub>H₂O</sub>/L<sub>CO</sub>(6–5) ratio of 1.10 ± 0.17 (González-Alfonso et al. 2010). The quasar J2310+1855 exhibits an L<sub>H₂O</sub>/L<sub>CO</sub>(6–5) ratio of 0.97 ± 0.09, which is comparable to Mrk 231 and NGC 1266 but is much higher than Orion bar and M82. These suggest that the molecular gas heating is not likely dominated by PDR. In addition, we also consider the high-redshift lensed SMGs for comparison (Omont et al. 2013; Yang et al. 2013, 2016, 2017). These lensed SMGs are all starbursts that are exposed to intense radiation field and shocks, and some have extremely luminous water emission lines that cannot be purely explained by PDR. For the majority of the SMGs, the L<sub>H₂O</sub>/L<sub>CO</sub>(6–5) ratio is found to be less than 0.8, while only quite a few objects have ratios of ~1, which is unlikely PDR. The even higher L<sub>H₂O</sub>/L<sub>CO</sub>(6–5) ratio of J2310+1855 compared to typical high-redshift lensed SMGs might suggest additional gas heating by the central luminous quasar. We also compare the H₂O/CO ratio of the quasar with other *z* ~ 6 quasars that are detected in water emission. J0439+1634 was detected in the H₂O(3₂,₁−3₁,₂) emission and suggests a L<sub>H₂O</sub>(3₂,₁−3₁,₂)/L<sub>CO</sub>(6–5) ratio of 1.23 ± 0.22. Adopting a mean H₂O(3₂,₁−3₁,₂)/H₂O(2₀,₂−I₁,₁) of 1.4 estimated from SMGs in
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Left panel: $L_{H_2O(2\text{0}_2\text{-}l_{1,1})}$ vs. $L_{IR}$ (from Figure 3 in Yang et al. 2016). The green squares are the local ULIRGs from Yang et al. (2013), the blue diamonds are the high-redshift U/HyLIRGs from Yang et al. (2016), Omont et al. (2013), and van der Werf et al. (2011), the magenta diamonds are the high-redshift (dusty) star-forming galaxies from Jarugula et al. (2019) and Apostolovski et al. (2019), and the red down-triangles represent local and high-redshift AGNs, namely Mrk 231 (González-Alfonso et al. 2010), SDP81 (Yang et al. 2016), and APM 08279+5255 (Bradford et al. 2011) from left to right. Note that all of the luminosities plotted are intrinsic luminosities that have been corrected for lensing. The gray and black squares mark two cases of J2310+1855: the former (gray square) is the water line to total infrared luminosity ratio with $L_{IR(total)} = (5.7 \pm 0.6) \times 10^{13} L_\odot$ (with contributions from both quasar and host galaxy), and the latter (black square) shows the water line to galaxy infrared luminosity ratio with $L_{IR(galaxy)} = (1.4 \pm 0.3) \times 10^{13} L_\odot$ (which is purely from the host galaxy). The black dashed line represents the best fit to the local and high z U/HyLIRGs (green squares and blue diamonds) with $L_{H_2O(2\text{0}_2\text{-}l_{1,1})} \sim L_{IR}^{1.06}$, see (Yang et al. 2016). Note that the definition of infrared luminosity is 8–1000 µm. Right panel: spectrum of CO(8–7) (magenta), (9–8) (green), and $H_2O(2\text{0}_2\text{-}l_{1,1})$ (blue) normalized to the peak flux densities.</figcaption></figure>

Yang et al. (2016), we estimate the $L_{H_2O}/L_{CO}$ (6–5) ratio in J0439+1634 of $0.86 \pm 0.17$, which is comparable to the value of J2310+1855.

We tentatively detected the $OH^+$ line at an S/N of 4. The spectra of the $OH^+$ line shows a P Cygni-like profile, with a possible absorption in the line center frequency. This might hint at possible outflows/inflows in this source, although further high-sensitivity observations are needed to confirm this. The ratios between $OH^+$, $H_2O^+$, and $H_3O^+$ reactive molecular ions are ideal tracers of the ionization rate and serve as the ionization source diagnostics. In J2310+1855, we are not able to constrain the ionization rate with the $OH^+$ detection alone. On the other hand, $OH^+$ and $H_2O$ are all oxygen-hydrogen species, and $OH^+$ can be formed by photodissociation of $H_2O$. Accordingly, we are expecting a higher $L_{OH^+}/L_{H_2O}$ ratio in the presence of cosmic-rays or X-rays. Mrk 231 is one of the best-studied AGNs that has been detected in a series of $OH^+$, $H_2O^+$, and $H_3O^+$ emission lines. González-Alfonso et al. (2018) detected very bright $OH^+$ emission with $L_{OH^+}/L_{H_2O}$ ratio of $0.37 \pm 0.13$. The ionization rate derived by making use of all of the molecular ions detected is very high, and cannot even be explained by its observed X-ray flux and requires ionization by cosmic-rays (i.e., the ionization rate produced by X-ray photons is ~1/10 the value required). The quasar J2310+1855 exhibits an $L_{OH^+}/L_{H_2O}$ ratio of $0.20 \pm 0.15$, comparable with that of Mrk 231. The observed X-ray luminosity in J2310+1855 is ~170× that observed in Mrk 231 (Vito et al. 2019). If we simply assume a similar ionization rate in J2310+1855 as that of Mrk 231 (as the $L_{OH^+}/L_{H_2O}$ ratios are comparable in these two AGNs), the X-ray photons from J2310+1855 are more than adequate to explain the observed $OH^+$ emission.

## 5.5. Molecular Gas Heating Mechanisms

In the CO SLED modeling of J2310+1855, we find the mid to high $J$ ($J \ge 5$) CO emission lines are dominated by a

“warm” and “dense” gas component with $T \gtrsim 150$ K, $\log(n(H_2)/\text{cm}^{-3}) \gtrsim 5$, and $\log(N(\text{CO})/\text{cm}^{-2}) \gtrsim 15.0$. Such “warm” and “dense” gas components are warmer and denser than typical values found in local (U)LIRGs (Papadopoulos et al. 2012, 2013), and comparable to the extreme cases, e.g., the NGC 1068 circumnuclear disk (Krips et al. 2011; Viti et al. 2014) and highly excited (U)LIRGs (Papadopoulos et al. 2012, 2013). Studies of local and high-redshift systems suggest that the molecular gas can be heated through (1) the UV heating from young massive stars or AGNs, (2) mechanical heating by shocks generated from supernovae or AGN outflows, (3) cosmic-ray heating from supernovae or AGNs, or (4) X-ray heating from the AGNs (Bradford et al. 2003, 2011; Spinoglio et al. 2012; Meijerink et al. 2013; Rosenberg et al. 2015; Uzgil et al. 2016). We inspect the most probable mechanism that contributes to the high CO excitation and luminous $H_2O$ and $OH^+$ emission lines observed in J2310+1855.

To investigate the PDR origin, we fit the CO SLED of J2310+1855 to the PDR grid (Meijerink et al. 2007), and the best-fitting result indicates an FUV flux of $G_0 = 1.0 \times 10^4$ and $n \sim 5.6 \times 10^5 \text{ cm}^{-3}$. This suggests a higher FUV flux than that reported in Carniani et al. (2019) due to the lack of information from the $J \ge 10$ CO lines in their study. Although the CO emission can be explained by a dense PDR illuminated by an intense FUV radiation field, the bright $H_2O$ and $OH^+$ emission lines (high $H_2O/\text{CO}$ and $OH^+/H_2O$ ratio) hint at different physical and chemical conditions than in a PDR, suggesting rather the presence of X-rays, shocks, or cosmic-rays that heat and ionize the molecular gas. In the nuclear region, both the star formation activity and the luminous quasar are capable of influencing the physical and chemical conditions of the molecular gas. The powerful quasar is able to participate in all of the possible gas heating mechanisms through the
X-rays it radiates, the shocks generated by AGN outflows, and
the cosmic-rays it produces.

6. Summary

We report new detections of CO(5–4), (6–5), (12–11) and
(13–12) with NOEMA and CO(8–7), (9–8), H₂O and OH⁺
with ALMA in the *z* = 6.003 quasar J2310+1855. This is the
most complete CO SLED ever obtained for a *z* ≥ 6 quasar. We
spatially resolved the CO(8–7), (9–8) and H₂O lines at similar
source sizes of ~2 kpc in FWHM, which are slightly more
compact than the [C II] 158 µm and CO(2–1) emission lines.
These suggest that the high *J* CO lines and the H₂O line are
probing the nuclear dense molecular regions closer to the
quasar. We analyze the physical conditions of the molecular
gas through CO SLED modeling, and we also compare the CO
emission lines from the quasar with local and high-redshift
starburst samples and some representative local and high-
redshift AGNs. The main results are summarized below.

1. The CO SLED of J2310+1855 at *J* ≥ 5 is dominated by a “warm” and “dense” gas component in the parameter range of *T*<sub>kin</sub> ≈ 167<sup>+153</sup><sub>-56</sub> K, log(*n*(*H*₂)/cm<sup>-3</sup>) ≈ 5.11<sup>+1.83</sup><sub>-0.58</sub>, and log(*N*(*CO*)/cm<sup>-2</sup>) ≈ 17.28<sup>+0.33</sup><sub>-0.42</sub> (in the one-component model) or *T*<sub>kin</sub> ≈ 306<sup>+263</sup><sub>-149</sub> K, log(*n*(*H*₂)/cm<sup>-3</sup>) ≈ 5.22<sup>+1.04</sup><sub>-0.49</sub>, and log(*N*(*CO*)/cm<sup>-2</sup>) ≈ 15.29<sup>+1.34</sup><sub>-1.17</sub> (in the two-component model). We are not able to rule out a “cold” component that dominates the molecular gas mass but barely contributes to the *J* ≥ 5 CO fluxes.

2. The CO SLED of J2310+1855 shows higher excitation compared to local/high-redshift starburst samples and local AGNs. Such high CO excitation is also found in other *z* ≳ 6 quasars (e.g., J1148+5251, J0100+2802, J0439+1634), and lensed high-redshift quasars (e.g., APM 08279+5255, the Cloverleaf).

3. The *L*<sub>H₂O</sub>/*L*<sub>IR(galaxy)</sub> ratio in this quasar is higher than local and high *z* (U/Hy) LIRGs. The luminous detections of H₂O and OH⁺ (high H₂O/CO and OH⁺/H₂O ratios) are suggesting other heating and ionization sources (e.g., cosmic-rays, shocks, and X-rays) in addition to PDR. In the nuclear region, the luminous quasar and the starburst activity are able to impact on the molecular gas through all of these possible mechanisms.

Complete measurements of the CO SLED of the quasar hosts at
*z* ≥ 6 are of great importance for our understanding of the physical
conditions and the heating mechanisms of the molecular gas in a
complex environment with both AGN and nuclear starburst
activity. It is also essential for higher-resolution observations to
map the distributions and kinematics of the highly excited
molecular gas around the AGN.

This work was supported by the National Science Foundation of
China (11721303, 11991052) and the National Key R&D program
of China (2016YFA0400702). We acknowledge Andres Asensio
Ramos for his substantial help on the code MOLPOP-CEP.
We thank Chentao Yang, Roberto Decarli, Toshiki Saito, and
Daizhong Liu for important discussions. R.W. acknowledges
support from the Thousand Youth Talents Program of China, and
NSFC grants No. 11533001, and 11473004. D.R. acknowledges
support from the National Science Foundation under grant
numbers AST-1614213 and AST-1910107 and from the Alex-
ander von Humboldt Foundation through a Humboldt Research

Fellowship for Experienced Researchers. Y.G.'s research is
supported by the National Key Basic Research and Development
Program of China (grant No. 2017YFA0402704), the National
Natural Science Foundation of China (grants No. 11861131007
and 11420101002), and the Chinese Academy of Sciences Key
Research Program of Frontier Sciences (grant No. QYZDJSSW-
SLH008). F.B. acknowledges support through the DFG Colla-
borative Research Centre 956. This paper is based on ALMA
observations: ADS/JAO.ALMA 2015.1.01265.S. ALMA is a
partnership of ESO (representing its member states), NSF (USA),
and NINS (Japan), together with NRC (Canada), MOST 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. This paper also makes
use of observations of IRAM NOEMA Interferometer: Project
number W18EE. IRAM is supported by INSU/CNRS (France),
MPG (Germany), and IGN (Spain).

ORCID iDs

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

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

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

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

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

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

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

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

Chris L. Carilli https://orcid.org/0000-0001-6647-3861

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

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

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

Desika Narayanan https://orcid.org/0000-0002-7064-4309

References

Apostolovski, Y., Aravena, M., Anguita, T., et al. 2019, A&A, 628, A23

Asensio Ramos, A., & Elitzur, M. 2018, A&A, 616, A131

Bañados, E., Decarli, R., Walter, F., et al. 2015, ApJL, 805, L8

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

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

Blake, G. A., Sutton, E. C., Masson, C. R., & Phillips, T. G. 1987, ApJ,
315, 621

Bothwell, M. S., Aguirre, J. E., Aravena, M., et al. 2017, MNRAS, 466, 2825

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. L., & Walter, F. 2013, ARA&A, 51, 105

Carniani, S., Gallerani, S., Vallini, L., et al. 2019, MNRAS, 489, 3939

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

Daddi, E., Dannerbauer, H., Liu, D., et al. 2015, A&A, 577, A46

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

Elitzur, M., & Asensio Ramos, A. 2006, MNRAS, 365, 779

Falgarone, E., Zwaan, M. A., Godard, B., et al. 2017, Natur, 548, 430

Feruglio, C., Fiore, F., Carniani, S., et al. 2018, A&A, 619, A39

Fixsen, D. J., Bennett, C. L., & Mather, J. C. 1999, ApJ, 526, 207

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP,
125, 306

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

González-Alfonso, E., Fischer, J., Aalto, S., & Falstad, N. 2014, A&A,
567, A91

González-Alfonso, E., Fischer, J., Bruderer, S., et al. 2018, ApJ, 857, 66

González-Alfonso, E., Fischer, J., Isaak, K., et al. 2010, A&A, 518, L43

Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and
Computational Science, 5, 65

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 217, ed. J. G. Mangum & S. E. J. Radford (San Francisco, CA: ASP), 299

Habart, E., Dartois, E., Abergel, A., et al. 2010, A&A, 518, L116

Hashimoto, T., Inoue, A. K., Tamura, Y., et al. 2018, arXiv:1811.00030

Jarugula, S., Vieira, J. D., Spilker, J. S., et al. 2019, ApJ, 880, 92

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

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

Jin, S., Daddi, E., Magdis, G. E., et al. 2019, ApJ, 887, 144

Kamenetzky, J., Glenn, J., Rangwala, N., et al. 2012, ApJ, 753, 70

Kollatschny, W., Weilbacher, P. M., Ochmann, M. W., et al. 2019, arXiv:1910.12813

Kormendy, J. 1993, in IAU Symp. 153, Galactic Bulges, ed. H. DeJonghe & H. Jan Habing (Dordrecht: Kluwer), 209

Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511

Krips, M., Martín, S., Eckart, A., et al. 2011, ApJ, 736, 37

Leipski, C., Meisenheimer, K., Walter, F., et al. 2013, ApJ, 772, 103

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

Maiolino, R., Cox, P., Caselli, P., et al. 2005, A&A, 440, L51

Maiolino, R., Gallerani, S., Neri, R., et al. 2012, MNRAS, 425, L66

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, adass XVI, 376, 127

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

Meijerink, R., Spaans, M., Loenen, A. F., et al. 2011, A&A, 525, A119

Neeleman, M., Banados, E., Walter, F., et al. 2019, ApJ, 882, 10

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

Odorico, D', Feruglio, V., & Ferrara, C. 2018, ApJL, 863, L29

Omont, A., Yang, C., Cox, P., et al. 2013, A&A, 551, A115

Panuzzo, P., Rangwala, N., Rykala, A., et al. 2010, A&A, 518, L37

Papadopoulos, P. P., van der Werf, P. P., Xilouris, E. M., et al. 2012, MNRAS, 426, 2601

Papadopoulos, P. P., Zhang, Z.-Y., Weiss, A., et al. 2013, in IAU Symp. 292, Molecular Gas, Dust, and Star Formation in Galaxies (Cambridge: Cambridge Univ. Press), 209

Pellegrini, E. W., Smith, J. D., Wolfire, M. G., et al. 2013, ApJL, 779, L19

Putaud, T., Michaut, X., Le Petit, F., et al. 2019, A&A, 632, A8

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., Leung, T. K. D., Ivison, R. J., et al. 2017, ApJ, 850, 1

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

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

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., Jones, G. C., et al. 2017, ApJ, 845, 138

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

Stanley, F., Jolly, J. B., König, S., et al. 2019, A&A, 631, A78

Uzgil, B. D., Bradford, C. M., Hailey-Dunsheath, S., Maloney, P. R., & Aguirre, J. E. 2016, ApJ, 832, 209

van der Werf, P. P., Berciano Alba, A., Spaans, M., et al. 2011, ApJL, 741, L38

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

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

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

Venemans, B. P., Walter, F., Decarli, R., et al. 2017b, ApJL, 851, L8

Viti, S., García-Burillo, S., Fuente, A., et al. 2014, A&A, 570, A28

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

Volonteri, M. 2012, Sci, 337, 544

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

Walter, F., Carilli, C., Bertoldi, F., et al. 2004, ApJL, 615, L17

Walter, F., Riechers, D., Cox, P., et al. 2009, Natur, 457, 699

Walter, F., Riechers, D., Novak, M., et al. 2018, ApJL, 869, L22

Walter, F., Riechers, D. A., Carilli, C. L., et al. 2007, in ASP Conf. Ser. 375, From Z-Machines to ALMA: (Sub)Millimeter Spectroscopy of Galaxies, ed. A. J. Baker et al. (San Francisco, CA: ASP), 182

Wang, F., Wang, R., Fan, X., et al. 2019a, 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. 2019b, ApJ, 887, 40

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

Wang, R., Wagg, J., Carilli, C. L., et al. 2011b, 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., Requena-Torres, M. A., Güsten, R., et al. 2010, A&A, 521, L1

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

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

Zhang, Z.-Y., Papadopoulos, P. P., Ivison, R. J., et al. 2016, RSOS, 3, 160025