A benchmark for extreme conditions of the multiphase interstellar medium in the most luminous hot dust-obscured galaxy at z = 4 . 6

WISE J224607.6–052634.9 (W2246–0526) is a hot dust-obscured galaxy at z = 4 . 601, and the most luminous obscured quasar known to date. W2246–0526 harbors a heavily obscured supermassive black hole that is most likely accreting above the Eddington limit. We present observations with the Atacama Large Millimeter / submillimeter Array (ALMA) in seven bands, including band 10, of the brightest far-infrared (FIR) ﬁne-structure emission lines of this galaxy: [OI] 63 µ m , [OIII] 88 µ m , [NII] 122 µ m , [OI] 145 µ m , [CII] 158 µ m , [NII] 205 µ m , [CI] 370 µ m , and [CI] 609 µ m . A comparison of the data to a large grid of C loudy radiative transfer models reveals that a high hydrogen density ( n H ∼ 3 × 10 3 cm − 3 ) and extinction ( A V ∼ 300mag), together with extreme ionization (log( U ) = − 0 . 5) and a high X-ray to UV ratio ( α ox ≥ − 0 . 8) are required to reproduce the observed nuclear line ratios. The values of α ox and U are among the largest found in the literature and imply the existence of an X-ray-dominated region (XDR). In fact, this component explains the a priori very surprising non-detection of the [OIII] 88 µ m emission line, which is actually suppressed, instead of boosted, in XDR environments. Interestingly, the best-ﬁtted model implies higher X-ray emission and lower CO content than what is detected observationally, suggesting the presence of a molecular gas component that should be further obscuring the X-ray emission over larger spatial scales than the central region that is being modeled. These results highlight the need for multiline infrared observations to characterize the multiphase gas in high redshift quasars and, in particular, W2246–0526 serves as an extreme benchmark for comparisons of interstellar medium conditions with other quasar populations at cosmic noon and beyond.


Introduction
"Cosmic noon" (z ∼ 1 − 3) and the epoch leading to it (z ∼ 3 − 6) are the two most important periods of cosmic history in terms of galaxy growth and supermassive black hole (SMBH) activity.It is during these periods that SMBHs and their hosts assembled most of their mass, as attested by the tight relationship between the stellar mass of galaxies and their central SMBH mass across cosmic time (e.g., Magorrian et al. 1998;Kormendy & Ho 2013).Understanding the properties of the interstellar medium (ISM) of actively star-forming galaxies with an actively accreting SMBH provides important insights into the co-evolution of a SMBH and their host galaxies.Additionally, observations of such extreme systems allow meaningful constraints to be placed on cosmo-logical simulations that aim to reproduce correlations between galaxy properties on spatial scales ranging from approximately parsecs (SMBH gas accretion) to tens of kiloparsecs (evolution of the host galaxy).
In particular, quasars, powered by accreting SMBHs at their centers, exert a profound influence on the ISM that surrounds them.The intense radiation and high-energy photons emitted by active galactic nuclei (AGNs) can significantly impact the physical and chemical conditions of the ISM, playing a crucial role in shaping the evolution of galaxies and regulating the growth of their central SMBHs (e.g., Di Matteo et al. 2005).Nevertheless, most of the accretion onto SMBHs is expected to be heavily obscured by dust and gas (Hickox & Alexander 2018), making the study of heavily obscured quasars key to understand the interaction between AGN and their surrounding gas.
Among the large variety of dust enshrouded quasars, hot, dust-obscured galaxies (Hot DOGs;Eisenhardt et al. 2012;Wu et al. 2012) are a population of hyperluminous quasars (L bol ≳ 10 13 L ⊙ ), firstly discovered by NASA's Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010).Hot DOGs are characterized by very high bolometric luminosities (L bol ≳ 10 13 L ⊙ ) and high dust temperatures (> 60 K, Wu et al. 2012), powered by highly obscured AGNs that dominate the spectral energy distributions (SEDs) beyond λ > 1 µm (Tsai et al. 2015), and Eddington-limited accretion (Wu et al. 2018).Signs of quasar feedback have been observed in the Hot DOG population, including the highly turbulent [CII] 158µm gas (Díaz-Santos et al. 2021) and the high-velocity ionized gas outflows detected in the rest-frame optical (Finnerty et al. 2020;Jun et al. 2020).Individual and statistical studies also suggest that Hot DOGs live in over-dense environments (Jones et al. 2014;Assef et al. 2015;Ginolfi et al. 2022;Luo et al. 2022;Zewdie et al. 2023) and their activity is likely merger-driven (Fan et al. 2016;Díaz-Santos et al. 2018).Moreover, these objects could be experiencing several recurrent episodes of extreme accretion and obscuration (Díaz-Santos et al. 2021).This scenario is in agreement with recent cosmological simulations of high redshift quasars (Anglés-Alcázar et al. 2021), but it is in tension with the paradigm for the formation of luminous quasars in the local Universe, that is, involving a single encounter of two massive, gas-rich spiral galaxies (e.g., Sanders et al. 1988;Hopkins et al. 2008;Treister et al. 2012).WISE J224607.6-052634.9(hereafter W2246-0526) is a Hot DOG at z [CII] = 4.601 (Díaz-Santos et al. 2016) with a bolometric luminosity of 3.6 × 10 14 L ⊙ , making it one of the most luminous galaxies known to date (Tsai et al. 2015(Tsai et al. , 2018)).It is part of a multiple merger (Díaz-Santos et al. 2018), where W2246-0526 is connected to at least three galaxy companions by dusty tidal streamers.W2246-0526 hosts an AGN with a SMBH mass of M • = 4.0 × 10 9 M ⊙ , and it radiates well above the Eddington limit, with an Eddington ratio of λ Edd = 2.8 (Tsai et al. 2018, assuming that the MgII 2799 Å emission line is being broadened by the SMBH gravity).W2246-0526 is a radio-quiet quasar (according to the classification of Kellermann et al. 1989), but the central SMBH powers parsec-scale radio activity (Fan et al. 2020).Given the extreme conditions, W2246-0526 provides a key laboratory to study the effect of the AGN on its surrounding ISM.Moreover, because of the extreme luminosity of W2246-0526, one might expect bright far-infrared (FIR) spectral lines that trace the multiphase ISM, and that are straightforward to detect with the Atacama Large Millimeter/submillimeter Array (ALMA).
ALMA covers a wavelength range where, at z ∼ 4.6, finestructure lines from [OI] 63µm to [CI] 609µm can be observed.Pre-vious observations of [CII] 158µm and CO(J = 2−1) (Díaz-Santos et al. 2016, 2018) suggest that the ionized, neutral, and molecular phases of the ISM are very turbulent, likely due to the feedback from the central AGN, out to scales of at least a few kiloparsecs.Therefore, we can expect all the fine-structure lines from all regions to be strongly affected by the central energy source.
From the photodissociation regions (PDRs) paradigm, we expect many of these lines (e.g., [CII], [OI], and [CI]) to arise from the surface of molecular clouds that are exposed to UV radiation.In the ionized gas phase outside the PDR, we find ionized nitrogen, and doubly ionized oxygen, which originates in highly ionized HII regions requiring radiation from early-type O stars or an AGN.[NII] and [OIII] lines will therefore trace this region.Tracking from the UV-exposed surface of the PDR, carbon is initially ionized, and oxygen is in atomic form.This region is therefore a strong emitter of both [CII] and [OI].Beyond, the far-UV field is extenuated to the extent that carbon is in atomic form.Deeper into the cloud, carbon is found in CO.The neutral carbon region, traced in its [CI] line emission, is relatively much thinner and cooler than the ionized carbon region.However, some AGNs can also generate X-ray-dominated regions (XDRs, e.g., Maloney et al. 1996), which significantly changes both the ionization structure and the heating of the molecular cloud, and therefore the regions from where finestructure FIR lines arise (see Fig. 1).Recent studies targeting these and other FIR lines, in combination with ISM models, have proven to be very successful in characterizing the properties and phases of the ISM of high-redshift quasars (Novak et al. 2019;Pensabene et al. 2021;Meyer et al. 2022;Decarli et al. 2023).
In this paper, we analyze ALMA observations of the and [CI] 609µm emission lines in W2246-0526.Modeling the AGN and ISM with the brightest FIR lines observed with ALMA provides an understanding of the role of quasar feedback in galaxy evolution.
The paper is structured as follows: Section 2 describes the observational setup, data reduction, and line measurements of the ALMA observations.In Sect. 3 we analyze individual line ratios that provide diagnostics of the ISM, describe the setup of the ISM modeling, and fit all the lines with the models.In Sect. 4 we discuss the results and place them into a wider context.Finally, in Sect. 5 we summarize and draw our conclusions.
Throughout the paper, we assume a flat ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 and Ω M = 0.3.For reference, in this cosmology, 1 ′′ is equivalent to 6.54 kpc at z = 4.601.We reduced and calibrated the data using the Common Astronomy Software Applications (CASA; McMullin et al. 2007).We processed the data with the ScriptForPI.pyprovided through the ALMA science archive to generate the measure-  .We note that n cr,e were extracted from the Cloudy database, computed for a gas temperature of 10000 K.The last two columns are taken from Draine (2011), with the n cr,H computed for a gas temperature of 100 K. ment sets, using the version of the pipeline used by the observatory to reduce the data from each Cycle.We note that the automatic pipeline was run without the need for manual intervention in all bands.We imaged the measurement sets to create spectral cubes for the emission lines using the tclean task of CASA v6.4.3.We adopted a uniform pixel scale of 0.05 ′′ to sample the synthesized beams (ranging from ∼ 0.3 ′′ to 0.5 ′′ in the different observations).We run the Hogbom cleaning algorithm (Högbom 1974) to a flux density threshold of two times the root mean square of each cube with the Briggs weighting mode set to a robustness = 2.0 (close to natural weighting), and using a mask with radius 1.25 ′′ centered at the peak luminosity pixel, corresponding to ∼ 8.2 kpc at z = 4.601, enough to cover the host galaxy.We created the line cubes by combining appropriate spectral windows to select the continuum from neighboring line-free channels and subtracting it using the CASA task uvcontsub.The intensity maps of all the lines as well as each of the corresponding spectra are presented in Fig. 2.

Observations, data reduction, and emission line measurements
W2246-0526 is surrounded by a complex structure of dust (Díaz-Santos et al. 2018), close companion galaxies (Díaz-Santos et al. 2016), and Lyman-break galaxies at larger scales (Zewdie et al. 2023).We note in Fig. 2 the resolved compan-ion galaxy northwest of W2246-0526 in the [CII] 158µm emission map (labeled as C1 in Díaz-Santos et al. 2018).The remaining lines are barely resolved or not resolved at all.W2246-0526 itself extends over ∼ 1 ′′ (∼ 7 kpc).The smallest maximum recoverable scale in our observations is ∼ 2.7 ′′ for ALMA Bands 9 and 10, and therefore we can be sure that all the flux of the host galaxy is being recovered.To perform the photometry and extract the spectra we used an aperture of diameter 0.7 ′′ centered at the peak emission pixel, enough to encircle the Hot DOG but not any of the known companion galaxies.This aperture size ensures that 99% of the line flux is included assuming a Gaussian distribution the size of the largest beam in our observations (0.517 ′′ × 0.440 ′′ ), meaning that effectively no aperture correction is needed.The continuum-subtracted spectra were then fitted with a single Gaussian (enough to describe all the line profiles) with the central wavelength of the line, peak intensity, and full width at half maximum (FWHM) as free parameters.To obtain the moment 0 map of the lines, the emission line cubes were collapsed using the channels within ±2.75σ of the line peak (where σ is the standard deviation of the fitted Gaussian), a compromise between including most flux of the Gaussian to a level of ∼ 99.5%, and avoiding including noise beyond the line.
Fig. 1: Sketches of the emission of the targeted lines in the different phases of the ISM, for a PDR on the top and an XDR on the bottom.The incident photons, far-UV for PDRs and X-rays for XDRs, enter from the left.Thick dashed lines separate the different hydrogen phases: ionized (H + ), neutral (H), and molecular (H 2 ).We do not show a scale regarding the depth of the cloud because it strongly depends on the ionization source and cloud properties.We note that in the XDR case, ionized transitions arise also where the gas is mostly neutral, and neutral carbon emission is distributed much further into the cloud since X-rays can penetrate much more deeply than far-UV photons.
As shown in Table 3, all lines have a similar FWHM around 600-700 km s −1 , with the exception of [NII] 122µm , [NII] 205µm , and [CI] 370µm , with ∼ 800 km s −1 , ∼ 900 km s −1 , and ∼ 350 km s −1 respectively.The range of line FWHM is close to the value previously found for [CII] 158µm in W2246-0526 and other Hot DOGs (Díaz-Santos et al. 2021).Interestingly, it is broader than the ∼ 200 km s −1 typically observed for CO(1 − 0) in Hot DOGs (Penney et al. 2020).In general, the line FWHM in W2246-0526 are higher than high-redshift star-forming galaxies (∼ 250 km s −1 in Béthermin et al. 2020), optical quasars (∼ 385 km s −1 in Decarli et al. 2018), or infrared quasars (from ∼ 200 km s −1 to ∼ 500 km s −1 in Trakhtenbrot et al. 2017).We do not detect any evidence of underlying, broad emission line components (which may indicate outflows), suggesting that outflows may only be clearly identified in optical spectra of highly ionized gas (Finnerty et al. 2020), or that W2246-0526 does not have significant outflow activity at the moment.
The line luminosities are calculated following Solomon et al. (1992): where S line ∆υ is the velocity integrated flux of the line over a range of ±2.75σ from the line peak, in Jy km s −1 ; ν rest is the rest frequency in GHz; and D L is the luminosity distance in Mpc.We also calculate the luminosity using the analytic expression of the Gaussian fit integral, especially for the cases in which the whole line profile is not covered by the observations.The fluxes and luminosities for all lines are presented, among other properties, in Table 3.The [CI] 609µm and [OIII] 88µm lines are not detected and we use them as upper limits in our analyses, calculated as the integrated luminosity of a Gaussian assuming a FWHM of 600 km s −1 and a peak of 1σ detection.To make sure the non-detections are not a product of the aperture choice or channel binning, we extracted the spectra using different aperture sizes and channel averaging.The lines were not detected in any case.We further discuss the non-detection in Sect.4.4.

Results
The targeted lines are among the brightest FIR fine-structure lines that can be observed in star-forming and AGN-dominated galaxies, and they trace different gas phases and their physical properties.When modeling and studying the ISM phases, a common division is made in the literature between PDRs and XDRs, depending on the dominant source of radiation (optical-UV photons or X-rays, respectively; see Fig. 1).X-ray photons can penetrate deeper into the ISM than UV photons, leading to a change in the ionization structure, chemistry, and physical conditions within the cloud.The difference between PDRs and XDRs, in terms of their effect on the gas phases and FIR lines emission, has been extensively discussed in the literature.For further details, we refer the reader to the review of Wolfire et al. (2022).Our Cloudy modeling indicates that the observed line ratios in W2246-0526 are reproducible with a single XDR component, which dominates the observed emission.We begin this section by discussing individual line ratios.Next we describe the setup of the CLOUDY models, followed by an explanation of the insights the models provide on the line ratio diagnostics, and conclude with a description of the best-fit model.

Individual line ratio diagnostics
Before interpreting the observations using the Cloudy models, we analyze key individual line luminosity ratios that are very useful to infer specific properties of the ISM with limited information.These line ratios are presented in Table 4, with the Gaussian fits to the observed data, and the best-fit model explained in Sect.3.4.We also report the observed line ratios when the FWHM is fixed to 600 km s −1 for all the lines, to discard a major influence in the ratios due to line width differences.The same line ratios are shown in Fig. 3.

Ionized carbon and nitrogen
Owing to their ionization potentials, the emission from the [CII] 158µm line in a PDR arises from both the neutral and the ionized gas medium, while [NII] 205µm is produced only in the ionized phase.Because of this, the ratio between [CII] 158µm and [NII] 205µm has been extensively used to disentangle the fraction of [CII] 158µm produced in PDRs (e.g., Abel 2006;Oberst et al. 2006;Goldsmith et al. 2015;Croxall et al. 2017;Díaz-Santos et al. 2017).Within this context, the observed [CII] 158µm /[NII] 205µm ratio of 20.5 ± 4.3 would correspond to an ∼ 80% of the [CII] 158µm line emission produced in the neutral ISM (following the model of Oberst et al. (2006) for n e ≥ 10 3 cm −3 ).However, as mentioned before, the presence of an XDR changes the chemistry of the ISM, and as described by Pereira-Santaella et al. ( 2017), X-ray radiation from an AGN     145µm 5.07 ± 0.44 6.64 ± 1.33 3.47 ± 0.30 4.54 ± 0.91 677 ± 117 0.410 × 0.330 -83.6 11.6 1.47 [CII] 158µm 10.5 ± 0.4 11.1 ± 0.5 6.60 ± 0.26 6.98 ± 0.34 590 ± 23 0.385 × 0.355 48.5 26.3 0.58 [NII] 205µm 0.68 ± 0.08 0.70 ± 0.14 0.33 ± 0.04 0.34 ± 0.07 894 ± 140 0.517 × 0.440 -80.5 7.2  0.12  [CI] 370µm 0.41 ± 0.15 0.45 ± 0.22 0.11 ± 0.04 0.12 ± 0.06 368 ± 128 0.422 × 0.290 -49.8 3 Note: Columns correspond to the line, flux and luminosity of both the data and the Gaussian fit, full width at half maximum (FWHM) of the Gaussian fit, size, and angle of the beam, signal-to-noise ratio (S/N) of the line, and depth (root mean square; r.m.s.) per channel.The channel width for all but [OIII] 88µm and [CI] 609µm FIR lines is 50 km s −1 while the channel width for those two lines (denoted by * ) is 100 km s −1 , in order to enhance the S/N and possible detection of the lines.The upper limits correspond to the integrated luminosity of a Gaussian assuming a FWHM of 600 km s −1 and a peak of 1σ detection.Note: Columns correspond to the line ratios measured directly from the data, the Gaussian fit of the lines, the ratio of the best-fit model described in Sect.3.4, and also the best-fit model ratio but when assuming a Gaussian with a fixed FWHM of 600 km s −1 for all the lines.For non-detections we report 1σ upper limits.
can cause nitrogen to be single ionized in regions where hydrogen is mainly neutral.Both [CII] 158µm and [NII] 205µm lines will arise from similar regions of the ISM, with their ratio strongly depending on the parameters of the model.This ratio is shown in the x-axis of the bottom left panel of Fig. 3.

Optically thick [OI] 63µm
The neutral oxygen ratio [OI] 145µm /[OI] 63µm provides information about the density and temperature of the dense neutral gas if both lines are optically thin, but the 63 µm line is usually optically thick in most of ISM conditions (Tielens & Hollenbach 1985a) An alternative scenario for high values of this ratio that does not require [OI] 63µm to be optically thick is that the line emission is absorbed by less excited (cooler or more dense) neutral oxygen along the line of sight (Kraemer et al. 1998), where even small amounts of cold foreground neutral oxygen in the ground state could cause this effect (Liseau et al. 2006).However, we do not observe narrow absorption components in the [OI] 63µm spectral profile of W2246-0526, and therefore the optically thick scenario remains as the most likely cause of the observed high ratio.
An additional effect that could increase the oxygen ratio is the masering of [OI] 145µm , which may occur via population inversion within warm giant molecular clouds (Elitzur & Asensio Ramos 2006;Liseau et al. 2006), although it is not expected to affect the line emission significantly.Finally, dust obscuration can also act as a source of opacity at lower wavelengths, dimming more effectively the [OI] 63µm emission line compared to [OI] 145µm .However, we can discard this effect as a potential mechanism for faint [OI] 63µm emission since we see that the ratio of [OI] 63µm to emission lines at longer wavelengths (e.g., [CII] 158µm ) is not particularly low when compared to other galaxies (see bottom left panel of Fig. 3).The neutral oxygen ratio is shown in the y-axis of the top right panel of Fig. 3.

XDR diagnostics
The neutral carbon [CI] 609µm /[CI] 370µm ratio primarily depends on the temperature of the cold neutral gas, with higher temperatures lowering the value (Meijerink et al. 2007;Papadopoulos et al. 2022).The non-detection of [CI] 609µm in W2246-0526 sets a 3σ upper limit of [CI] 609µm /[CI] 370µm < 0.3.The presence of intense X-ray emission is expected for this range of ratios, ≲ 0.19 (Fernández-Ontiveros et al. 2016), as X-rays penetrate deeper in the cloud than UV radiation and can be an important source of heating in the neutral medium.The neutral carbon ratio is shown in the x-axis of the top right panel of Fig. 3.
Since in XDRs neutral carbon emission is not limited to a thin layer (as in PDRs) but distributed more smoothly through the ISM, the [CII] 158µm /[CI] 370µm ratio is also used to discern between PDR and XDR models (Venemans et al. 2017;Novak et al. 2019;Pensabene et al. 2021;Decarli et al. 2022Decarli et al. , 2023)).Ratios of [CII] 158µm /[CI] 370µm > 20 discard XDR models in these previous studies.Our observation of [CII] 158µm /[CI] 370µm = 58±29 for W2246-0526 favors PDR models, which appears to contradict the interpretation of the low [CI] 609µm /[CI] 370µm and [CI] 370µm /CO(7 − 6) ratios.Given that these models cover a range of parameters that may not be applicable to our source, we address this apparent tension together with the results of the Cloudy modeling in Sect.3.3.The [CII] 158µm /[CI] 370µm ratio is shown in the y-axis of the top left panel of Fig. 3.

Extreme electron density
The [NII] 122µm and [NII] 205µm emission lines originate from fine-structure transitions within the same ionization species.As such, the ratio between the two is insensitive to the metallicity of the gas and the intensity or hardness of the incident radiation field.Instead, it primarily depends on the electron density and weakly on the electron temperature of the low density ionized gas.Above a density of n e ∼ 300 cm −3 both population levels are beyond their critical densities (shown in Table 1) and become saturated.This places an upper limit to the [NII] 122µm /[NII] 205µm ratio of ∼ 10 ( Oberst et al. 2006;Spinoglio et al. 2015;Pereira-Santaella et al. 2017).
Our result of [NII] 122µm /[NII] 205µm = 13.8 ± 4.2 indicates that the warm ionized gas in W2246-0526 is most likely in the high density limit.This ratio is shown in the x-axis of the bottom right panel of Fig. 3.

High gas temperatures
The [OI] 63µm /[CII] 158µm ratio depends mainly on the density of the gas due to the different critical densities of both species.In addition, it is also sensitive to the gas temperature, as discussed in Guillard et al. (2015).A comparison of their ISM models with our observed ratio of 4.7±1.0 in W2246-0526 suggests gas temperatures of a few hundred K, which is also expected in XDRs (Wolfire et al. 2022).High [OI] 63µm /[CII] 158µm ratios can also be found in high density and temperature PDRs exposed to the high radiation fields from OB stars (Stacey et al. 1983(Stacey et al. , 1993)).As discussed in Sect.3.1.2,[OI] 63µm is likely optically thick, so the ratio with [CII] 158µm and therefore the implied gas temperatures may be even higher.This ratio is shown in the y-axis of the bottom left panel of Fig. 3.

Radiation field intensity
Recently, Harikane et al. (2020) and Algera et al. (2024) reported low [OIII] 88µm /[CII] 158µm ratios in galaxies with low dust temperatures, associating the low ratios to low ionization parameters.The [OIII] 88µm /[NII] 122µm ratio has also been proven useful as a radiation field intensity indicator (Ferkinhoff et al. 2011) and as a gas-phase metallicity diagnostic for highredshift starburst galaxies (Rigopoulou et al. 2018;Harikane et al. 2020), with the value increasing with increasing ionization parameter and decreasing with increasing metallicity.Consequently, our 3σ upper limits of [OIII] 88µm /[CII] 158µm < 0.69 and [OIII] 88µm /[NII] 122µm < 1.02 for W2246-0526 may indicate low ionization parameters, but the models used in these previous works have a more restricted parameter space than our Cloudy models.We discuss these observed ratios thoroughly in Sect.
From previous studies, a pronounced deficit (reduced lineto-FIR ratio with increasing L FIR ) is expected for [NII], [CI] and [CII] emission lines, while a less significant deficit is found for [OI] and [OIII] emission lines (e.g., Rosenberg et al. 2015;Díaz-Santos et al. 2017;Walter et al. 2018;Li et al. 2020;Pensabene et al. 2021).Figure 4 shows a clear tendency of lower line-to-FIR ratio with increasing L FIR for [CII] 158µm , [NII] 205µm , [CI] 370µm , and [CI] 609µm , with some ratios in W2246-0526 at least two orders of magnitude smaller than lower luminosity galaxies, but well in line with high-z quasars (blue points in Fig. 4).The decreasing tendencies are weaker or not significant in [OI] 63µm , [OIII] 88µm , [NII] 122µm , and [OI] 145µm .The origin of such FIR line deficits is still under study, and we do not address them in the current work.

ISM models
In order to determine the physical conditions of the ISM, we simultaneously compare all observed FIR line ratios of W2246-0526 with grids of photoionization models built using the spectral synthesis code Cloudy v17.00 (Ferland et al. 2017).To simulate astrophysical environments, Cloudy requires as an input the SED and intensity of a radiation field, as well as the chemical and dust compositions, and the geometry of a cloud of gas.Cloudy solves the equations of statistical and thermal equilibrium of the gas cloud, obtaining the physical conditions (density, temperature, ionization) across the cloud and its resulting spectrum.
With Cloudy, we model five physical parameters of the cloud and incident radiation, namely: the extinction, A V ; the metallicity of the gas, Z; its volume gas density, n H ; the X-ray to UV ratio of the input SED, α ox ; and the ionization parameter, U.The range and step of the model grid are presented in Table 5. Due to the high obscuration present in W2246-0526, we adopted a spherical model with a covering factor of unity, which corresponds to a semi-infinite plane-parallel slab with a closed geometry.
Hot DOGs are expected to be powered by the central AGN, with little contribution from their star-forming host (Jones et al. 2014;Tsai et al. 2015;Finnerty et al. 2020).Specifically for W2246-0526, the best-fit SED longward of ∼ 1µm and the total bolometric luminosity is dominated by the central quasar (Díaz-Santos et al. 2018;Fan et al. 2018;Tsai et al. 2018).Exhaustively exploring the vast landscape of Cloudy setups with fine grids is computationally expensive and beyond the scope of this paper.Given the amount of ancillary data and literature supporting that W2246-0526 is an extremely obscured quasar, we thus selected a pure AGN as the source of radiation for the Cloudy models.This corresponds to a continuum SED consisting of a sum of two components: the "Big Bump" component, a rising power law peaking at a temperature T BB , with a low-energy slope parameterized by α uv , and an infrared and a high-energy exponential cutoffs; and the X-ray component, parameterized by a slope α x between 1.36 eV and 100 keV, and decreasing with increasing frequency as f ν ∝ ν −2 .The two components are normalized by the X-ray to UV ratio α ox , defined as follows: This multicomponent continuum of the AGN template in Cloudy can be shaped by varying T BB , α uv , α x , and α ox .To introduce potential enhanced X-ray emission from the AGN as a parameter in our models we decided to vary α ox between −1.4 and −0.4 in steps of 0.2.We note that α ox becomes less negative as the ratio of X-ray to UV emission increases.The other parameters are degenerate with α ox and/or the ionization parameter, as they vary the X-ray to UV ratio and the total amount of ionization photons.Therefore, we decided to keep their default values (T BB = 10 6 K, α uv = −0.5, and α x = −1) in order to not increase the complexity and degeneracies of the fitting.An example of the AGN SED used by Cloudy is shown in Fig. 5.We run models with different values of the ionization parameter, which is defined as follows: Article number, page 9 of 20 A&A proofs: manuscript no.aanda where Q(H) is the number of hydrogen-ionizing photons, r 0 is the distance between the ionizing source and the inner part of the cloud, n(H) is the hydrogen density of the cloud, c is the speed of light, and Φ(H) is the flux of hydrogen-ionizing photons at r 0 .Increasing U, therefore, increases the ratio between hydrogenionizing photons and total hydrogen density.In our grid of models, we vary U between 10 −3.0 and 10 0.5 in steps of 10 0.25 , a range similar to other models in the literature (e.g., Moy & Rocca-Volmerange 2002;Groves et al. 2004;Nagao et al. 2006).
In addition to the intensity and SED of the radiation field, Cloudy requires the chemical composition and the density profile of the cloud of gas and dust to be specified.We assume a constant hydrogen density through the cloud, ranging between 10 2.0 and 10 5.0 cm −3 in steps of 0.25 dex, similar to other works in the literature (e.g., Meijerink et al. 2007;Feltre et al. 2016;Pensabene et al. 2021).We let both the gas and dust metallicity vary and calculate models for 0.05, 0.10, 0.25, 0.5, 1 and 2 Z ⊙ .We do not vary, however, relative elemental abundances,  2008).We note that some relative elemental abundances, specifically the N/O ratio, due to the secondary nature of nitrogen production in stellar evolution, depend on the overall, absolute metallicity (e.g., O/H) of the gas, especially at values below 0.25 Z ⊙ (Nicholls et al. 2017;Chartab et al. 2022).However, as we subsequently see, the best-fit metallicity for W2246-0526 is 0.5 Z ⊙ , still sufficiently high for relative abundances not to be a dominant factor in the models.Finally, we also vary the depth of the cloud in the form of different A V values: 10 1.0 , 10 1.5 , 10 2.0 , 10 2.5 , 10 3.0 , 10 3.5 , 10 4.0 mag.The code stops when the different A V values are reached.A cosmic ray background for the redshift of W2246-0526 is also included as it plays an important role in deep-heating the molecular cloud, as well as at low temperatures (Gaches et al. 2022).Additional modeling options were investigated and are detailed in the Appendix B, supporting our choice of parameters and model grid.
The temperature of the coldest dust component in W2246-0526 is ∼ 70 K (Díaz-Santos et al. 2016, Tsai et al. in prep.), and the cosmic microwave background (CMB) temperature is around ∼ 15 K at z ∼ 4.6.Therefore, following da Cunha et al. ( 2013), we can neglect the CMB heating in the models.
We also note that, in order to avoid introducing additional complexity into the models, we do not attempt to fit the continuum emission.Modeling the CO SLED of W2246-0526 is beyond the scope of this paper as well, and a separate analysis of the available CO line measurements will be addressed in a future publication.

Line ratio diagnostics within Cloudy
In Sect.3.1.2we proposed optically thick [OI] 63µm as the most likely cause of the observed [OI] 145µm /[OI] 63µm ratio.The comparison with our Cloudy models reinforce this scenario, as we can only reproduce the observed value with column densities above 10 23 cm −2 .
As mentioned in Sect.3.1.3ratios of [CII] 158µm /[CI] 370µm > 20 discard XDR models in previous studies, but we find in our Cloudy models that this applies only to log(U) < −1.0, while higher ratios are reproducible with higher ionization parameters.Our observation of [CII] 158µm /[CI] 370µm = 58 ± 29 for W2246-0526 points toward an extreme ionization of the ISM in XDR.Another mechanism that can enhance the [CII] 158µm /[CI] 370µm ratio is [CII] excitation by shocks and turbulence driven by jets (Appleton et al. 2018), which can boost the ratio up to ∼ 40%.The observed [CII] 158µm kinematics in W2246-0526 show high turbulence, suggesting the central quasar might be blowing isotropically its ISM in large-scale outflows (Díaz-Santos et al. 2016).We do not address this case as it falls beyond the scope of the current study.
The neutral carbon [CI] 609µm /[CI] 370µm ratio supports the XDR scenario, as we see in our Cloudy models that increasing the X-ray to UV ratio decreases the neutral carbon line ratio (Fig. 6).We find that it also decreases with increasing density and U, and with decreasing A V , and it does not depend on the metallicity.Given our upper limit, the models favor higher α ox , U and densities, and lower A V .Regarding radiation field intensity indicators involving the [OIII] 88µm emission line (Sect.3.1.6),previous studies find lower [OIII] 88µm /[CII] 158µm and lower [OIII] 88µm /[NII] 122µm for lower ionization parameters.In our Cloudy models, however, this dependence for [OIII] 88µm /[CII] 122µm is seen only in the log(U) range ∼ −3 to ∼ −2 that they cover in their models, while low ratios are also found in our models for much higher ionization intensities (and therefore much higher dust temperatures, as is the case for Hot DOGs).Lower [OIII] 88µm /[NII] 122µm ratios for lower ionization parameters are seen in our models only in low density (log(n H ) ≤ 3 cm −3 ) and low α ox (≤ −1.0).In other regimes, the ratio mainly decreases with increasing α ox , increasing density, and increasing ionization parameter.Perhaps counter-intuitively, the emitting region where [OIII] is produced is actually smaller at high radiation field intensities and X-ray luminosities, as oxygen is ionized to even higher levels (see middle panel of Fig. 9), while the emitting region of [NII] becomes larger with increasing X-ray ionization due to the change in the ISM structure.Our 3σ upper limit [OIII] 88µm /[NII] 122µm < 1.02 for W2246-0526 strongly indicates again a high-density ISM being affected by high X-ray ionization.

Best-fit model
The line ratios described above are independent, powerful tools to derive specific properties of a particular phase of the ISM gas.Modeling various ISM phases together can be challenging, as they may not be, for instance, uniform in density, homogeneous or isotropic, and the emission lines can be produced in different regions in the galaxy nucleus or at different physical scales.Nonetheless, even when assumptions on the geometry and spatial distribution of the media are simplified (sometimes unavoidably due to a lack of observational information), models can still be used to provide insights into the average conditions of the multiphase gas that is exposed to the radiation field of the central quasar.
To this end, we compare all the observed lines with our Cloudy model grid using a Bayesian framework from which we can estimate the marginalized 1D and 2D posteriors of the probability distribution functions (P = e −χ 2 /2 ) for each of the considered parameters.An uncertainty of 15% is assumed for the model ratios to account for the sparsity of the parameter grid.As we see in the corner plot shown in Fig. 7, the most probable solution (and in parenthesis the posterior mean with 68% confidence intervals) corresponds to a hydrogen density log(n H ) = 3.5 cm −3 (3.44 ± 0.37), an ionization parameter log(U) = −0.5 (−0.47 ± 0.62), an extinction A V = 316 mag (222 +546 −158 ), an X-ray to UV ratio α ox = −0.8(−0.78 ± 0.29), and a metallicity Z = 0.5 Z ⊙ (0.49 +0.67  −0.29 ).The best fit between model and observations has a reduced chi-square χ 2 ν = 3.34, and it is shown in Fig. B.1.For fixed FWHM, the ratios are very similar and the reduced chi-square is χ 2 ν = 4.12.The average gas density is well constrained and comparable to values obtained for high-z quasars (Wang et al. 2016;Venemans et al. 2017;Novak et al. 2019;Yang et al. 2019;Pensabene et al. 2021;Meyer et al. 2022).A density of log(n H ) = 3.5 cm −3 is above the critical densities of the lines arising from ionized gas (Table 1), indicating they are all near or above the high density limit.Densities higher than n cr saturate the emission, decreasing thus the efficiency of the lines as coolants of the ISM.The extinction has a very high value; an A V of 10 2.5 mag corresponds to a column density close to 10 24 cm −2 , in agreement with the column densities expected for Hot DOGs (Stern et al. 2014;Piconcelli et al. 2015;Assef et al. 2015Assef et al. , 2016Assef et al. , 2020)), as well as for obscured high-z quasars (Vito et al. 2019;Pensabene et al. 2021;Meyer et al. 2022).A metallicity of 0.5 Z ⊙ suggests a fairly evolved and enriched ISM at z ∼ 4.6, and falls within the range of metallicities (0.2 to 2.0 Z ⊙ ) reported by Carvalho et al. (2020) for AGNs at both low and high redshifts (z = 0 − 7).The ionization parameter is also well constrained and significantly higher than other high-z AGNs (−2.2 < log(U) < −1.4 in Nagao et al. 2006).The X-ray to UV ratio α ox is again extreme among the population of quasars.Miller et al. (2011) finds values of α ox varying between −2.0 and −0.8.We can see in the corner plot that U is degenerate with α ox , in the sense that the model is able to reproduce the line ratios with almost the same probability by increasing α ox if U is reduced.Because of this, α ox is unbound at the high end and we can only put a lower limit of α ox ≥ −0.8.
The lower limit in α ox and the high value of U correspond to intense X-ray emission from the central AGN.Observations with XMM-Newton (Vito et al. 2018) set an upper limit to the X-rays that is a factor of ∼ 2 lower than what is predicted by the best model.We address and discuss this result in Sect.4.2.2.
In Fig. 8, we present the emission of each FIR fine-structure line at each layer of the cloud of gas for the best model.The structure of the gas cloud, also seen in Fig. 9 for the different ion abundances, corresponds to that of an XDR: the [NII] lines continue to arise even where hydrogen is mostly neutral, and the neutral carbon emission is distributed through the ISM and not limited only to a thin layer as expected in PDRs.We find that ∼90% of [CII] emission is emitted after the ionization front (i.e., in the non-ionized gas).While this fraction depends on the parameters of the model, the difference with the estimation in Sect.3.1.1(80%) clearly shows that in sources with important contribution from XDRs, using a scaling based on PDR models leads to different results.

Discussion
In Sect.3.4 we presented the best-fit model parameters and their posterior distributions and found that, even with the limitations imposed by the simple geometry used to model the ISM in the vicinity of W2246-0526's quasar, Cloudy is able to broadly reproduce the measured emission line ratios relatively well.However, we also find tensions between some of the predictions of the best model and other independent observations available for W2246-0526.In this section, we discuss potential explanations for those discrepancies, address the implications for the nondetections and multiphase gas budget, and put our results in the context of other populations of quasars at high redshift.

Comparison with other high-z quasars and local AGNs and ULIRGs
We compare the FIR emission line ratios measured for W2246-0526 to those found in local AGNs and ULIRGs (ultra-luminous infrared galaxies) as well as high-redshift quasars reported in the literature.This compilation of sources is shown in Fig. 3, together with the values of W2246-0526 and the Cloudy model grids for varying n H and U.
We note that the model grids shown in Fig. 3 are only a small part of the parameter space we covered by the Cloudy models.We only use the best-model values obtained for the metallicity, α ox , and A V , which are representative only of models with high X-ray to UV emission and high extinctions.Figure B.3 in the Appendix illustrates the effect of varying these two parameters (α ox and A V ).We also caution that our modeling considers an AGN as the only source of radiation, while in some quasar populations there may be an additional significant contribution to the ionization from star formation in the host galaxy.Therefore, the main goal of including these galaxies along with the grids of models is not to characterize them, but to put W2246-0526 in the context of other AGN-powered galaxies.
From a modeling point of view, W2246-0526 stands out in terms of the large X-ray and radiation field intensities needed to reproduce the observed line ratios when compared to other AGNs and quasars at both low and high redshift.This is reflected in some specific individual line diagnostics, where we find that W2246-0526 is at the edge of the distribution of other populations, showing extremely high values of the [OI] 63µm /[CII] 158µm and [NII] 122µm /[NII] 205µm ratios, and very low values of the [CI] 609µm /[CI] 370µm ratio, all of them independently suggesting a very dense and highly ionized ISM with X-ray dominated heating.Therefore, the unusual ISM conditions in W2246-0526 can serve as a benchmark for comparisons with other quasar populations at cosmic noon and beyond that may be less extreme.radiative pressure will create a pressure gradient in ionized regions, that can make U to saturate, implying that values of U higher than 10 0.5 (the upper value in our grid of Cloudy models) are highly unlikely.In fact, such high (or higher) values have never been reported in the literature when modeling sources.

X-ray emission
Since hyper-luminous AGNs and quasars are expected to be powerful X-ray emitters (e.g., Lusso et al. 2010), we introduced α ox as an additional parameter in order to search for a more physical solution.The result is that α ox and U are degenerate, as shown in the 2D marginalized posterior distributions presented in Fig. 7.In this case, however, the value of α ox is unbound, with a lower limit of ∼ −0.8, but the value of U is more constrained and physical, with the best fit peaking at ∼ 10 −0.5 .
In order to quantitatively assess the model assumptions that better reproduce the data, that is, that with fixed α ox = −1.4(henceforth model A), or that with α ox as a free parameter (henceforth model B), we use the χ 2 statistic defined as , where d i and ∆d i are the observed ratios and their uncertainties, and m i and ∆m i are the model ratios and their uncertainties, which account for the sparsity of the parameter grid.We obtain χ 2 (A) = 15.2, χ 2 (B) = 6.7, χ 2 ν (A) = 5.1, and χ 2 ν (B) = 3.3, with χ 2 ν defined as χ 2 ν = χ 2 n−k , n being the number of data points (number of ratios), and k the number of free parameters of each model.Both χ 2 and χ 2 ν are smaller for model B. To investigate if the difference is significant, we use the Bayesian Information Criteria (BIC), defined as BIC = log(n)k − 2log( L), with log( L) being the maximum log-likelihood of the model.The BIC statistic takes into account the performance of a model penalizing its complexity.A lower BIC means less information is lost, and hence the model is better.Following Burnham & Anderson (2004), we can compute the posterior model probability as p i = exp(−1/2∆BIC i ) i (exp(−1/2∆BIC i )) , with ∆BIC i being the difference of the BIC of model i with the minimum BIC of all models.The probability of model A over B is p(A) = 0.15, and of model B over A is p(B) = 0.85.Therefore, model B is statistically better than model A, indicating that it is preferable to vary α ox .

X-ray non-detection
Vito et al. ( 2018) reported a non-detection of X-rays in W2246-0526, as part of a study of X-ray properties in Hot DOGs.To compare with our best-fit model we extract, from its associated input SED (Fig. 5), the X-ray flux between 2-10 keV transmitted through the cloud and the bolometric incident flux, with a ratio of f (2−10 keV) f bol ≃ 10 −2 .The 3σ upper limit of X-ray ≃ 0.5 × 10 −2 .Comparing the two ratios, the best model predicts a factor of two more X-rays compared to the upper limit of the observations.Given the conversion factor A V /N H provided by Cloudy and the best fit value for A V we estimated a hydrogen column density N H = (1.2+3.1 −0.8 ) × 10 24 cm −2 .To be able to obscure the model X-rays below the observed limit (a factor of two less), we would need an extra N H = 4.6 × 10 23 cm −2 along the line of sight, which is within 0.5σ from the best-fit model given the uncertainties in N H . Therefore, the predictions of our models indicate that we may be close to detecting the Xray emission from W2246-0526, highlighting the importance of future, deeper X-ray observations.
A non-detection of hard X-rays in heavily obscured quasars is not surprising (see, e.g., Andonie et al. 2022).In Piconcelli et al. (2015), the Hot DOG W1835+4355 is also extremely obscured in X-rays and in the Compton thick regime.Stern (2015) presented a mid-infrared (MIR) X-ray relation that broadly agrees with our best model predictions estimate.Namely, given the MIR luminosity of W2246-0526, νL ν (6µm) ≃ 2 × 10 14 L ⊙ (see Fig.  2017) is predicted to have around ∼ 3 times more X-ray emission than observed, when using the Stern (2015) relation.To reconcile the difference, they suggest either X-ray weakness or significantly Fig. 9: Ionized, neutral, and molecular hydrogen and carbon (top), oxygen (middle), and nitrogen (bottom) number density fractions as a function of the cloud depth, measured in hydrogen column density, for the Cloudy model that best fits the W2246-0526 observations.We note the negligible fraction of molecular gas and CO (top), even at depths as large as 10 24 cm −2 .more obscuration than estimated.Given our model prediction, and the fact that a large fraction of MIR-selected AGNs are Compton thick (Carroll et al. 2023), we favor the second scenario for W2246-0526 and propose that the reason for a nondetection of W2246-0526 in X-rays from the literature is likely due to additional obscuring columns of gas along the line of sight, not associated with the cloud we are modeling.

Cold molecular gas
The amount of CO predicted by the models is extremely low (Fig. 9), as expected from the high ionization and especially the intense X-ray emission implied by the best-fit model.Xray photons cause the CO to dissociate much deeper in the gas cloud than in PDRs where photo-dissociation is dominated by far-UV photons.AGN X-ray dissociation of CO has been observed, as reported in Kawamuro et al. (2020).The amount of single and double ionized species is also very high in the presence of intense ionization and X-rays, and ions can also destroy CO molecules via charge transfer (Ferland et al. 2013).
Interestingly, Díaz-Santos et al. ( 2018) used JVLA to observe the CO(2-1) emission line transition and reported the existence of large amounts of molecular gas in W2246-0526, which is in contrast with our model predictions.In particular, the observational CO(2-1)/[CII] 158µm luminosity ratio in W2246-0526 is ∼ 5 × 10 −3 , versus ∼ 3 × 10 −12 from the best-fit model.Moreover, the measured CO(7-6)/CO(2-1) luminosity ratio is ∼ 17, but the best-fit model yields a value of ∼ 100.This agrees with the findings of Lee et al. (in prep.), who studied a sample of ten Hot DOGs and identified that molecular gas is, in general, remarkably excited and abundant.However, it is key to note that the angular resolution of the JVLA CO(2-1) observations of W2246-0526 is much lower, 2.47 ′′ × 2.01 ′′ , than the aperture we use to extract the fluxes of the fine-structure lines and CO(7-6), suggesting that the cold molecular gas component detected in W2246-0526 may be located and spread over a region at much larger scales than those being radiatively affected by the central quasar.Critically, this would be also in agreement with our findings regarding the non-detection of X-ray emission, as we discussed in the previous section.We note, however, that the presence of multiple gas phases and the likely not very homogeneous ISM distribution are not accounted by the relative simplicity of the Cloudy models and could also affect the CO luminosity and X-ray absorption.A detailed discussion regarding the properties of the molecular gas component of the Hot DOG population will be addressed in a future work via CO SLED modeling.

Physical size and total gas mass of the cloud
Díaz-Santos et al. ( 2018) calculated the mass of molecular hydrogen, M(H 2 ), in W2246-0526 using the CO(2-1) emission line.They estimate a M(H 2 ) = 7 ± 1.5 × 10 10 M ⊙ , assuming an α CO conversion factor for ultra luminous infrared galaxies (ULIRGs; L IR ≳ 10 12 L ⊙ ) of α CO,ULIRG = 0.8 M ⊙ [K km s −1 pc 2 ] −1 .Assuming that the geometry we have selected for the Cloudy modeling is a reasonable approximation, we can estimate the hydrogen mass (ionized + neutral, since the fraction of H 2 is negligible; see Fig. 9) contained in the cloud using its volume and the best-fit value of n H .To calculate the volume, we derive the inner radius of the cloud by comparing the observed bolometric luminosity of W2246-0526 with the intrinsic bolometric flux (black solid line in Fig. 5) obtained from the best-fit model.We use the observed, tens-of-kiloparsecs integrated bolometric luminosity as an upper limit for our modeled SED, as the latest only considers the emission from the central quasar.This puts an upper limit on the inner radius of the cloud of r 0 ≲ 670 pc.Cloudy also reports the size of the cloud, which is ∼ 120 pc.Using the volume of this shell and n H , we estimate an upper limit for the hydrogen mass of the cloud (again, ionized + neutral) of M H < 0.6 × 10 11 M ⊙ , which is smaller than the amount of molecular gas derived from the CO(2-1) line.

Implications for the [OIII] 88µm and [CI] 609µm
non-detections Intense [OIII] 88µm line emission has been successfully detected in high-z galaxies up to z ∼ 9 (e.g., Harikane et al. 2020;Witstok et al. 2022), which in principle would make the non-detection of the line in W2246-0526 surprising.However, [OIII] 88µm emission arises most effectively from mildly ionized, diffuse gas (Lebouteiller et al. 2012;Witstok et al. 2022) while instead, our best-fit model favors very high hydrogen densities, U, and α ox .
In these XDR conditions most of the oxygen in the ionized gas is in higher ionization states (see middle panel on Fig. 9) and [OIII] 88µm emission becomes highly suppressed, thus likely explaining the non-detection.This may have critical implications for the detection of this line in sources at the epoch of reionization, where the average intensity of the radiation field in galaxies should be stronger than at later cosmic times (Hashimoto et al. 2019;Algera et al. 2024).
Regarding the non-detection of the [CI] 609µm line, the Cloudy models clearly show that the combination of high ionization and X-rays leads to less [CI] 609µm emission, and an overall lower abundance of neutral carbon with respect to ionized carbon (see top panel of Fig. 9).Following Papadopoulos & Greve (2004), we can also estimate the expected flux of [CI] 609µm using the value of M(H 2 ) derived from the CO(2-1) observations mentioned above.Based on this, we obtain a flux of S ([CI] 609µm ) ∼ 0.2 Jy km s −1 , which is consistent with the 3σ observational upper limit (< 0.24 Jy km s −1 ).Given these estimations, deeper observations with ALMA of [CI] 609µm should achieve this sensitivity and successfully detect the emission line.

Summary and conclusions
In this work we have presented multiline ALMA observations of the extremely luminous, hot, dust-obscured galaxy W2246-0526 at redshift z = 4.6.We detected the FIR fine-structure lines [OI]  With the goal of characterizing the average physical conditions of the ISM of W2246-0526, we created a large grid of models with the photoionization code Cloudy.The main results are summarized as follows: -The model that best fits the observations indicates that W2246-0526 stands out from other quasars at high and low redshift in terms of its very high X-ray to UV ratio (α ox ≥ −0.8) and ionization parameter (U = 10 −0.5 ).The average density of the cloud is comparable to other high-z quasars (∼ 3 × 10 3 cm −3 ), and with a metallicity close to 0.5 Z ⊙ .The extinction is very high (A V ∼ 300 mag), with hydrogen column densities in the Compton thick regime.-The observed emission line ratios require high gas temperature, intense ionization, and X-ray emission in the ISM, conditions typical of XDRs.This also explains the nondetections of [OIII] 88µm and [CI] 609µm emission lines, which may be critical in the observation of these lines in other highz AGNs and galaxies in general.-The models predict an X-ray flux from W2246-0526 larger (although compatible within the model uncertainties) than the upper limit set by a non-detection with XMM-Newton.Also, the models predict almost no CO emission, which is in contrast with previous works reporting a strong JVLA detection of CO(2-1) on scales of a few tens of kiloparsecs.When combined, these two pieces of evidence suggest the presence of an additional molecular gas component located farther away from the central quasar, distributed over spatial scales much larger than the region we are modeling with Cloudy.
This study demonstrates the need for multiline FIR observations to properly characterize the multiphase ISM properties of galaxies and quasars at high redshift.In particular, W2246-0526 serves as a benchmark for the extreme conditions that may be potentially present in dust-obscured quasar populations at and beyond cosmic noon.them, the metallicity is set to 1 Z ⊙ , and the A V is 100 mag.The choice of density profiles in ISM modeling can highly affect the molecular gas and CO emission predictions.We selected a constant density profile, but other options may be more appropriate for modeling the molecular phase as discussed in Popping et al. (2019).We investigated the possibility of using a constant gas pressure through the cloud instead of a constant density, with an AGN with α ox = −1.4 as the ionization source (bottom right in

Fig. 2 :
Fig. 2: Spectra (left) and intensity map (right) for each observed line in W2246-0526.The blue dashed line indicates the zero flux density in each spectrum, the red line is the Gaussian fit for each emission line, and the two orange lines indicate the 2.75σ of each Gaussian fit that is used to integrate the flux of the emission line directly on the spectra.In the intensity maps, the blue dot-dashed line shows the aperture used to extract each spectrum, and the clean beam is shown at the bottom left.The zero-flux intensity level is shown as a white, dotted contour.Solid contours indicate [2.5, 3, 2 (4+n)/2 ]σ levels (with σ being the rms of the map, and n = [0, 1, 2, ...]).Dashed contours show negative flux at the same absolute levels.

Fig. 3 :
Fig. 3: Line ratio diagrams from the Cloudy models for the metallicity, α ox and A V of the model with the lowest χ 2 ν described in Sect.3.4.The observed line ratios for W2246-0526 are indicated with a red square, and a red triangle for upper limits.The color grid spans the full range of values explored in this work for the ionization parameter and the density, as shown in the legend.The red open circle indicates the best-fit value for our source (i.e., lowest χ 2 ν model).Literature values for high-z quasars (Weiß et al. 2003; Uzgil et al. 2016; Venemans et al. 2017; Walter et al. 2018; Lee et al. 2019; Novak et al. 2019; Yang et al. 2019; Li et al. 2020; Pensabene et al. 2021; Meyer et al. 2022; Decarli et al. 2022, 2023) are indicated with blue stars, or solid blue lines with shadow uncertainties if only one ratio is available.Literature values for local AGNs are indicated with open green crosses for Fernández-Ontiveros et al. (2016) and open magenta plus signs for Spinoglio et al. (2015), while local ULIRGs from Rosenberg et al. (2015) are indicated with open brown diamonds.
3.3.The [OIII] 88µm /[CII] 158µm and [OIII] 88µm /[NII] 122µm ratios are shown, respectively, in the x-axis of the top left panel and in the y-axis of the bottom right panel of Fig. 3.

Fig. 5 :
Fig. 5: Input AGN SED (black solid line) and reprocessed radiation by the cloud of gas and dust (red dotted line), for the best-fit model for W2246-0526 described in 3.4 with log(U) = −0.5 and α ox = −0.8.

Fig. 6 :
Fig. 6: Neutral carbon [CI] 609µm /[CI] 370µm ratio from our Cloudy models as a function of the ionization parameter.Each color represents varying X-ray to UV ratios.The metallicity is fixed to 1Z ⊙ , the A V to 10 mag, and the hydrogen density to 10 3 cm −3 .The gray dashed line (denoted by * in the legend) represents α ox = −0.4but with an A V = 100 mag.

Fig. 7 :
Fig. 7: Corner plot showing the marginalized 1D and 2D posterior probability distributions (PPDs) for the parameters in the Cloudy models used to fit the emission line ratios of W2246-0526.The orange solid lines indicate the parameter values of the best model (maximum likelihood estimation, MLE), the blue dashed line the median of each parameter marginalized over the rest, and the blue shaded region the 68% confidence interval.The y-axes of the 1D PPDs show the probabilities with their sum normalized to one for each parameter.

Fig. 8 :
Fig. 8: Cumulative line intensity (top panel) and line emissivity (bottom panel) at each depth in the cloud, for the Cloudy model that best fits the W2246-0526 observations.Both panels show each modeled FIR finestructure line as a function of the cloud depth, measured in hydrogen column density.The HII fraction is included to visualize the ionization front, where it drops sharply.

Fig
Fig. B.3: Line ratio diagrams from Cloudy AGN models for solar metallicity, α ox = −1.4 with A V = 10 1 mag (labeled as A, grid with dashed lines), and α ox = −0.8 with A V = 10 2.5 mag (labeled as B, grid with dash-dotted lines).Both grids span the full range of values explored in this work for the ionization parameter and the density, as shown in the legend and Fig. 3.The observed line ratios for W2246-0526 are indicated with a red square, or a red triangle for upper limits.
Fig. B.4).In this case, the density of the cloud n H is the initial density at the inner part of the cloud.We also explored three different stellar SEDs as the input source instead of an AGN, using single stellar populations from Mollá et al. (2009) of ages 2 Myr, 7 Myr, and 20 Myr (with initial mass function from Chabrier 2003, and solar metallicity.)The density for the stellar populations in these cases is constant through the cloud.As seen in Fig.B.4, none of the test grids fall in the region of the parameter space where the observed emission ratios lay, a finding that favors a constant density modeling with an AGN as a radiation source.

Fig
Fig. B.4: Same line ratios as in Fig. 3, but for smaller Cloudy grids varying the ionization parameter and the density of the cloud to test stellar populations (a single 2 Myr population in the top left, 7 Myr in the top right, and 20 Myr in the bottom left) as a source of radiation, instead of an AGN.In the bottom right we test a constant pressure density profile instead of the constant density one.

Table 1 :
Fundamental physical parameters for the targeted lines Columns correspond to the emission line, transition, central wavelength, central frequency, formation potential (FP), electron (n cr,e ), and Hydrogen (n cr,H ) critical densities, and energy difference between upper and lower level(E ul

Table 2 :
Details of the ALMA observations