Diverse Molecular Gas Excitations in Quasar Host Galaxies at z ∼ 6

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

Active star formation has been detected in the host galaxies of quasars from the local to the distant universe, revealing co-evolution of the supermassive black holes (SMBHs) and galaxies (e.g., Kormendy & Ho 2013;Wang et al. 2013Wang et al. , 2016;;Willott et al. 2015;Venemans et al. 2016;Decarli et al. 2017;Pensabene et al. 2020).The cold (∼ 100 K) interstellar medium (ISM) in the quasar host galaxies fuels both the active galactic nucleus (AGN) and the nuclear star formation.Obser-Li et al.
vations of the molecular and atomic lines at submillimeter/millimeter wavelengths provide rich information about the physical conditions and heating power of the cold ISM, which are the keys to understanding the evolution of the young quasar hosts (e.g., Wang et al. 2011Wang et al. , 2013Wang et al. , 2016;;Decarli et al. 2017Decarli et al. , 2018Decarli et al. , 2023;;Venemans et al. 2018;Shao et al. 2019Shao et al. , 2022;;Li et al. 2020a,b;Meyer et al. 2022).
The CO molecule, as the second most abundant species of the molecular ISM, gives the brightest molecular line emission in the (sub)millimeter regime.The CO spectral line energy distributions (SLEDs), i.e., the CO line fluxes as a function of rotational transition quantum number J, are sensitive diagnostics of the physical properties of the molecular gas (e.g., temperature, density, and radiation field strength).The low-J (J ≤ 3) CO emission lines directly constrain the total molecular gas mass of galaxies (e.g., Walter et al. 2003;Riechers et al. 2006Riechers et al. , 2011a;;Wang et al. 2011;Shao et al. 2019;Ramos Almeida et al. 2022;Montoya Arroyave et al. 2023).Low-J CO lines are found to be sublinearly correlated with far-infrared (FIR) luminosity (e.g., Riechers et al. 2006;Riechers 2011b;Kamenetzky et al. 2016).The mid-J (4 ≤ J ≤ 8) CO lines are reported to be linearly correlated with FIR luminosity (e.g., Greve et al. 2014;Lu et al. 2014;Liu et al. 2015;Kamenetzky et al. 2016), which is consistent with the gas heated by the far-ultraviolet (FUV; 6 eV < hν < 13.6 eV) photons from young-massive stars (Carilli & Walter 2013).The high-J (J ≥ 9) CO lines require high temperature and density to be excited, and processes such as X-ray (∼ 1−100 keV) heating from the AGN, mechanical heating from shocks, and cosmic-ray heating are frequently proposed to explain their excitation (e.g., Weiß et al. 2007;Spinoglio et al. 2012;Meijerink et al. 2013;Li et al. 2020a).Thus, the CO SLEDs also serve as probes for discriminating between different gas heating scenarios (Meijerink & Spaans 2005;Meijerink et al. 2007;Spaans 2008), e.g., the photo-dissociation region (PDR), where the gas is heated by the far-ultraviolet (FUV) photons from young massive stars and the AGN-heated X-ray dominated region (XDR).
Searching for direct evidence of the influence of the AGN activity on the ISM excitation is observationally challenging even in local/low-z well-studied samples.Valentino et al. (2021) find that for a sample of normal z ∼ 1 − 1.7 galaxies that host AGNs, marginal correlations have been found between the presence and strength of AGNs and the CO excitation on galaxy scales.On nuclear scales of 250 pc, Esposito et al. (2022) report weak correlations between the CO excitation and either the Xray or FUV flux in a sample of 35 local active galaxies, which suggests that neither the AGN nor star formation has a dominant effect on the gas excitation in the centers of these galaxies.Studies of CO SLEDs of individual local AGN host galaxies suggest a mix of a PDR component that contributes to the low-to-mid J CO fluxes and an XDR that dominates the high-J CO emission (e.g., van der Werf et al. 2010;Spinoglio et al. 2012;Pozzi et al. 2017;Mingozzi et al. 2018).Other mechanisms, such as mechanical heating by shocks, dense PDRs, and cosmic-ray heating, are also reported to dominate the high-J CO excitation in some local AGN hosts and starburst galaxies (e.g., van der Werf et al. 2010;Spinoglio et al. 2012;Pellegrini et al. 2013).
Luminous quasars at high-z represent the most luminous AGNs across redshifts.They are reported to be responsible for their high CO excitation at high-J.The lensed quasar APM 08279 at z = 3.9, is currently the highest excited CO SLED ever published across redshifts.Its CO emission is detected within the central ∼ 550 pc scale (Riechers et al. 2009), and the high CO excitation is best explained by the X-ray heating from the AGN (Weiß et al. 2007;Bradford et al. 2011).Another example is the Cloverleaf quasar at z = 2.6, for which the CO SLED is best fitted with a combination of a PDR and an XDR component (Uzgil et al. 2016).
Observations of quasars at z ∼ 6 from optical to (sub-) millimeter wavelengths reveal that ∼ 1/3 of them reside in extremely active environments with rapid SMBH accretion at rates of Ṁ BH > 10 M ⊙ yr −1 and intense star formation at rates up to a few ∼ 1000 M ⊙ yr −1 in the host galaxy (e.g., Wang et al. 2008;Carilli & Walter 2013;De Rosa et al. 2014;Venemans et al. 2020).In such systems, both the star formation and the SMBH activity lead to the injection of vast amounts of energy into the ISM of the host galaxy, making its ISM an ideal target to study the early co-evolution between the SMBHs and their host galaxies.Well-sampled CO SLEDs (with at least four transitions detected) have been obtained for five quasars at z ∼ 6 (Bertoldi et al. 2003;Walter et al. 2003;Riechers et al. 2009;Gallerani et al. 2014;Carniani et al. 2019;Wang et al. 2019;Yang et al. 2019;Li et al. 2020a;Pensabene et al. 2021).Notably, all of them show highly-excited CO SLEDs with much brighter high-J CO lines compared to that of the local starburst galaxies and AGNs, suggesting that the high-J CO emission is likely dominated by an XDR component (Gallerani et al. 2014;Wang et al. 2019;Yang et al. 2019;Li et al. 2020a;Pensabene et al. 2021).However, the shapes of the CO SLEDs are diverse from object to object.For example, the CO SLED of J0100+2802 is double-peaked, suggesting more than one gas component (Wang et al. 2019).In contrast, the CO SLED of J2310+1855 reveals one dominant gas component for the mid and high-J CO transitions (Li et al. 2020a).
In this paper, we use Northern Extended Millimeter Array (NOEMA) operated by the Institut de Radioastronomie Millimétrique (IRAM) to sample the CO SLEDs of the quasar P215−16 at z = 5.78 and the quasar J1429+5447 at z = 6.18, each with at least five transitions ranging from CO(5 − 4) to CO(10 − 9).P215−16 is among the optically brightest quasars at z ∼ 6, first identified in the Pan-STARRS1 survey (Morganson et al. 2012).It is also the source with the brightest 450 and 850 µm flux densities in our SCUBA2 High rEdshift bRight quasaR surveY (SHERRY) of 5.6 < z < 6.9 quasars (Li et al. 2020).J1429+5447 is among the most distant radio-loud quasars at z > 6, discovered in the Canadian-French-Hawaii Quasar Survey (CFHQS; Willott et al. 2010).It is also among the most X-ray luminous quasar observed by eROSITA (Medvedev et al. 2020;Khorunzhev et al. 2021;Migliori et al. 2023), and the second (sub)millimeter brightest and the [C II] 158µm brightest radio-loud quasar published at z > 6 ( Khusanova et al. 2022).
The structure of this paper is as follows.We describe observations and data reduction in Section 2 and observational results in Section 3. In Section 4, we derive the physical conditions of the dust and the gaseous ISM using dust spectral energy distribution models and radiative transfer models, respectively.Discussions of the gas properties and comparison with other galaxy samples are presented in Section 5. We adopt a flat ΛCDM cosmology model with H 0 = 70 km s −1 Mpc −1 and Ω M = 0.3, so that 1 ′′ corresponds to 5.8 kpc and 5.6 kpc at z = 5.78 and z = 6.18, respectively.
We reduced the data using the CLIC and MAP-PING packages that are part of the Grenoble Image and Line Data Analysis System (GILDAS) software (Guilloteau & Lucas 2000).We obtained the continuum uv tables by averaging all the line-free channels using the UV AVERAGE task in MAPPING.To generate the continuum subtracted spectral line uv tables, we adopted the continuum uv table obtained above as a continuum model, and use the UV SUBTRACT task to subtract it from the original uv tables.The uv tables of the spectral lines and the continuum were cleaned and imaged with natural weighting to maximize the S/N ratio (SNR) for point source.The host galaxies are spatially unresolved in beam sizes of ∼ 3 ′′ -6 ′′ × 1 ′′ -2 ′′ for P215−16 and 2 ′′ -7 ′′ × 2 ′′ -5 ′′ for J1429+5447, respec- -4000 -3000 -2000 -1000 0 1000 2000 Velocity (km s 1 ) tively.The spectra were finally binned to 46-78 km s −1 for the quasar P215−16 and 45-75 km s −1 for the quasar J1429+5447.The rms noise for the continuum emission in all bands is 0.02-0.05mJy beam −1 for both sources.Detailed information on the targeted lines and derived parameters can be found in Table 1.

RESULTS
We detect the CO(5 − 4), CO(6 − 5), CO(8 − 7), CO(9 − 8), CO(10 − 9), and the H 2 O(3 21 − 3 12 ) lines from the host galaxy of P215−16 at SNR > 5, and report a non-detection of the OH + (1 1 − 0 1 ) line (Fig. 1).A tentative detection (∼ 2.3σ) of the H 2 O(3 12 − 2 21 ) line blended with the CO(10 − 9) line is also obtained for this object.The continuum emission is detected at SNR ratio > 10 in all the observed bands (Table 1).We extract the spectra from the peak pixels and fit the spectra of CO(5 − 4), CO(6 − 5), CO(8 − 7), and CO(9 − 8) lines separately, each with a Gaussian profile.The redshifts and line widths derived for different CO lines are consistent with each other within the uncertainties.The CO(10−9), H 2 O(3 12 −2 21 ), and H 2 O(3 21 −3 12 ) lines are covered in the same sideband.We fit these three lines simultaneously, using a Gaussian model, with the line centers fixed to the redshift derived from the CO(9 − 8) line and an identical line width for all of the three spectral lines (assuming that these three lines are originating from the same gas component).The spectra and the Gaussian fitting results are shown in Fig. 1.Adopting the line width of the CO(9 − 8) line, we obtain a velocity-integrated map centered at the OH + (1 1 − 0 1 ) line frequency and measure the 3σ detection upper limit for the OH + (1 1 − 0 1 ) line.The measurements for both spectral lines and dust continuum are presented in Table 1.
For the quasar J1429+5447, we detect the CO(5 − 4), CO(6 − 5), CO(7 − 6), and CO(9 − 8) lines at SNR > 4, and tentatively detect the CO(10 − 9) and the H 2 O(3 21 − 3 12 ) lines at the 3.5 and 2.5σ levels, respectively.The [C I](2−1), H 2 O(3 12 −2 21 ), and OH + (1 1 −0 1 ) lines are undetected (Fig. 1).The CO lines detected in J1429+5447 has lower S/N compared to those detected in P215−16.The SNR of the CO(7 − 6) line spectrum is low.We fix the Gaussian line width and redshift to that of the CO(6 − 5) line and fit the line flux.The resulting line widths, redshifts, fluxes, and continuum fluxes are listed in Table 1.The 3σ upper limits of the [C I](2 − 1), OH + (1 1 − 0 1 ), and H 2 O(3 12 − 2 21 ) lines are estimated based on the rms of the moment 0 maps, which are integrated over channels expected for these lines assuming the same line widths and redshifts as that of their neighboring CO lines.The derived redshifts and line widths of the CO lines are consistent with each other within the uncertainties.These are also consistent with the fitting results of the W component of the CO(2 − 1) line (the western peak that coincides with the optical quasar) reported in Wang et al. (2011).The continuum emissions are all detected at a > 4σ level.More details about the measured dust and spectral line properties are listed in Table 1.

Dust SED fitting
To derive the dust properties of the two quasars, we analyze their dust spectral energy distributions (SEDs).Optically thick dust model has been suggested to explain the dust emission in a few high-z starburst systems, especially at rest-frame wavelengths <∼ 200 µm (e.g., Riechers et al. 2013;Wang et al. 2019).Since the quasars in this work are all far-infrared bright sources, it is reasonable to consider that the dust can be optically thick close to the peak of the dust SED.We thus employ the general opacity dust formula published in da Cunha et al. ( 2021) to fit the dust SEDs of the quasars in this work.In the general opacity regime, the dust emission can be described by: taking into account the effect of the cosmic microwave background (CMB).Here B ν [T dust (z)] or B ν [T CMB (z)] is the Planck function at the temperature of the dust or CMB at the quasar's redshift, and the dust optical depth τ ν is proportional to the dust surface density Σ dust via: where κ ν is the dust opacity.Its dependence on frequency is: where β is the dust emissivity index and κ 0 is the emissivity of dust grains per unit mass at frequency ν 0 .If we assume a simple spherical geometry, then the dust mass is: where R is the dust radius of a galaxy.In da Cunha et al. ( 2021), an additional parameter λ thick is introduced, which is defined when the optical depth reaches 1,  Accordingly, the dust emission can be rewritten as: where c is the speed of light.
Complemented with SCUBA 2 observations of the dust continuum at 450 µm (666.2GHz) and 850 µm (352.7 GHz) from Li et al. (2020), we use the general opacity model (Eq.6) to fit the observed dust SED of P215−16, normalized to the continuum flux density at 666.2 GHz.Under the general opacity assumption, the dust SED is determined by three parameters: namely T dust , β, and λ thick .We consider three different dust models.To break the degeneracy between different parameters in dust SED fitting, we fix λ thick to (1) the average value of SMGs (λ thick = 100 µm; da Cunha et al. 2021) (Model 1), and (2) the average value of QSOs at z ∼ 6 (λ thick = 264 µm; derived from Gilli et al. 2022) (Model 2), while enabling T dust and β to vary in these two models.In a third model, we enable all three param-eters λ thick , T dust , and β to vary (Model 3).To explore the posterior distributions of the parameters, we use the Markov Chain Monte Carlo (MCMC) python package emcee to fit the data to the model (Foreman-Mackey et al. 2013), adopting the best-fit values obtained for the least square estimation as educated initial guesses of the parameters.
J1429+5447 is a radio-loud quasar, we thus also include a power-law radio emission component f ν ∝ ν α , with a power-law index of α to fit the data in addition to the thermal dust SED model (Eq.6).Besides the continuum emission detected in this paper, we include published continuum data at 1.4, 3 GHz from the Very Large Array Sky Survey (VLASS) and Faint Images of the Radio Sky at Twenty Centimeters (FIRST) (Villarreal Hernández & Andernach 2018; Becker et al. 1995), and 32 and 250 GHz obtained using the Karl Jansky Very Large Array (VLA), and NOEMA (Wang et al. 2011;Khusanova et al. 2022), respectively.The observational data cannot constrain all the model parameters for the radio power-law + general opacity model.We thus fix λ thick = 100µm (typical of SMGs) and explore two cases with different assumptions of T dust : namely, 50 K (Model 1) and 100 K (Model 2) representing a typical "cold" and "warm" dust, respectively1 .
The derived dust properties of P215−16 and J1429+5447 are presented in Table 2. Future highfrequency continuum observations, possibly sampling the peak of the dust SED and/or the Wien tail will be critical to further constrain the dust temperature as well as the IR luminosity.

ISM properties
CO SLEDs are sensitive probes of the physical conditions (e.g., temperature, density, and radiation field) of the molecular gas.To derive the physical properties of the molecular gas in the host galaxies of these two quasars, we employ the CO SLED grid models published in Pensabene et al. (2021).These models consider two different gas heating mechanisms, namely a PDR heated by the FUV photons from young massive stars and an XDR heated by the X-ray photons from the AGNs.In the PDR model, the shape of the CO SLED is described by two physical parameters: the gas volume density n(H 2 ) and the FUV radiation field strength G 0 in Habing field units (1 Habing field corresponds to 1.6 × 10 −3 erg s −1 cm −2 ; Habing 1968).The CO SLED in the XDR model is determined by the gas volume density n(H 2 ) and the strength of the X-ray radiation field F x in the unit of erg s −1 cm −2 .To fully sample the molecular region which has a typical column density of >2 × 10 22 cm −2 (McKee & Ostriker 2007), we fix the column density of molecular hydrogen to 10 23 cm −2 for both sets of the PDR and XDR models.The adopted column density is also consistent with the values derived/adopted for other z∼6 quasars (e.g., Meyer et al. 2022;Decarli et al. 2023).Ranges of n(H 2 ), G 0 , and F x for typical molecular clouds explored in the model grids are 10 2 −10 6 cm −3 , 10 1 −10 6 and 10 −2 −10 2 erg s −1 cm −2 (Pensabene et al. 2021).We consider two cases for the CO SLED fitting: 1) a single PDR component and 2) an XDR component.The fitting procedure for each case is as follows.Firstly, we use the least square estimation for a first-order estimation of the best-fit values of the parameters.Secondly, we employ the MCMC package emcee to sample the posterior distribution of the parameters while using the best-fit values derived above as educated initial guesses for the parameters (Foreman-Mackey et al. 2013).
We show the fitting results in Fig. 4. The CO SLED of P215−16 cannot be fitted with any of our PDR models with n(H 2 ) and G 0 in the range of 10 2 − 10 6 cm −3 and 10 1 − 10 6 (Fig. 4      the PDR regions in star-forming galaxies, starbursts and AGNs (e.g., Gallerani et al. 2014;Pereira-Santaella et al. 2014;Pensabene et al. 2021;Esposito et al. 2022).To further explore if the data can be fitted with an extreme PDR model, we extend the range of n(H 2 ) and G 0 parameters of Pensabene et al. (2021) to 10 2 − 10 8 cm −3 and 10 1 − 10 8 , respectively, while using the same set of assumptions in CLOUDY.The results suggest that the data can be fitted with an extreme PDR of n(H 2 ) = 10 7.1 cm −3 and G 0 = 10 7.8 (Fig. 4(b)).The CO SLED of P215−16 can also be fitted with an XDR model.The least square approximation finds a dense XDR component of n(H 2 ) = 10 4.3 cm −3 illuminated by a radiation field of F x = 10 1.1 erg s −1 cm −2 .The MCMC analysis finds n(H 2 ) = 10 4.7 +0.8 −0.6 cm −3 and F x = 10 0.87 +0.6 −0.4 erg s −1 cm −2 for the XDR model, where values and uncertainties are estimated based on the 50th, 16th, and 84th percentiles of the samples in the marginalized distributions of the parameters.It is also possible that the CO SLED of P215−16 is contributed by both PDR and XDR components.However, the current CO SLED measurements are insufficient to constrain the parameters for a PDR+XDR model.Future observations of the high-J (J ≥ 11) CO lines will be critical to discriminate between different models.As for J1429+5447, the combination of available  7−4.7 adopting an L IR value of 6.89 +0.18  −0.18 ×10 12 L ⊙ − 3.68 +0.42  −0.42 × 10 13 L ⊙ .Detections of the CO emission lines complemented with previous CO(2 − 1) observations (Wang et al. 2011) enable us to constrain, for the first time, the molecular ISM physical conditions of a radio-loud quasar at z > 6.
We analyze the CO SLED of J1429+5447 following a similar procedure as that of P215−16.A single "dense" PDR model illuminated by an intense FUV radiation field is able to reproduce the observational data, with parameters of n(H 2 ) = 10 6.0 cm −3 , G 0 = 10 5.4 .A single XDR model is also able to reproduce the CO SLED with best-fit parameters of n(H 2 ) = 10 3.7 cm −3 , F x = 10 1.5 erg s −1 cm −2 .Using the gas conditions derived from [C II] 158µm , [C I](2 − 1) and L IR as a constraint of the PDR component, we also try to fit the CO SLED of J1429+5447 with a PDR + XDR model.This leads to a best-fit result of n(H 2 ) = 10 4.3 cm −3 and F x = 10 2.0 erg s −1 cm −2 for the XDR component, and n(H 2 ) = 10 5.1 cm −3 and G 0 = 10 4.2 for the PDR component.The PDR component contributes to 96% of the total molecular gas mass in this case.The MCMC method finds bestfit values in the range of n(H 2 ) = 10 5.8 +0.7 −0.3 cm −3 and G 0 = 10 5.6 +0.4 −0.7 (PDR model), n(H 2 ) = 10 3.8 +1.1 −0.3 cm −3 and F x = 10 0.33 +1.4 −0.7 erg s −1 cm −2 (XDR model), and n(H 2 ) = 10 5.1 +0.7 −1.1 cm −3 and G 0 = 10 4.2 +0.4 −0.5 , and n(H 2 ) = 10 4.4 +0.9 −0.7 cm −3 and F x = 10 1.3 Results of the gas properties derived from CO SLED modeling are shown in Table 2.In short, both the extreme-PDR and the XDR model could explain the observed CO SLEDs (up to J = 10) of P215−16.However, we also notice that the extreme-PDR model for P215−16 requires extremely high gas density and FUV radiation field intensity which is rarely seen in starburst systems and may suggest AGN contribution.For J1429+5447, the CO SLED could be fitted either with a PDR, XDR, or PDR + XDR model.More discussions about the excitation mechanism of the molecular gas can be found in Section 5.
In and CO in the molecular ISM phase, traces obscured star formation in dense molecular clouds.Linear relations were reported between the mid/high-J H 2 O and the IR luminosities in local and high-z infrared bright galaxies and AGNs (e.g., Yang et al. 2013Yang et al. , 2016;;Pensabene et al. 2022).The detections and upper limits of H 2 O(3 12 − 2 21 ) and H 2 O(3 21 − 3 12 ) for P215−16 and J1429+5447 follow the H 2 O−IR trend of other z ∼ 6 quasars (Yang et al. 2019;Pensabene et al. 2022;Riechers et al. in prep).

DISCUSSION: THE HEATING MECHANISMS OF THE MOLECULAR GAS
With the NOEMA observations of the multiple CO transitions, we analyze the CO SLEDs of two z > 5.7 quasars.Both of them are amongst the most infrared luminous objects found so far in the early Universe.The IR luminosities are of order of a few 10 12 to 10 13 L ⊙ , revealing massive star formation in a dusty ISM.J1429+5447 is a radio-loud quasar, which enables us to study the CO excitation close to a young and radio-loud AGN for the first time.
The CO SLEDs of local and high-z starburst galaxy samples typically peak at J ≲ 6, and drop rapidly at high J, resulting in low line fluxes at J > 10 (Fig. 6a).This is consistent with the gas heated by the FUV photons from young massive stars, indicating a PDR origin.On the other hand, the local AGNs and high-z quasars tend to have higher CO excitation with the CO SLEDs peaking at J ≳ 6.The bright high-J (J ≳ 9) CO transitions are usually associated with other heating mechanisms, e.g., X-ray heating from AGNs, cosmic ray heating, and mechanical heating by shocks in addition to the FUV heating from young massive stars (Fig. 6b  and c; e.g., van der Werf et al. 2010;Spinoglio et al. 2012).As is suggested by Fig. 6 and radiative transfer model predictions, the CO SLEDs of the starburst galaxies and AGNs/quasars are distinguishable through observations of the high-J (e.g., J ≳ 9) CO transitions (e.g., Vallini et al. 2019).
P215−16 exhibits a highly excited CO SLED indicating an excitation level that is much higher than that of all local/high-z starburst galaxies (Fig. 6a and b) and local AGNs (Fig. 6b), with no turn-over up to J = 10.It has the highest CO(10 − 9)/CO(6 − 5) ratio among all published quasars at z ∼ 6 (Fig. 6c) and the second highest excited CO ever across redshifts (less excited compared to APM 08279).We have demonstrated in Section 4.2 that either an XDR component with high gas density and intense radiation (n(H 2 ) = 10 4.3 cm −3 , F x = 10 1.1 erg s −1 cm −2 ) or an extreme PDR of n(H 2 ) = 10 7.1 cm −3 and G 0 = 10 7.8 could explain such a highly excited CO SLED (Fig. 4).Mechanical heating by shocks can also result in highly excited molecular CO which was reported in some local galaxies (Rosenberg et al. 2015).Shocks are relatively more efficient in heating the gas compared to the dust, and a typical signature of shock-heated gas is a high CO/IR luminosity ratio of a few times 10 −4 (e.g., Pellegrini et al. 2013).The CO/IR ratio (even for the brightest CO(10 − 9) transition) of P215−16 is ∼ 10 −5 , which disfavors the shock heating scenario.Enhanced cosmic ray heating could also be responsible for high-J CO line excitation which is sometimes indistinguishable from XDR. Vallini et al. (2019) investigated several CO excitation mechanisms in high-z galaxies and reported that the most extreme cosmic-ray heated CO SLED peaked at J = 9.The high CO excitation of P215−16, which peaks at J ≥ 10, is best explained by an XDR or extreme-PDR component.In either case, AGN could be the source powering the high CO excitation, by contributing to the X-ray emission or intense FUV emission as suggested by the results of the XDR and extreme-PDR model.
To the contrary, J1429+5447 exhibits a relatively lower CO excitation, which is comparable to that measured for local and high-z starburst galaxies (Fig. 6a  and b), and z ≳ 6 starburst galaxies, e.g., SPT 0311 W, SPT 0311 E, and HFLS3 (Fig. 6d).It has the lowest CO(J −J-1) [J≥7] /CO(6−5) line ratios among the z ∼ 6 quasars with published CO SLEDs (Fig. 6d).The CO SLED could be explained by a single PDR model with extreme gas densities of n(H 2 ) = 10 6.0 cm −3 and radiation of G 0 = 10 5.4 .
We also explore if the CO SLED of J1429+5447 can be fitted with an XDR model included, given the luminous X-ray detection in this quasar.The CO SLED of J1429+5447 fitted by an XDR or a PDR+XDR model prefers an XDR component with moderate density (n(H 2 ) = 10 3.7−4.3cm −3 ) and illuminated by an intense X-ray radiation field (F x = 10 1.5−2.0erg s −1 cm −2 ).This is at the lower edge of the molecular gas density range found in z ∼ 6 quasars (e.g., Li et al. 2020a;Pensabene et al. 2021).As J1429+5447 is a radio-loud object, it is possible that the radio jet expels the gas out of the host galaxy resulting in lower gas density.This is consistent with the indication of a possible outflow found in previous [C II] 158µm observations (Khusanova et al. 2022).The outflowing gas mass constrained from the broad [C II] 158µm component is 1.2 × 10 10 M ⊙ assuming an α [CII] = 30 M ⊙ L −1 ⊙ (Zanella et al. 2018).The OH + (1 1 − 0 1 ) absorption is a sensitive tracer of outflows in galaxies and quasars (e.g., Shao et al. 2022).Unfortunately, the low S/N of our OH + (1 1 − 0 1 ) spectrum prevents us from deriving useful constraints on the outflowing molecular gas mass and deeper observations will be needed to obtain an independent estimate.In this work, we only include a simple plane parallel geometry and include no dust torus attenuation.A low CO excitation can be produced in a sophisticated geometry XDR model of high dust attenuation even when the X-ray flux and gas density are high (e.g., Vallini et al. 2019).Future CO SLED observations of more radio-loud quasars will possibly enable us to study if the low excitation is rare/common in the radio-loud population.
Adopting a general opacity dust SED model, we derive dust temperatures in the range of ∼ 50 − 110 K for the two quasar host galaxies.The derived infrared luminosities of P215−16 and J1429+5447 differ by factors of ∼ 2.5 and 5.3, respectively, depending on the best-fit models used.
The quasar P215−16 exhibits highly excited molecular CO with a CO SLED peaking at J ≥ 10.It has the highest CO(J − J-1) [J≥7] /CO(6 − 5) ratio among the z ∼ 6 quasars with published CO SLEDs known to date.Analysis of the CO SLED with models suggests that a "dense" (n(H 2 ) = 10 4.3 cm −3 ) XDR component illuminated by an intense X-ray radiation field (F x = 10 1.1 ) or a "'dense" PDR component with n(H 2 ) = 10 7.1 cm −3 and G 0 = 10 7.8 reproduce the observed CO SLED.In either case, the central AGN is likely to contribute to the radiation field to excite the high−J CO transitions.
The CO SLED of the radio-loud quasar J1429+5447 reveals the lowest CO excitation among all published quasar CO SLEDs at z ∼ 6.This is the first reported CO SLED of a radio-loud quasar at this redshift and is comparable to those of local AGNs and starburst galaxies at z ≳ 6.The CO SLED could be explained by a PDR component.In addition, we also propose two possibilities where XDR is considered to explain the low CO excitation: 1) a luminous X-ray radiation field illuminating a moderate-density gas, which may be attributed to a result of the gas expelled by an outflow, and 2) a sophisticated geometry XDR model taking into account the X-ray attenuation by a dust torus.In this case, low CO excitation can be reproduced in a high dust attenuation model even when the X-ray fluxes and gas densities are high.
In addition, we detect bright H 2 O(3 21 − 3 12 ) emission in P215−16.This places P215−16 exactly on the linear relation between the H 2 O(3 21 − 3 12 ) and IR luminosities recently found for other quasars at z ∼ 6 and local/highz infrared bright galaxies.
Future observations of the CO SLEDs towards more radio-loud quasars at z ≳ 6, possibly up to (J ≳ 10) will be critical to systematically reveal the physical conditions and heating mechanisms of the molecular gas in the host galaxies of the radio-loud quasar population in

Figure 1 .
Figure 1.Spectra of the CO, H2O, OH + (11 − 01), [C I](2 − 1) lines for P215−16 (left) and J1429+5447 (right).The continuum emission has been subtracted for each of the spectra.The data are shown as yellow histograms.Gaussian profiles fitted to the line emission are shown as solid or dashed lines.The expected positions of the undetected lines are indicated as dashed blue vertical lines.In the spectra where more than one line is detected, the solid black lines represent the sum of different Gaussian components.

Figure 2 .
Figure 2. (a) Dust SED fitting results for P215−16.The data are shown as black squares.The Top and bottom axes show the observed and rest continuum frequency.The dust SED is fitted with three general opacity models, and the solid lines represent the models with best-fit values (see Sect. 4.1 for further details).(b-d) Posterior distributions of the model parameters for Model 3, 1 and 2, with the same color as the best-fit models in (a).The 2D contours show the [1, 2]σ confidence intervals in the marginalized distributions.

Figure 3 .
Figure 3. (a) Dust SED fitting results for J1429+5447.The data are shown as black squares.The Top and bottom axes show the observed and rest continuum frequency.The dust SED is fitted with the general opacity dust SED model + radio power-law emission (see Sect. 4.1 for further details).The solid/dashed blue (red) lines represent the best-fit general opacity/radio powerlaw models for Model 1(2).(b) and (c) Posterior distributions of the model parameters for P215−16 (c to e) and J1429+5447, respectively.The color scheme for different models is the same as in panel (a).The 2D contours show the [1, 2]σ confidence intervals in the marginalized distributions.
dust SED and the CO SLED fitting results."LS Best-fit" represents the best-fit result of the least square method."MCMC" represents the fitting result of the MCMC method.Model(λ thick−fixed ,T dust−fixed ) represents the model where λ thick or T dust is fixed to λ thick−fixed or T dust−fixed .As for the CO SLED modeling, n(H 2 ),G 0 , and Fx are in units of cm −3 , habing field, and erg s −1 cm −2 .
[C II] 158µm data with the information on the nondetection of the [C I](2 − 1) line from this work enables us to constrain the physical conditions of the atomic gas.The [C II] 158µm /[C I](2 − 1) luminosity ratio of > 23 argues for a PDR origin of the line excitation, as the XDR scenario usually results in a lower line ratio of ≲ 10 (Pensabene et al. 2021).Predic-tions of the [C II] 158µm /[C I](2 − 1) ratio in the PDR model are shown in Fig. 5 (Pensabene et al. 2021).The [C II] 158µm /[C I](2 − 1), [C II] 158µm /L IR , and [C I](2 − 1)/L IR ratios indicate a PDR component with parameters in the range of n(H 2 )= 10 3−6 cm −3 and G 0 = 10 3.

Figure 4 .
Figure 4. CO SLED fitting results of the PDR (left), extreme-PDR (middle), and XDR model (right) for P215−16 (top) and PDR (left), XDR (middle), and PDR+XDR (right) model for J1429+5447 (bottom).The fitting results are shown in solid blue (PDR) and red (XDR) lines, and the data are shown in black squares and diamonds for P215−16 and J1429+5447, respectively.The sum of the PDR+XDR model is shown as a solid black line.

Table 2 .
Dust SED and CO SLED modeling results +1.4  −1.3erg s −1 cm −2 (PDR + XDR model).Analysis of the gas properties using either the [C II] 158µm , [C I](2 − 1) and IR ratios or the CO SLED suggests that the [C II] 158µm , [C I](2−1), CO and IR emission can be explained with one PDR component.However, we are not able to rule out a single XDR model or a composite PDR + XDR model, which reproduces the observed CO SLED.In the case of either XDR or PDR + XDR model, a "moderate density" XDR component illuminated by an intense radiation field is suggested.