A comprehensive view of the interstellar medium in a quasar host galaxy at z~6.4

Characterizing the physical conditions (density, temperature, ionization state, metallicity, etc) of the interstellar medium is critical to our understanding of the formation and evolution of galaxies. Here we present a multi-line study of the interstellar medium in the host galaxy of a quasar at z~6.4, i.e., when the universe was 840 Myr old. This galaxy is one of the most active and massive objects emerging from the dark ages, and therefore represents a benchmark for models of the early formation of massive galaxies. We used the Atacama Large Millimeter Array to target an ensemble of tracers of ionized, neutral, and molecular gas, namely the fine-structure lines: [OIII] 88$\mu$m, [NII] 122$\mu$m, [CII] 158$\mu$m, and [CI] 370$\mu$m and the rotational transitions of CO(7-6), CO(15-14), CO(16-15), and CO(19-18); OH 163.1$\mu$m and 163.4$\mu$m; and H$_2$O 3(0,3)-2(1,2), 3(3,1)-4(0,4), 3(3,1)-3(2,2), 4(0,4)-3(1,3), 4(3,2)-4(2,3). All the targeted fine-structure lines are detected, as are half of the targeted molecular transitions. By combining the associated line luminosities, the constraints on the dust temperature from the underlying continuum emission, and predictions from photoionization models of the interstellar medium, we find that the ionized phase accounts for about one third of the total gaseous mass budget, and is responsible for half of the total [CII] emission. It is characterized by high density (n~180 cm$^{-3}$), typical of HII regions. The spectral energy distribution of the photoionizing radiation is comparable to that emitted by B-type stars. Star formation also appears to drive the excitation of the molecular medium. We find marginal evidence for outflow-related shocks in the dense molecular phase, but not in other gas phases. This study showcases the power of multi-line investigations in unveiling the properties of the star-forming medium in galaxies at cosmic dawn.


Introduction
Quasar host galaxies are among the first massive galaxies emerging from the dark ages (z ∼ > 6).They host black holes with masses that can exceed 10 9 M (e.g., Willott et al. 2003;De Rosa et al. 2011;Wu et al. 2015;Bañados et al. 2018).They often form stars at high rates (SFR=100-1,000 M yr −1 ; Bertoldi et al. 2003a;Walter et al. 2004Walter et al. , 2009a;;Wang et al. 2008aWang et al. ,b, 2013;;Leipski et al. 2014;Drouart et al. 2014;Venemans et al. 2018;Decarli et al. 2018).Immense gaseous reservoirs (M H2 ∼ > 10 10 M ) fuel this intense nuclear and star formation activity (e.g., Walter et al. 2003;Bertoldi et al. 2003b;Wang et al. 2008a;Venemans et al. 2017a).The luminosity of emission lines associated with heavy elements and ions in the broad-line region (Kurk et al. 2007;De Rosa et al. 2011, 2014;Schindler et al. 2020), as well as the presence of large reservoirs of dust in the interstellar medium (Wang et al. 2008a;Venemans et al. 2018) indicate that these galaxies are already highly metal-enriched.The mass, star formation rate, and metallicity of these early quasar host galaxies exceed by orders of magnitudes those of typical star-forming galaxies at z > 6 (see, e.g., Vanzella et al. 2014;Bouwens et al. 2015;Oesch et al. 2016;Harinake et al. 2017;Bañados et al. 2019;Ota et al. 2018;Salmon et al. 2018;Curti et al. 2022).In this respect, quasar host galaxies are arguably the most active and some of the most evolved objects emerging from the dark ages.They represent a challenge for models of early black hole and galaxy coexistence (Volonteri 2012;Habouzit et al. 2016Habouzit et al. , 2017;;Agarwal et al. 2017;Yue et al. 2017;Lupi et al. 2019Lupi et al. , 2022;;Romano 2023).Characterizing their early growth is there-fore a mandatory step to understand how first massive galaxies formed.
Direct observations of the starlight in quasar host galaxies at z ∼ > 6 are hindered by the unfavorable contrast against the bright nuclear emission from the quasar, and by a combination of redshift and dust obscuration that suppresses the observed restframe UV emission and shifts the rest-frame optical emission into the mid-infrared bands, where the atmosphere is opaque.These limitations also affect nebular gas tracers (the hydrogen Balmer lines: Hα, Hβ, Hγ, etc; ionized and neutral oxygen: [O iii] at 5008 Å, [O ii] at 3727 Å, [O i] at 6300 Å; ionized nitrogen [N ii] at 6584 Å; etc) and will only be effectively overcome by the James Webb Space Telescope.On the other hand, a wealth of information on the star-forming medium at z ∼ > 6 is accessible nowadays at sub-mm and mm wavelengths.Rest-frame farinfrared (IR) gas tracers detected in distant quasar host galaxies include various rotational molecular transitions from, e.g., carbon monoxide, CO; water, H 2 O; hydroxyl, OH (see, e.g., Walter et al. 2003;Weiß et al. 2007;Wang et al. 2011;van der Werf et al. 2011;Omont et al. 2011Omont et al. , 2013;;Riechers et al. 2013Riechers et al. , 2014;;Venemans et al. 2017a;Wang et al. 2019;Yang et al. 2019;Li et al. 2020;Pensabene et al. 2021;Decarli et al. 2022).In addition, many key elements and ions have fine-structure transitions at these wavelengths, such as ionized carbon [C ii] at 158 µm; neutral carbon [C i] at 370 µm and 609 µm; ionized nitrogen [N ii] at 122 µm and 205 µm; ionized oxygen [O iii] at 52 and 88 µm; and neutral oxygen [O i] at 63 µm and 146 µm.Many of these lines have been detected in high-redshift galaxies over the last few years (e.g., Maiolino et al. 2005Maiolino et al. , 2009;;Walter et al. 2009aWalter et al. , 2011Walter et al. , 2012;;Ferkinhoff et al. 2011Ferkinhoff et al. , 2015;;Coppin et al. 2012;Combes et al. 2012;Nagao et al. 2012;Decarli et al. 2012Decarli et al. , 2014;;Venemans et al. 2012Venemans et al. , 2017a,b;,b;Carilli & Walter 2013;Brisbin et al. 2015;Gullberg et al. 2015;Capak et al. 2015;Pavesi et al. 2016;Uzgil et al. 2016;Bothwell et al. 2017;Trakhtenbrot et al. 2017;Willott et al. 2017;Lamarche et al. 2017Lamarche et al. , 2019;;Carniani et al. 2018;Hashimoto et al. 2019;Novak et al. 2019;Tadaki et al. 2019;Boogaard et al. 2020;Valentino et al. 2020;Sugahara et al. 2021;Harrington et al. 2021;Meyer et al. 2022).Combinations of these emission lines can shed light on the gas mass in various phases of the interstellar medium (ISM; see, e.g., Ferkinhoff et al. 2015;Zanella et al. 2018;Dunne et al. 2021 and the review of Carilli & Walter 2013); on the star formation rate (SFR; see, e.g., De Looze et al. 2014;Herrera-Camus et al. 2015, 2018); on metallicity (Z; see, e.g., Nagao et al. 2012;Peng et al. 2021;Lamarche et al. 2022); on the gas density (n), on the origin of the excitation mechanism, and on the strength and hardness of the incident flux (for a recent review, see Wolfire et al. 2022).While mapping individual clouds in multiple ISM tracers at high redshift is impossible even with the excellent sensitivity and angular resolution provided by the Atacama Large Millimeter Array (ALMA), galaxy-averaged observations of these tracers can shed light on global properties of the ISM (e.g., on the mass budget in the ionized, neutral, and molecular phases; on the contribution of the central quasar to the ISM excitation; and on the total metallicity of the galaxy; see, e.g., Romano 2023).
In this work, we present a multi-line study of the starforming medium in the host galaxy of the quasar PSO J183.1124+05.0926 (hereafter, PJ183+05;R.A.: 12:12:26.97, Dec.: +05:05:33.4)at z = 6.4386.This quasar was discovered by Bañados et al. (2016) using color-color selections from the Pan-STARRS1 database (Chambers et al. 2016), and follow-up photometric and spectroscopic observations.The [C ii] luminosity in this source is the highest among 27 quasars at z > 6 surveyed in Decarli et al. (2018), and one of the highest from non-lensed sources currently known at any redshift (e.g., Pavesi et al. 2016;Tadaki et al. 2019;Andika et al. 2020;Mitsuhashi et al. 2021).The quasar resides at a redshift that is convenient for observations of various far-infrared emission lines that fall in transparent windows of the atmosphere at z ≈ 6.4.We targeted a suite of lines sampling the molecular, neutral and ionized components of the ISM in PJ183+05.The structure of the paper is as follows: In Sec. 2 we present the observations and the data processing.In Sec. 3, we analyze the dataset.In Sec. 4 we infer physical quantities from the observed spectra, and discuss our findings.Finally, we draw our conclusions in Sec. 5.
Throughout the paper we adopt a flat ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 , Ω m = 0.3 and Ω Λ = 0.7 (consistent with the measurements by the Planck Collaboration 2015).Within this cosmological framework, z = 6.4386 corresponds to a luminosity distance of D L = 62671 Mpc and an angular scale of 5.491 kpc per arcsec.

Observations and data reduction
The dataset used in this project consists of the [C ii] 158 µm observations of PJ183+05 presented in Decarli et al. (2018) (program ID: 2015.1.01115.S, PI: Walter); as well as data from a dedicated ALMA program (ID: 2016.1.00226.S, PI: Decarli).Observations were carried out in compact array configurations.Table 1 summarizes the observations.Integration times ranged between 8 and 80 min in each frequency setting.The quasars J1229+0203 and J1222+0413 served as bandpass / pointing and phase calibrators, respectively.Ganymede was observed for flux calibration, with the exception of the band 8 observations for which the quasar J1229+0203 served as flux calibrator.
Data were reduced and calibrated with the official ALMA pipeline in CASA (McMullin et al. 2007, version 4.7.2).We imaged the measurement sets using natural weighting, in order to maximize the signal-to-noise ratio of line detections.The [C ii] data presented in Decarli et al. (2018) was binned in 30 km s −1 channels.The other line observations presented here are binned in 90 km s −1 wide channels.The expected line width (∼ 375 km s −1 from the analysis of the [C ii] data) is thus sampled in ∼ 4 independent channels, while at the same time maximizing the signal-to-noise ratio in the lines.From each frequency setting, we create two cubes, corresponding to the lower and upper side bands.The cubes are cleaned to the 2-σ level, using cleaning masks on the quasar.
Two-dimensional channel-by-channel fits of the quasar emission at 88 µm in the rest-frame revealed that the emission is spatially unresolved or only marginally resolved at > 0.7 resolution (see also the analysis of the [C ii] emission presented in Venemans et al. 2020).Therefore, we extracted the spectra of the dust continuum and targeted lines in PJ183+05 on a singlepixel basis.The extracted spectra are shown in Fig. 1.

Spectral fits
We model the observed spectra as a flat continuum plus a Gaussian profile for each line.This is shown to provide a good description of the typical line profiles observed in z > 6 quasars (Decarli et al. 2018).For the fitting, we used our custom Monte Carlo Markov Chain tool to sample the posterior probability of each fitted parameter.Because of the limited S/N of some of the lines in our new data, we opt to assume a fixed Gaussian profile for each transition, with a line full width at half maximum of 375 km s −1 as observed for the [C ii] line (Decarli et al. 2018).The redshift of the line, its normalization, and the underlying continuum flux density are left free in the fit.As priors, we used a Gaussian with a standard deviation of 0.002 around the [C ii] redshift, corresponding to ∼ 80 km s −1 ; a Gaussian distribution for the line normalization (centered on a rough estimate of the line flux based on the line peak and assumed width); and a Maxwellian distribution for the continuum flux density, with a scale width based on the median value of the spectrum in each frequency tuning.
Figure 1 shows the spectra and their fits for all the frequency settings used in this work.We list the results of the fits in Tables 2  and 3.The dust continuum is clearly detected in all the frequency settings.We consider lines detected if the best fit of the line flux (as inferred from the median value of the posterior distribution) exceeds its 3-σ confidence level.For non-detections, we use a line flux upper limit based on the 3-σ confidence level of the posterior distribution, under the assumption of a fixed line width (375 km s −1 , based on [C ii]).All the fine structure lines in our study, as well as the CO(7-6) line, the OH doublet, and the H 2 O 3(3,1)-3(2,2) and 3(0,3)-2(1,2) transitions match this criterion, whereas the J up ≥ 15 transitions from CO, and the transitions involving the J up ≥ 4 levels from H 2 O remain undetected, and will be considered as 3-σ upper limits in our analysis.We note that, for the sake of internal consistency, we refit the [C ii] line with the same approach as all of the other lines analyzed in this paper, rather than referring to values reported in the literature.
We derive line luminosities as: where F line is the line flux, as measured by integrating over the fitted Gaussian profile, ν 0 is the rest-frame frequency of the transition, and D L is the luminosity distance; and Carilli & Walter 2013, for a discussion of line luminosity definitions).

Dust continuum modeling
A modified black body provides us with a good description of the dust continuum in PJ183+05.Under the assumption of a single dust temperature, T dust , the flux density at a given observed frequency ν = ν 0 /(1 + z) results from the black body emissivity, B(ν 0 , T dust ) observed against the cosmic microwave background (CMB) at the quasar's redshift, integrated over the source apparent area, Ω s , and corrected for the effects of radiative transfer: where: is the black body emissivity law, h is the Planck constant, c is the speed of light, and k b is the Boltzmann constant.The temperature of the CMB is T CMB = 2.725 (1+z) K = 20.27K at the quasar's redshift.The radiative transfer correction consists of a transmissive term and an absorbed term, where τ ν 0 is the optical depth: Here M dust is the dust mass seen within a solid angle Ω s ; D A = D L (1+z) −2 is the angular diameter distance; and κ(ν 0 ) is the dust emissivity law, for which we adopt an interpolation of the values reported in Table 5 of Draine (2003) at ν 0 > 1500 GHz, and a power-law extrapolation κ(ν 0 ) = 6.37 (ν 0 /1500 GHz) β cm 2 g −1 at lower frequencies.In this framework, the overall dust continuum is uniquely determined by a combination of four parameters: T dust , Ω s , M dust , and β, which we determine via a Monte Carlo Markov Chain fit.As priors, we adopt a log-normal distribution in Tdust = T dust − T CMB with scale parameter Tref =100 K −T CMB for the dust temperature (for reference, typical dust temperatures in quasar host galaxies at high redshifts are in the range T dust = 40 − 100 K; see, e.g., Walter et al. 2022).For the angular size, Ω s , we adopt a Gaussian distribution around the spatial extent reported in Venemans et al. (2020) (deconvolved dust sizes: 0.45 × 0.35 ≈ 2.4 × 1.9 kpc 2 ).The width of the prior distribution is set to 30% of the central value.We also adopt a log-normal distribution for the prior on the dust mass, centered around M dust = 10 8 M (based on the values typically found in z > 6 quasars, see, e.g., Venemans et al. 2018) and with a 1-σ width of 1 dex; and a log-normal distribution centered on β = 1.5 with a 1-σ width of 0.1 dex (see, e.g., Beelen et al. 2006).We find T dust = 47.0 +1.5 −2.0 K, Ω s = 0.155 +0.029 −0.022 arcsec 2 (≈ 4.67 +0.87  −0.66 kpc 2 ), log M dust /M = 8.94 +0.06 −0.05 , and β = 1.84 0.15 −0.16 .Fig. 2 shows the best-fit model.Fig. 3 shows the posterior distributions of the fitted parameters.The most identifiable correlation relates the dust mass and temperature, as expected in the optically-thin regime (from eq. 3 and 5, for τ 1).
We also infer posterior distributions for derived parameters, namely the optical depth at the [C ii] frequency, τ 1900 GHz , and the IR luminosity, L IR .Eq. 5 gives us τ 1900 GHz .The IR luminosity is computed by integrating the dust emission model in the range IR luminosity in the 8-1000 µm wavelength range.We find an optical depth at the frequency of [C ii] of τ 1900 GHz = 0.48 ± 0.04, and an IR luminosity of log L IR [L ]=12.98 +0.04  −0.03 .The posterior distributions for L IR and τ 1900 GHz are also shown in Fig. 3.

Emission line modeling
In order to put our observations in the context of the physical properties of the interstellar medium in PJ183+05, we derive analytical predictions based on first principles, and complement them with grids of predicted line ratios using the pho-toionization code Cloudy (Ferland et al. 1998(Ferland et al. , 2013(Ferland et al. , 2017)).We treat all the emission lines as arising from three distinct components (ionized, neutral, molecular), each described with a simple parametrization (e.g., constant density throughout each cloud).This is clearly a crude approximation that does not capture the complex structure of the ISM in galaxies.E.g., both spatial  emission in high-z galaxies (e.g., Carniani et al. 2018).However, these approximations suffice in order to infer luminosityweighted properties of the galaxy as a whole.

Analytical prescriptions from first principles
Analytical prescriptions and empirical scaling relations can provide us with a back-of-the-envelope approach to the radiative transfer problem.They rely on a number of simplifications and assumptions (e.g., photoionization equilibrium, local thermal equilibrium) that are likely not universally valid when one stud- ies the ISM of a galaxy in detail.However, this method allows for a straightforward interpretation of the observed line ratios in terms of physical quantities.Analytical implementations of the radiative transfer in ionized bubbles can work remarkably well, when compared with more sophisticated Cloudy models (see, e.g., Yang & Lidz 2020;Lamarche et al. 2017Lamarche et al. , 2022)).
In order to predict the luminosity of the ionized gas transitions in our study, we assume that a fully-ionized gas (i.e., volume density n = n e ) is in local thermodynamical and photoionization equilibrium, with T gas setting the energy distribution of the electrons following a Maxwell distribution.Free electrons are the main collision partners in a fully-ionized gas cloud.Bound electrons are excited to energy levels beyond ground state via collisions, and are de-excited via collisions and radiation.The population of level2 i is thus set by solving a system of linear equations of the form: Here R ul are the rates at which electrons are de-excited from a higher-energy level u to a lower-energy level l: We assume local thermal equilibrium, thus the number of collisions depends only on the gas temperature T gas and density n.The vast majority of electrons reside in the 3 P J levels.
where n γ is the photon occupation fraction, n γ = c 3 (8π hν 3 ) −1 u ν , and u ν is the energy volume density of the radiation field.On the other hand, R lu are the rates at which electrons are excited from l to u: Here A ul are the Einstein coefficients, g i are the statistical weights of the levels, k ul (T gas ) are the collision rates: and Ω(l, u) are the collision strengths, which we take from Appendix F of Draine (2011).The energy difference between two energy levels is E ul = hν ul .We assume that n γ is negligible.For 6-electron ions (C, N + , O ++ ), we consider a five-level structure where electrons populate the 3 P 0 (ground), 3 P 1 , 3 P 2 , 1 D 2 , or 1 S 0 levels, corresponding to g i =(2 J+1)=1, 3, 5, 5, and 1, respectively.For 5-electron ions (specifically, C + ), we only consider the two lower energy levels, 2 P 1/2 (ground) and 2 P 3/2 , as the next levels ( 4 P 1/2 , 4 P 3/2 , 4 P 5/2 ) are only significantly collisionally populated at T gas ∼ > 50000 K.
From our analysis, we infer the critical densities of the various levels, i.e., the densities at which collisional de-excitation equals radiative de-excitation: Tab. 4 lists the critical densities thus computed for the relevant transitions.
Solving the system of equations 6 as a function of n and T gas allows us to determine the fraction of the electron population in each level.This is shown in Fig. 4. The line luminosity associated with an optically-thin transition is then: Table 4. Critical densities for the main far-IR fine-structure lines from the ionized gas considered in our work, computed for different gas temperatures T gas .For [C ii], we consider the ionized component only.
Ion Line T gas [K] 5,000 10,000 20,000 where V is the total volume in Hii regions.Compared to the analysis of, e.g., Yang & Lidz (2020), we assume that all the ionizing photons contribute to the photoionization equilibrium, and that the sizes of the bubbles where O ++ , N + and H + reside are comparable.In sec. 4 we will discuss the impact of these assumptions.
We assume that photoionization sets the abundance of ions.For a given template of the ionizing source, we can derive the number of photons emitted by the photoionizing source per unit time that are energetic enough to ionize different elements and ions.Namely, for a species X N with N marking the ionization state (0=neutral, 1=first ionization, 2=second ionization, etc), we have that: and we define its net equivalent as the difference between the number of photons emitted per unit time that can photoionize X N into X N+1 minus the photons that can photoionize X N+1 into X N+2 : Fig. 5 compares the templates used in our radiative transfer analysis, and their yield in terms of flux of photons energetic enough to singly and doubly ionize nitrogen, and doubly and triply ionize oxygen.As the source of photoionization we consider either a black body with varying temperature T * (mimicking the impact of individual massive stars), an AGN, or a single stellar population.For the AGN case, we adopt the default template in Cloudy (based on Zamorani et al. 1981;Francis et al. 1993;Elvis et al. 1994): with the IR cutoff determined by k b T IR = 0.136 eV or T IR =1580 K; the Blue Bump temperature set to T BB = 1.5 × 10 5 K; and a = 2.4 × 10 6 is a coefficient set to define the relative strength of the Blue Bump and the X-ray corona, which in our case is constrained by the requirement that α ox = −1.4.For the single stellar population model, we refer to Bruzual & Charlot (2003) models with solar metallicity and for a Salpeter initial stellar mass function, computed at various time steps after the initial burst.The hardness of the radiation field from these templates, gauged by Q(O + )/Q(N), increases with T * for the black body templates, reaching the AGN value for T * ≈ 50000 K.The single stellar population shows a non-monotonic behavior, due to the rapid evolution of the most massive stars from the main sequence to the asymptotic giant branch phase; but we find an overall trend, with the Q(O + )/Q(N) ratio slowly decreasing in the first ∼ 10 Myr since the burst, followed by a rapid decline as all massive (O, B-type) stars evolve out of the main sequence.
The analytical model used for the ionized medium cannot be trivially expanded to the PDR/XDR regime.In this regime, collision excitation/de-excitation terms need to account for multiple collision partners (electrons, neutral H atoms, para-H 2 , ortho-H 2 , He, etc), the abundance of which depends on the ionization and dissociation conditions in different layers of the clouds.Furthermore, the chemistry itself of the neutral and molecular phases is more complicated (e.g., photo-dissociation of carbon monoxide molecules alters the abundance of neutral carbon atoms in the clouds).Opacity is often non-negligible (both for lines and continuum).Finally, many parameters are interconnected (e.g., the kinetic temperature of the gas, T gas , might be partially coupled with the dust temperature, T d ).For instance, Harrington et al. (2021) build on the formalism outlined in Weiß et al. (2007) to solve the radiative transfer problem simultaneously for the dust and the molecular gas.Their modeling involves seven free parameters (gas density, n; gas kinetic temperature, T gas ; dust temperature, T dust ; turbulence velocity; virial velocity gradient; carbon abundance, [C/H 2 ]; scale size of the emitting region) for each gas component.Other parameters, such as the [CO/H 2 ] abundance of the gas-to-dust mass ratio, are set independently.Their approach allows them to characterize various physical properties of the dust and molecular gas in a sample of Planck-selected sub-mm galaxies.While the method by Harrington et al. (2021) can in principle be applied here, our ob-servational constraints (in particular concerning the CO spectral line energy distribution) are too sparse to lead to informative results.For the analysis of the neutral and molecular regimes, we thus opt to stick to the Cloudy models presented in Pensabene et al. (2021) and summarized in the next section.

Cloudy modeling
We based our Cloudy analysis on the models presented in Pensabene et al. ( 2021) and Meyer et al. (2022).Namely, we model the ISM as a homogeneous plane-parallel slab of gas exposed to a radiation field.We independently model the ionized vs. neutral and molecular phases.
For all the input templates of photoionization sources, we fix the the ionization parameter U, which represents the density of ionizing photons per hydrogen atom: where Q(H) is the number of hydrogen-ionizing photons emitted per unit time by the photoionizing source, and r is the distance between the photoionizing source and the cloud.Rigopoulou et al. (2018) use the dust color C(88/122) as a proxy for the ionization parameter.Following their approach, the measured C(88/122) ≈ 1.4 corresponds to log U = [−1.9,−2.5] for a gas density n=[10, 1000] cm −3 , respectively (see their Fig. 4).As we aim to compare the emission of lines with different ionization energies close to or above the ionization energy of hydrogen, for this part of the analysis it is more critical to study the hardness of the radiation field at energies higher than the ionization energy of hydrogen.We thus use a fixed log U=−2 in our analysis.
We assume that the gas volume density is constant within the cloud, and we sample the line emissivity over a range of log n [cm while we use the default ISM abundances in Cloudy and a linear scaling with [O/H] for all the other elements.
We also include a CMB background computed at the source's redshift, and a turbulence with an RMS of 100 km s −1 .We stop the integration when we reach a ionized fraction of x HI = 0.05.
For the ionized medium, we run Cloudy for all the templates shown in Fig. 5.For the neutral and molecular phases, we adopt the models of Photon-Dominated Regions (PDRs) and X-ray Dominated Regions (XDRs) presented in Pensabene et al. (2021).In brief, we compute the radiative transfer using Cloudy for a grid of 15 × 18 × 3 combinations of gas volume density, intensity of the radiation field, and column density.For the former, we adopt log n [cm −3 ] = [2, 6] (0.29 dex spacing).For the latter, we adopt a range of log G 0 = [1, 6] (0.29 dex spacing) for PDRs, and log F X [erg s −1 cm −2 ] = [-2.0,2.0] (0.24 dex spacing) for XDRs.Contrary to what we have done for the ionized medium, we here fix the shape of the stellar template using a black body with T * =50,000 K, as the intensity of the radiation field is more critical in setting the PDR conditions.For the AGN, we adopt the same template as in eq.14.We adopted the default Cloudy recipes for the cosmic ray ionization rate background.
The total cloud column density lies in the range log N H [cm −2 ]= [22,24] (1 dex spacing).A higher column density N H > 5 × 10 23 cm −2 is required to properly model H 2 O and OH emission (e.g., Goicoechea et al. 2005Goicoechea et al. , 2006;;González-Alfonso et al. 2014;Pensabene et al. 2022).However, because of the sparsity of the water and hydroxyl transitions studied here, an extensive analysis of the radiative properties of these molecules is beyond the scope of this work.

Results
In this section, we use the observed line luminosities to infer physical properties of the ISM.

Opacity
If the line-emitting clouds are interspersed within the dusty medium, line emission should be corrected for opacity.In sec.3.2, we inferred a non-negligible opacity, τ ν ≈ 0.47 at the frequency of [C ii].For this opacity, the flux density emerging at ν 0 ∼ 1900 GHz is 2/3 of the intrinsic value.The higher the frequency, the larger the correction (see eq. 5), with the largest value amounting to 4.6 for [O iii].The extinction values computed at the frequencies of the lines studied in this work are listed in Tab. 2.
If we correct for opacity by scaling all the line luminosities by a factor e τ ν , we find that the intrinsic [N ii] 122 µm /[N ii] 205 µm ratio would be 1.77× higher than the observed value, while the intrinsic [O iii] 88 µm /[N ii] 122 µm would be 2.05× higher than the observed value.Our measurements for the molecular phase would be virtually unaltered, given that they mostly rely on lowerfrequency lines.As we will show in the following subsections, applying these corrections would result in a higher estimate of the electron density (n > 200 cm −3 ) and a slightly harder radiation field.However, the actual values of such corrections strongly depend on how the line-emitting gas clouds and the dust are distributed along the line of sight, which is unknown.Because of this uncertainty, and due to the modest impact that this correction would have on our final results, we opt not to apply the opacity correction to our fiducial measurements.

Inferred [N ii] 205µm and the origin of [C ii]
We estimate the luminosity of the [N ii] 205 µm line following the empirical trend presented in Lu et al. (2015) for a sample of local luminous IR galaxies: where the [N ii] 205 µm to [C ii] 158 µm line ratio is tied to the color of the dust emission, C(60/100), computed as the ratio between the rest-frame flux density at 60 and 100 µm.In the case of PJ183+05, C(60/100) = 0.84, yielding log L [NII] [L ]≈ 8.51 ± 0.23, where the uncertainty is dominated by the scatter in the relation from Lu et al. (2015).The [C ii] 158 µm and [N ii] 205 µm lines both have relatively low critical density (see Tab. 4, assuming electrons as collision partners) and have similar ionization energies (14.5 eV for nitrogen, 11.2 eV for carbon).In the ionized medium, their luminosity ratio is thus determined almost exclusively by the abundance ratio, [C + /N + ].This is confirmed by both our Cloudy modeling and our analytical prescription, as shown in Fig. A.1.We find that, for the ionized component alone, the expected [C ii] ion 158 µm /[N ii] 205 µm luminosity ratio is 6-12 at low densities, and 10-20 at high densities, with the spread completely dominated by the relative abundance.Romano (2023) has recently reviewed the evolution of carbon, nitrogen and oxygen in galaxies.Here we refer to the Nicholls et al. (2017) calibration, which yields [C + /H]=2.6 × 10 −4 and [N + /H]=6.2× 10 −5 at solar metallicities, or [C + /N + ]=4.28.On the other hand, in the literature it is common to follow Oberst et al. (2006), who rely on the abundance estimates in Savage & Sembach (1996), [C + /H]=1.4× 10 −4 and [N + /H]=7.9× 10 −5 .These slightly lower (higher) carbon (nitrogen) abundances result in a ∼ 2.
for a gas density of n ∼ 180 cm −3 and solar metallicities.We point out that this value is highly sensitive to the assumed relative abundances.If we adopt the traditional values from Savage & Sembach (1996), we find a much higher f ([CII] PDR ) = 0.81.

Electron density and size of the Hii regions
The ratio between the two far-IR fine-structure lines of [N ii] at 122 µm and 205 µm is independent of metallicity (as they are two transitions associated with the same ion).It is also insensitive to the hardness of the radiation field and to the gas temperature, due to the low energy of the 3 P 1 and 3 P 2 levels (E/k b = 70 and 188 K, respectively) compared to the typical electron temperatures in the ionized medium (T gas = 5000-20000 K).On the other hand, the two transitions have different critical densities (see Tab. 4).Thus, their luminosity ratio is sensitive to the electron density in the range 3-3000 cm −3 .This is shown in Fig. 6, where again we show both the results from our Cloudy models and from our analytical prescriptions.The two approaches lead to consistent results throughout the range of interest.Using the measured [N ii] 122 µm luminosity from Tab. 2, and the inferred [N ii] 205 µm from the previous section, we find a ratio of log[N ii] 122 µm /[N ii] 205 µm =0.56 +0.23  −0.16 , corresponding to a gas density of log n [cm −3 ] = 2.26 ± 0.35.
Our density estimate allows us to compute a fiducial volume occupied by the ionized gas, by comparing the observed line luminosity with the population fraction and gas density estimated via models.For T gas = 10, 000 K, n=180 cm −3 , and Z=Z , using These models are color-coded by the hardness of the impinging radiation field, parameterized as the effective temperature of black bodies T * .We also show the predicted ratios for our analytical prescriptions as solid lines, color-coded by gas temperature T gas .The observed ratio is marked with green shading.The predicted ratio is practically independent of T * or T gas , and is solely determined by the gas density n.
eq. 11 in the case of the [N ii] 122 µm transition, we obtain an effective volume of the ionized gas, V eff HII ∼ 0.2 kpc 3 .Assuming this volume is organized in a collection of N HII spherical Hii regions, the Stromgren radius would be ∼ 360 (N HII ) −1/3 pc.

Hardness of the radiation field and metallicity
The first ionization energies of carbon, nitrogen and oxygen are 11.26, 14.53 and 13.61 eV, respectively.Their second ionization happens at 24.38 eV, 29.60 eV, and 35.11 eV, respectively.Finally, the third ionization energy of oxygen is 54.93 eV.Thus, the ratios [O iii]/[C ii] ion and [O iii]/[N ii] are sensitive to the hardness of the photoionizing radiation, i.e., the relative number of photons produced per unit time in the band 35-55 eV vs. the band 14-25 eV (see Fig. 5).In addition, these line ratios also depend on the relative abundance of the elements, and on the gas density (see Tab. 4).Assuming the best-fit value for n from sec.4.3, we can use the Cloudy models and analytical prescriptions described in sec.Our predictions are shown in Figs.7 and 8.We find that the [O iii] 88 µm /[N ii] 122 µm ratio shows only a mild dependence on the gas metallicity.The evolution of the [N/O] relative abundance set by eq. 17 suggests that the ratio is constant at low Fig. 7.The [O iii] 88 µm /[N ii] 122 µm luminosity ratio as a function of metallicity Z and hardness of the photoionizing radiation, parameterized as the temperature of an equivalent black body, T * , computed from Cloudy models assuming a gas density of n = 180 cm −3 .The shaded area marks the observed line ratio and its 1-σ confidence interval.The line ratio is nearly independent of metallicity at Z ∼ < 0.5 Z (see eq. 17), while it is strongly sensitive to the hardness of the ionizing radiation.Our observations point to a value of T * =21,000-32,000 K. metallicities (Z<0.1 Z , in the regime of primary abundances, where the enrichment is dominated by core-collapse supernovae) and becomes nearly linear with Z only at Z>Z (where the secondary abundances arise due to delayed nucleosynthesis in intermediate-mass stars; see Nicholls et al. 2017).On the other hand, the [O iii] 88 µm /[N ii] 122 µm luminosity ratio is strongly sensitive to the hardness of the photoionizing radiation.In the case the input source is a black body, we find that that the expected line ratio changes by 2-3 dex for T * ranging between 20,000 and 40,000 K.The AGN template yields a very high [O iii] 88 µm /[N ii] 122 µm luminosity ratio, ∼ 50 times higher than observed.Finally, the single stellar population templates predict very high [O iii] 88 µm /[N ii] 122 µm luminosity ratios in the first ∼ 5 Myr, followed by a rapid decrease, with virtually no significant [O iii] 88 µm expected ∼ 10 Myr after the burst.The observed line ratio points to a T * ≈ 25, 000 K for Z=Z , or a burst age of ≈ 9 Myr.
We notice that our Cloudy models and analytical prescriptions agree on qualitative but not quantitative terms (see appendix A).Because the two approaches are consistent in the other diagnostics discussed so far (involving n, T gas , and Z), we argue that the discrepancies arise due to the simplistic assumptions made in the analytical prescriptions concerning the budget of ionizing photons Q(N) and Q(O + ).That is, our prescription assumes that all of the photons contribute to the photoionization of nitrogen and oxygen.A more realistic description would consider the competing role of hydrogen photoionization within different layers of the cloud, and the frequency dependence of the cross-section of the process (∝ ν −3 ).A more refined analytical description of the internal ionization structure of Hii regions Cloudy, assuming n = 180 cm −3 and Z=Z , for the various templates adopted in our study: A black body radiation of various temperatures, T * ; an AGN as described in eq.14; and a single stellar population with different burst ages, based on Bruzual & Charlot (2003).The observed line ratio is marked with a light green shading.The observed line ratio is consistent with the expectations for a black body of T * ≈ 25, 000 K or with a single stellar population with age ≈ 9 Myr.The AGN template, as well as templates involving either hotter black body temperatures, or lower stellar population ages, appear to over-predict the [O iii]/[N ii] luminosity ratio.
is possible (see, e.g., Yang & Lidz 2020).However, due to the steep dependence of the observed line ratio on T * , our simple model for the Hii regions results in a T * constraint that is only ∼ 1.7 times higher than the one derived from the much more refined Cloudy simulations -a discrepancy that has very little impact on the conclusions of our work.Thus, we argue that a more sophisticated model is unnecessary for the present work.
The observed [O iii]/[N ii] luminosity ratio points to B-type stars as the main driver of the photoionization budget.These stars have a stellar mass of ∼ 15 M , a luminosity of ∼ 10, 000 L , and a main-sequence lifetime of ∼ 10 Myr.A quasar contribution to the photoionization budget, while likely present, is not required to explain the observed gas properties.
The non-detections of the very high-J CO transitions in our sample only allow us to exclude the most extreme scenarios of pure XDR with strong incident radiation field, F X ∼ > 100 erg s −1 cm −2 and n > 10 6 cm −3 .
The combination of [C i] 370µm and [C ii] 158µm provides a powerful diagnostic of the neutral and molecular medium.Models of XDRs predict that the harder X-ray radiation should be able to penetrate deeper into the cloud, and enhance the heating and the photo-dissociation of carbon monoxide compared to the predictions for models where the source of radiation is young stars.This should manifest as a low [C ii]/[C i] ratio: < 10 for a column density N H = 10 23 cm −2 , over a wide range of cloud density (see, e.g., Pensabene et al. 2021).Decarli et al. (2022) examined the ratio in a sample of z ∼ 6 quasars, and found [C ii]/[C i] ratios in the range 20-100, strongly pointing toward PDR rather than XDR regimes being the norm in the host galaxies of high-z quasars.Here we measure a ratio of log [C ii]/[C i] = 1.55 +0.17 −0.13 .Even if we consider only the component of [C ii] that is not associated with the ionized medium, [C ii] PDR , the ratio drops to ∼ 15, still significantly higher than the predictions for XDRs.
In summary, all these lines of evidence disfavor a scenario in which the neutral and molecular medium in PJ183+05 is dominated by XDR conditions, but leave room for a combination of PDR and XDR, or a PDR-dominated regime.A more extensive census of both lower-and higher-J CO transitions is required to further address the relative importance of PDR and XDR conditions in the host of PJ183+05.
level introduces "Λ"-doubling splitting of each energy level.In addition, "hyperfine" splitting occurs due to nuclear and electron spin coupling.The 163 µm transition thus actually splits into a total of 6 transitions, at 1834.7350 GHz (N=2 + -1 − , F=1-1), 1834.7469GHz (N=2 + -1 − , F=2-1), 1834.7499GHz (N=2 + -1 − , F=1-0), 1837.7461GHz (N=2 − -1 + , F=1-1), 1837.8163GHz (N=2 − -1 + , F=2-1), and 1837.8365GHz (N=2 − -1 + , F=1-0).The hyper-fine structure splitting is negligible for typical extragalactic observations ( ∼ < 0.1 GHz at rest frame, or ∼ < 15 km s −1 , i.e., much smaller than the width of the lines).Thus, the Λdoubling implies that the 163 µm transition appears as a doublet, spaced by ≈ 500 km s −1 .Fig. 10, right compares the line luminosity of the H 2 O and OH transitions in PJ183+05, normalized to the total IR luminosity, to the values measured in other z > 6 quasars (data from Pensabene et al. 2021Pensabene et al. , 2022) ) as a function of the upper energy level of the transitions.Strikingly, all of the detected line transitions point to a high H 2 O-or OH-to-IR luminosity ratio, about 3 times higher than what is observed in other quasars at similar redshifts.The modest number of transitions observed in our study makes it impossible to draw any robust conclusions on the origin of this discrepancy.We speculate on two scenarios: 1) The IR luminosity estimate is correct, in which case the excess may be associated with a low contribution of radiative excitation (as reflected by the non-detection of high-J H 2 O lines) in PJ183+05 compared to other quasars.Because of the high energies of the involved levels, this scenario would favor the presence of shocks (e.g., van Dishoeck et al. 2011).Alternatively, 2) the IR luminosity estimate in PJ183+05 is underestimated.The excellent sampling of the dust continuum shown in Fig. 2 pins down the rest-frame far-IR side of the dust emission pretty accurately, but leaves room for an additional hot dust component at mid-IR frequencies (such as an AGN torus), where various H 2 O and OH transitions relevant for the IR pumping mechanism reside.

Line profiles and relative velocities
Because of their relatively high energies, water levels at J up > 2 and hydroxyl levels above the ground state are expected to be significantly populated via collisions only in the presence of shocks.Furthermore, P-Cygni profiles have been observed in various water and hydroxyl transitions in lower-redshift, IRluminous galaxies, thus demonstrating that these species trace dense molecular outflows (e.g., Goicoechea et al. 2005Goicoechea et al. , 2006)).Similarly, high-ionization lines such as optical [O iii] commonly show blue-shifted wings and broadened spectral profiles that reveal the presence of outflows of the ionized medium (e.g., Nesvadba et al. 2008;Carniani et al. 2015;Bischetti et al. 2017).
In Fig. 11 we compare the spectral profiles of all the transitions studied in our work.We include in the comparison the three J up ≤ 4 water transitions that remain undetected in our study.We compare each transition to the [C ii] line, which has the highest S/N.The velocity profiles of all the detected transitions appear in good agreement at the S/N levels of our observations.All the detected lines show negligible velocity shifts along the line of sight compared to the rest frame of the galaxy as defined by the [C ii] line, with all the velocity shifts, ∆v = c δz/(1 + z), significantly smaller than the adopted channel width (90 km s −1 ).The line width is to first order identical in all the transitions.We do not find any clear evidence of P-Cygni profiles.We note the presence of a tentative absorption feature shifted by about -750 km s −1 in the H 2 O 3(3,1)-3(2,2), 3(0,3)-2(1,2), and possibly 4(3,2)-4(2,3) lines.Noticeably, Butler et al. (2022) found a similar absorption feature in the OH 119 µm transition for PJ183+05, suggesting that the H 2 O absorption found here is real.We note that this absorbing component does not match the velocity of the proximate damped-Lyα absorber reported in Bañados et al. (2019) (which shows a velocity offset of ≈ 1400 km s −1 ).

Masses
Having constrained the excitation conditions of the ionized species examined in our work (Fig. 7), we can use the information on the relative population of the energy levels to infer the total mass in each ion, and thus the total gas mass via our abundance assumptions.We show these gas estimates in Fig. 12, where we compare them with the mass estimates inferred from the dust in sec.3.2, scaled assuming δ g/d = 100 (Berta et al. 2016); from CO(7-6), assuming r 71 = 0.38 and α CO = 0.8 M (K km s −1 pc 2 ) −1 (Carilli & Walter 2013); and from [C i] and [C ii] PDR , following Weiß et al. (2003Weiß et al. ( , 2005) ) and Venemans et al. (2017b), where we assume T ex = T gas .We find that carbon fine-structure lines suggest the presence of ∼ 10 7 M in both neutral and singly-ionized carbon.On the other hand, from the observed luminosities of [N ii] and [O iii], we infer ∼ 2 × 10 6 M in each ion.Unsurprisingly, the mass estimates based on [C i] and [C ii] PDR are sensitive to the adopted kinetic energy of the collision partners at T gas ∼ < 100 K, i.e., comparable to the 3 P 1,2 and 3 P 3/2 energy levels, respectively.We scale these mass estimates by the abundances computed following eqs.16-17 for Z=Z (Nicholls et al. 2017).We find that the low-ionization species ([C ii], [N ii]) yield very consistent ionized gas mass estimates, M HII = (2 − 3) × 10 10 M .The mass estimate based on [O iii] is ∼ > 5× lower, likely due to the fact that only a small fraction of oxygen atoms are highly ionized.As for the molecular medium, the mass estimates based on singlyionized and neutral carbon are comparable.We sum them to infer the total mass associated with carbon atoms and ions.We find a corresponding molecular gas mass M H2,[CI]+[CII] ≈ 5×10 10 M at T gas ∼ > 100 K.This is in excellent agreement with the CO-based mass estimate, M H2,CO = 4.5 × 10 10 M , and a factor ∼ 2 lower than the dust-based estimate, M H2,dust = 8.7 × 10 10 M .While all the mass estimates rely on uncertain conversion factors (see discussions in, e.g., Dunne et al. 2021Dunne et al. , 2022;;Decarli et al. 2022), all the methods independently converge towards a massive gaseous reservoir in both molecular [(0.5−1)×10 11 M ] and ionized [(2 − 3) × 10 10 M ] phase.

Star formation rate estimates
Assuming that star formation is the main driver of the gas excitation and photoionization, and of the dust heating, we can convert some of the observed line and continuum luminosities into SFR estimates, using empirical scaling relations.In particular, following the scaling values reported in Tab.The discrepancies in the SFR estimates are much larger than formal measurement uncertainties.The scatter in the adopted SFR-luminosity relations is ∼ < 0.3 dex, insufficient to account for a ∼ 1 dex discrepancy between the [O iii]-and IR-based SFR estimates.Hence, the differences must have a more fundamental origin.One way to do so is by correcting for opacity as in Sec.4.1.Doing so, we find SFR [OIII] =500 M yr −1 and SFR [CII] =1190 M yr −1 .Another way to reconcile the different SFR estimates is to account for the different star-formation time scales probed by [C ii] and dust (∼ 100 Myr) and doubly-ionized oxygen (∼ 10 Myr; see Kennicutt & Evans 2012).Finally, other mechanisms might be at play.For example, if star formation is heavily dust-enshrouded, dust grains may absorb and reprocess the far-UV photons required to doubly-ionize oxygen atoms, thus suppressing [O iii] emission.In addition to these caveats, the uncertainties in the relative abundances, metallicity, and excitation properties discussed in the previous subsections should serve as a warning on the applicability of SFR prescriptions, especially when only one or two independent tracers are available.

Conclusions
We presented a multi-line study at FIR wavelengths of the luminous quasar PJ183+05 at z = 6.4386.We reported the detection of , OH 163µm , and two H 2 O transitions, and we place limits on the line emission of various high-J CO lines and other water transitions.At the same time, we sampled the dust continuum over a large fraction of the FIR band.We compared our measurements with Cloudy models and analytical prescriptions, with the goal to infer an insight of the physics of the interstellar medium in the host galaxy of a luminous quasar at cosmic dawn.We find that: i-The dust emission is well described by a modified black body with T dust = 47.0 +1.5 −2.0 K, an emissivity index β = 1.84 +0.15 −0.16 , a dust mass of log M dust /M =8.94 +0.06  −0.05 , and a size of the emitting region of Ω s = 0.155 +0.029 −0.022 arcsec 2 .This implies straints on the gas metallicity, Z, but tight constraints on the hardness of the photoionization.Cloudy models suggests that the photoionization source has a spectrum comparable to a black body with T * ≈ 25, 000 K or a single stellar population observed ≈ 9 Myr after the burst (assuming for Z=Z ).This implies that B-type stars with a luminosity of ∼ 10, 000 L and a main-sequence lifetime of ∼ 10 Myr are the main drivers of the photoionization.v-The observed [C ii]/[C i] line ratio, and the non-detection of high-J CO transitions, suggest that the neutral/molecular medium is not dominated by X-ray excitation, although we cannot exclude the presence of an X-ray Dominated Region.vi-Water and hydroxyl line emission appears brighter than what one would expect based on the IR luminosity.This is indicative of the presence of shocks that collisionally excite the higher energy levels of the H 2 O and OH molecules, or of a warm dust component that is not apparent from the dust continuum sampling in the far-infrared range presented here.vii-The line profiles of all the line transitions studied in this work are consistent at the level of S/N of our observations, irrespective of the gas phase they arise from (ionized, neutral, or molecular).A tentative absorption feature at ∆v ∼ −750 km s −1 is seen in the spectra of two water transitions.Further support to the robustness of this feature comes from a similar absorption reported in the OH 119µm line from Butler et al. (2022).This study demonstrates the power of multi-line observations in understanding the dominant radiative processes and astrophysical mechanisms at play in galaxies at cosmic dawn.These studies are possible in particular thanks to ALMA's exceptional sensitivity.For instance, other observations of molecules and atoms in PJ183+05 could further expand our understanding of the physical conditions in the star-forming ISM in this quasar.Expanding multi-line investigations to other sources will enable to test whether other quasars show the properties found in PJ183+05.The James Webb Space Telescope is poised to further expand this field of research, by complementing the study of the rest-frame FIR wavelengths with key diagnostics based on restframe optical tracers of the interstellar medium.

Fig. 1 .
Fig. 1.Observed ALMA spectra of PJ183+05 (histograms).The fitted continuum+line models are shown as thick red lines.Main transitions are marked with vertical dotted lines.

Fig. 2 .
Fig. 2. Observed dust spectral energy distribution in PJ183+05.We plot the observed spectra and best-fit continuum flux density estimates from the spectral analysis as black lines and red points with error bars, respectively.The best fit model is shown in blue, while a random subset of the models within the 1-σ confidence level are shown in gray for reference.The broad range of frequencies sampled in our study allows us to precisely pin down the dust temperature in this source.

Fig. 3 .
Fig.3.Corner-plot posterior distributions for the parameters used in the dust SED modeling: the dust temperature, T dust ; the observed solid angle of the emitting region, Ω s ; the mass of dust, M dust ; and the dust emissivity index, β.In the corner plots, darker points show values corresponding to higher posterior probabilities.We also show the marginalized posterior distributions for the dust optical depth at the frequency of [C ii], τ 1900 GHz ; and of the infrared luminosity, L IR .In all the marginalized distributions, the median values are marked with a dashed line, and the values corresponding to the 1-σ confidence intervals are marked with dotted vertical lines.

Fig. 4 .
Fig. 4. Results from the analytical model of the fine-structure line emission associated with O ++ , N + and C + ions in the ionized medium, as described in Sec.3.3.2.Energy levels are populated by collisions with electrons, and depopulated by collisions and by radiative de-excitation.We assume local thermal equilibrium, thus the number of collisions depends only on the gas temperature T gas and density n.The vast majority of electrons reside in the 3 P J levels.

Fig. 5 .
Fig. 5. Template comparison for the photoionizing source.Top panels: The rest-frame flux density of the AGN template (solid black lines), of black bodies with different temperatures T * (dotted lines in the lefthand panel) and single stellar populations with different ages (dotted lines in the right-hand panel).The frequency threshold of photons responsible for the ionization of hydrogen, nitrogen, oxygen relevant to this study are marked for reference.Bottom panels: The ratio of O + to N-ionizing photons produced by different templates, as a function of T * (black bodies; left-hand panel) and burst age (right-hand panel).Empty circles refer to values integrated at any ν 0 > ν ion , while filled circles refer to the 'net' values (see eq. 13).

Fig. 6 .
Fig.6.Predicted [N ii] 122 µm /[N ii] 205 µm luminosity ratio, as a function of the gas density.The results from our Cloudy models are shown as colored symbols, with small offsets introduced for the sake of clarity.These models are color-coded by the hardness of the impinging radiation field, parameterized as the effective temperature of black bodies T * .We also show the predicted ratios for our analytical prescriptions as solid lines, color-coded by gas temperature T gas .The observed ratio is marked with green shading.The predicted ratio is practically independent of T * or T gas , and is solely determined by the gas density n.
3.3.2 to predict the [O iii] 88 µm /[N ii] 122 µm ratio as a function of the hardness of the radiation field, parameterized as the effective temperature of a black body, T * , and of metallicity Z.We opt to focus on the [O iii]/[N ii] ratio rather than on the [O iii]/[C ii] ion ratio because of the uncertainties on f ([C ii] PDR ) discussed in sec.4.2.

Fig. 8 .
Fig. 8.The [O iii] 88 µm /[N ii] 122 µm luminosity ratio, computed withCloudy, assuming n = 180 cm −3 and Z=Z , for the various templates adopted in our study: A black body radiation of various temperatures, T * ; an AGN as described in eq.14; and a single stellar population with different burst ages, based onBruzual & Charlot (2003).The observed line ratio is marked with a light green shading.The observed line ratio is consistent with the expectations for a black body of T * ≈ 25, 000 K or with a single stellar population with age ≈ 9 Myr.The AGN template, as well as templates involving either hotter black body temperatures, or lower stellar population ages, appear to over-predict the [O iii]/[N ii] luminosity ratio.

Fig. 9 .
Fig. 9. Observed constraints on the CO spectral line energy distribution in PJ183+05.Our measurements are shown as black filled circles.Top and bottom panels show the models for PDRs and XDRs, respectively.We mark models, taken from Pensabene et al. (2021), with increasingly darker colors at increasing intensity of the radiation field.Dotted, shortdashed and long-dashed lines refer to models with gas density log n [cm −3 ] = 2, 4, 6 respectively.All models are normalized to the observed CO(7-6) line flux.We plot the observed line fluxes or upper limits as black circles.In addition, we show the predicted CO(1-0) flux, as derived from the dust continuum, as described in sec.4.5.Our constraints exclude a strong XDR environment, but leave room for a combination of PDR and XDR environments in the excitation conditions of the molecular ISM in PJ183+05.

Fig. 10 .
Fig. 10.Left: Scheme of the J≤4 energy levels of the H 2 O molecule, and for the Ω=1/2 and 3/2 levels of the OH molecule.Solid lines mark the transitions detected in our study, dashed lines show the transitions that we observed but not detected.The channels of IR pumping relevant for this study are also marked with dotted arrows.Right: Comparison between the luminosity of H 2 O and OH transitions, normalized to the dust IR luminosity, in PJ183+05 and in other z > 6 quasars from Pensabene et al. (2021, 2022).All of the water and hydroxyl transitions detected in our work show line-to-IR luminosities that are ∼ 3× higher than the values observed in other quasars at the same redshift.

Fig. 11 .
Fig. 11.Comparison between the line profiles of the main transitions in our study (shaded histograms) and the [C ii] line profile (black histograms).The grey shading masks regions that are not covered within our observations or are contaminated by other emission lines.To first order, all the lines appear to have a consistent bulk velocity and width.A tentative H 2 O absorption is detected at -750 km s −1 in the 3(3,1)-3(2,2) and 3(0,3)-2(1,2) transitions, suggesting the presence of an outflow.
viii-The mass estimates in ionized and molecular gas, computed based on different tracers, point to a massive gaseous reservoir in both phases, with (2 − 3) × 10 10 M in ionized gas, and (0.5 − 1.0) × 10 11 M in the molecular phase.ix-The host galaxy of PJ183+05 is forming stars at high rates, SFR IR = 1330 M yr −1 .Estimates based on either [C ii] and [O iii] point to lower values, SFR [CII] = 750 M yr −1 and SFR [OIII] = 110 M yr −1 .The discrepancy is partially mitigated if we account for the dust opacity suppressing the observed line luminosities.In addition, different time scales probed by [O iii], [C ii] and dust, as well as other processes (e.g., the '[C ii] deficit') might be responsible for the SFR discrepancies.

Table 1 .
Summary of the observations.The RMS refers to the observed noise per channel in the cleaned data cube, for a channel width of 90 km s −1 , except for [C ii] (setup 4) which is computed over channels of 30 km s −1 (see text).
and velocity shifts have been reported among [O iii] and [C ii]