© The Authors 2024

1. 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 cosmological 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¹³ 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¹³ 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¹⁴ L_⊙, making it one of the most luminous galaxies known to date (Tsai et al. 2015, 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⁹ 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, fine-structure lines from [OI]_63 μm to [CI]_609 μm can be observed. Previous 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 fine-structure 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).

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 (H2). 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.

In this paper, we analyze ALMA observations of the [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 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: Sect. 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₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_M = 0.3. For reference, in this cosmology, 1″ is equivalent to 6.54 kpc at z = 4.601.

2. Observations, data reduction, and emission line measurements

This work uses data from six different ALMA projects observed during Cycles 3 through 9, with a total on-source observation time of 19 h with the ALMA main 12-m array. These observations were conducted between June 2016 and October 2022, and targeted the following FIR emission lines: [OI]_63 μm (band 10, one of the few detections in the literature to date), [OIII]_88 μm (band 9), [NII]_122 μm (band 8), [OI]_145 μm (band 7), [CII]_158 μm (band 7), [NII]_205 μm (band 6), [CI]_370 μm (band 4), and [CI]_609 μm (band 3). The properties of each transition are listed in Table 1, and the details of the observations are listed in Table 2.

We reduced and calibrated the data using the Common Astronomy Software Applications (CASA; McMullin et al. 2007). We processed the data with the ScriptForPI.py provided through the ALMA science archive to generate the measurement 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.

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. 2. continued.

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 companion 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.

As shown in Table 3, all lines have a similar FWHM around 600–700 km s⁻¹, with the exception of [NII]_122 μm, [NII]_205 μm, and [CI]_370 μm, with ∼800 km s⁻¹, ∼900 km s⁻¹, and ∼350 km s⁻¹ 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⁻¹ 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⁻¹ in Béthermin et al. 2020), optical quasars (∼385 km s⁻¹ in Decarli et al. 2018), or infrared quasars (from ∼200 km s⁻¹ to ∼500 km s⁻¹ 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):

$$\begin{array}{c}\hfill L_{\mathrm{line}}[L_{\odot }]=1.04\times {10}^{-3}\times S_{\mathrm{line}}\mathrm{\Delta }\upsilon {\nu }_{\mathrm{rest}}{(1+z)}^{-1}{D}_{\mathrm{L}}^2,\end{array}$$(1)

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⁻¹; ν_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⁻¹ 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.

3. 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.

3.1. 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⁻¹ 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.

Fig. 3. Line ratio diagrams from the CLOUDY models for the metallicity, αox and AV of the model with the lowest ${\chi }_{\nu }^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 ${\chi }_{\nu }^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.1.1. 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³ cm⁻³). 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 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.

3.1.2. 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). For W2246–0526, we obtain [OI]_145 μm/[OI]_63 μm = 0.14 ± 0.04. The [OI]_63 μm line can become optically thick at N_H > 2 × 10²¹ cm⁻² (Tielens & Hollenbach 1985b), implying that a value of the [OI]_145 μm/[OI]_63 μm ratio of ∼0.10 or higher is indicative of optically thick [OI]_63 μm emission.

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.

3.1.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.

Together with the [CI]_370 μm observations, the CO J = 7 − 6 transition is also detected in W2246–0526 (Fig. 2), with a luminosity of 0.59 ± 0.04 × 10⁹ L_⊙. CO(7–6) is a warm dense gas tracer (e.g., Carilli & Walter 2013). Mid-J CO transitions such as J = 7 − 6 are more intense in XDRs (e.g., Meijerink et al. 2007; Vallini et al. 2019), and indicative of AGN heating. Comparing with the compiled galaxies and models of Valentino et al. (2020) and Hagimoto et al. (2023), the very low ratio of [CI]_370 μm/CO(7 − 6) = 0.19 ± 0.07 and low neutral carbon luminosity for W2246–0526 points toward a high-density and highly ionized ISM, compatible with the existence of XDRs.

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. 2022, 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.

3.1.4. 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⁻³ 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.

3.1.5. 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, 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.

3.1.6. 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 high-redshift 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. 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.

3.1.7. FIR line deficits

In Fig. 4 we present the line-to-FIR_{[42 − 122 μm]} ratio as a function of the FIR_{[42 − 122 μm]} luminosity (L_FIR) for all the targeted fine-structure emission lines in W2246–0526. In addition, for comparison we show a sample of high-z quasars gathered from the literature (Weiß et al. 2003; Uzgil et al. 2016; Lu et al. 2017; Venemans et al. 2017; Walter et al. 2018; Lee et al. 2019; Novak et al. 2019; Li et al. 2020; Pensabene et al. 2021; Meyer et al. 2022), local AGNs and ULIRGs from Farrah et al. (2013), Pereira-Santaella et al. (2013, 2017), Rosenberg et al. (2015), Herrera-Camus et al. (2018), and local LIRGs from Díaz-Santos et al. (2017). We note that for galaxies without a measurement of the FIR luminosity (Farrah et al. 2013; Pereira-Santaella et al. 2013, 2017), we assumed L_{IR[8 − 1000 μm]} = 6.3 × L_{FIR[42 − 122 μm]}, the scaling used for W2246–0526 in Díaz-Santos et al. (2016).

Fig. 4. Line-to-FIR[42 − 122 μm] luminosity ratio as a function of FIR[42 − 122 μm] luminosity for [OI]63 μm, [OIII]88 μm, [NII]122 μm (top row), [OI]145 μm, [CII]158 μm, [NII]205 μm (middle row), [CI]370 μm, and [CI]609 μm (bottom row). W2246–0526 is indicated with a red square. The high-z quasars are taken from Weiß et al. (2003), Uzgil et al. (2016), Lu et al. (2017), Venemans et al. (2017), Walter et al. (2018), Lee et al. (2019), Novak et al. (2019), Li et al. (2020), Pensabene et al. (2021), Meyer et al. (2022). Upper limits for high-z quasars and W2246–0526 are shown with downward arrows. We note that the only high-z quasar shown for [OI]63 μm and [CI]609 μm corresponds to the Cloverleaf quasar (Weiß et al. 2003; Uzgil et al. 2016).

From previous studies, a pronounced deficit (reduced line-to-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.

3.2. 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_ν ∝ ν⁻². The two components are normalized by the X-ray to UV ratio α_ox, defined as follows:

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{ox}}=\frac{log(L_{2\mathrm{keV}}/L_{2500\mathrm{\AA }})}{2.605}·\end{array}$$(2)

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⁶ 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.

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 Sect. 3.4 with log(U) = − 0.5 and αox = −0.8.

We run models with different values of the ionization parameter, which is defined as follows:

$$\begin{array}{c}\hfill U\equiv \frac{Q(\mathrm{H})}{4\pi r_0^{2}n(\mathrm{H})c}\equiv \frac{\mathrm{\Phi }(\mathrm{H})}{n(\mathrm{H})c},\end{array}$$(3)

where Q(H) is the number of hydrogen-ionizing photons, r₀ 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₀. Increasing U, therefore, increases the ratio between hydrogen-ionizing 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⁻³ 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, and assume the default gas-phase values provided by CLOUDY for the ISM (log(N/O) = − 0.6), with the default grains and polycyclic aromatic hydrocarbons (PAH) size distribution from Abel et al. (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.

3.3. Line ratio diagnostics within CLOUDY

In Sect. 3.1.2 we 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²³ cm⁻².

As mentioned in Sect. 3.1.3 ratios 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.

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 1 Z⊙, the AV to 10 mag, and the hydrogen density to 103 cm−3. The gray dashed line (denoted by * in the legend) represents αox = −0.4 but with an AV = 100 mag.

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⁻³) 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. 7), 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.

Fig. 7. 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 1024 cm−2.

3.4. 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. 8, 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.44 ± 0.37), an ionization parameter log(U) = − 0.5 (−0.47 ± 0.62), an extinction $A_{\mathrm{V}}=316\mathrm{mag}({222}_{-158}^{+546})$, an X-ray to UV ratio α_ox = −0.8 (−0.78 ± 0.29), and a metallicity $Z=0.5Z_{\odot }(0.{49}_{-0.29}^{+0.67})$. The best fit between model and observations has a reduced chi-square ${\chi }_{\nu }^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 ${\chi }_{\nu }^2=4.12$.

Fig. 8. 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.

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⁻³ 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²⁴ cm⁻², in agreement with the column densities expected for Hot DOGs (Stern et al. 2014; Piconcelli et al. 2015; Assef et al. 2015, 2016, 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–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. 9, 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. 7 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.

Fig. 9. 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 fine-structure 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.

4. 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 non-detections and multiphase gas budget, and put our results in the context of other populations of quasars at high redshift.

4.1. 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.

4.2. X-ray emission

4.2.1. α_ox as a free parameter

With the default CLOUDY value of α_ox = −1.4, the ionization parameter, U, for the best model is extremely high and unconstrained (Fig. B.2). According to Yeh & Matzner (2012), high 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.

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. 8. 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 χ² statistic defined as ${\chi }^2=\sum \frac{{(d_i-m_i)}^2}{\mathrm{\Delta }{d}_i^2+\mathrm{\Delta }{m}_i^2}$, 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 χ²(A) = 15.2, χ²(B) = 6.7, ${\chi }_{\nu }^2$(A) = 5.1, and ${\chi }_{\nu }^2$(B) = 3.3, with ${\chi }_{\nu }^2$ defined as ${\chi }_{\nu }^2=\frac{{\chi }^2}{n-k}$, n being the number of data points (number of ratios), and k the number of free parameters of each model. Both χ² and ${\chi }_{\nu }^2$ are smaller for model B. To investigate if the difference is significant, we use the Bayesian Information Criteria (BIC), defined as $\mathrm{BIC}=log(n)k-2log(\widehat{L})$, with $log(\widehat{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=\frac{\exp (-1/2\mathrm{\Delta }{\mathrm{BIC}}_i)}{{\sum }_i(\exp (-1/2\mathrm{\Delta }{\mathrm{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.

4.2.2. 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 $\frac{f_{(2-10\mathrm{keV})}}{f_{\mathrm{bol}}}\simeq {10}^{-2}$. The 3σ upper limit of X-ray emission between 2–10 keV reported by Vito et al. (2018) is < 6.9 × 10⁴⁵ erg s⁻¹. The bolometric luminosity of W2246–0526 is 3.6 × 10¹⁴ L_⊙ = 1.4 × 10⁴⁸ erg s⁻¹, which translates into a ratio of $\frac{L_{(2-10\mathrm{keV})}}{L_{\mathrm{bol}}}\simeq 0.5\times {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_{\mathrm{H}}=(1.2_{-0.8}^{+3.1})\times {10}^{24}$ cm⁻². 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²³ cm⁻² 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 X-ray 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¹⁴ L_⊙ (see Fig. S4 in Díaz-Santos et al. 2018), the relation predicts an X-ray luminosity of L(2 − 10 keV) = 6.4 × 10⁴⁵ erg s⁻¹, compatible with the observational upper value. Interestingly, the Hot DOG WISE J1036+0449 analyzed by Ricci et al. (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 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 non-detection 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.

4.3. Multiphase gas budgets

4.3.1. Cold molecular gas

The amount of CO predicted by the models is extremely low (Fig. 7), as expected from the high ionization and especially the intense X-ray emission implied by the best-fit model. X-ray 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⁻³, versus ∼3 × 10⁻¹² 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.

4.3.2. Physical size and total gas mass of the cloud

Díaz-Santos et al. (2018) calculated the mass of molecular hydrogen, M(H₂), in W2246–0526 using the CO(2–1) emission line. They estimate a M(H₂) = 7 ± 1.5 × 10¹⁰ M_⊙, assuming an α_CO conversion factor for ultra luminous infrared galaxies (ULIRGs; L_IR ≳ 10¹² L_⊙) of α_CO, ULIRG = 0.8 M_⊙ [K km s⁻¹ pc²]⁻¹. 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₂ is negligible; see Fig. 7) 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₀ ≲ 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¹¹ M_⊙, which is smaller than the amount of molecular gas derived from the CO(2–1) line.

4.3.3. Neutral carbon mass

Following Weiß et al. (2003), we can derive the neutral carbon mass from the [CI]_370 μm line luminosity via

$$\begin{array}{c}\hfill M_{\mathrm{CI}}=4.566\times {10}^{-4}Q(T_{\mathrm{ex}})\frac{1}{5}{\mathrm{e}}^{62.5/T_{\mathrm{ex}}}{L}_{[\mathrm{CI}]370\mathrm{\mu }\mathrm{m}}^{\prime },\end{array}$$(4)

where Q(T_ex) = 1 + 3e^−23.6/T_ex + 5e^−62.5/T_ex is the CI partition function, and $L^{\prime }=3.25\times {10}^7{S}_{\mathrm{line}}\mathrm{\Delta }\upsilon {\nu }_{\mathrm{obs}}^{-2}{(1+z)}^{-3}{D}_L^2$ is the luminosity expressed in units of K km s⁻¹ pc² (Solomon et al. 1992). With the ratio $R_{\text{CI}}=L_{[\text{CI}]370\mathrm{\mu }\text{m}}^{\prime }/L_{[\text{CI}]609\mathrm{\mu }\text{m}}^{\prime }$ we can calculate the excitation temperature T_ex = 38.8/ln(2.11/R_CI), assuming both lines are optically thin. Given our observational 3σ upper limit on [CI]_609 μm emission, we estimate T_ex > 37 K, and therefore M_CI < 1.3 × 10⁷ M_⊙. This upper limit of neutral carbon mass is in the range of estimations for other high-z quasars (Walter et al. 2011; Pensabene et al. 2021).

4.4. Implications for the [OIII]_{88 μ[sans]m} and [CI]_{609 μ[sans]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. 7) 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. 7). Following Papadopoulos & Greve (2004), we can also estimate the expected flux of [CI]_609 μm using the value of M(H₂) 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⁻¹, which is consistent with the 3σ observational upper limit (< 0.24 Jy km s⁻¹). Given these estimations, deeper observations with ALMA of [CI]_609 μm should achieve this sensitivity and successfully detect the emission line.

5. 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]_63 μm (one of the few band 10 detections in the literature to date), [NII]_122 μm, [OI]_145 μm, [CII]_158 μm, [NII]_205 μm and [CI]_370 μm, and we measure upper limits for [OIII]_88 μm and [CI]_609 μm.

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³ cm⁻³), 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 non-detections of [OIII]_88 μm and [CI]_609 μm emission lines, which may be critical in the observation of these lines in other high-z 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.

Acknowledgments

We thank the referee for their constructive comments and suggestions to improve the paper. R.F.A. thanks ESO for the hospitality during the period of the 2022/2023 studentship, and Anna Feltre, Luke Maud, and Dirk Petry for the advice during the ALMA data reduction and the CLOUDY modeling. R.F.A. acknowledges support from the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “First Call for H.F.R.I. Research Projects to support Faculty members and Researchers and the procurement of high-cost research equipment grant” (Project 1552 CIRCE). R.J.A. was supported by the ANID BASAL project FB210003 and by FONDECYT grant and 1231718. This research was carried out in part at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. M.A. acknowledges support from FONDECYT grant 1211951 and ANID BASAL project FB210003. A.P. acknowledges support from Fondazione Cariplo grant no. 2020-0902. G.J.S. acknowledges support from the US National Science Foundation through grants AAG 1910107 and AAG 1716229. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00576.S, ADS/JAO.ALMA#2015.1.00883.S, ADS/JAO.ALMA#2016.1.00668.S, ADS/JAO.ALMA#2017.1.00899.S, ADS/JAO.ALMA#2018.1.00119.S, ADS/JAO.ALMA#2021.1.00726.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.

References

Appendix A: ALMA data reduction

When using the CASA task tclean, a Briggs weighting was chosen, together with a robust parameter of 2.0 (close to natural weighting). Other options were also explored, shown in Fig. A.1 for the [CII] line. For example, we can see that extended emission is missed when using a robust parameter of r = 0.5 (between uniform and natural weighting). On the other hand, with a r = 2 and applying a tapering of the visibilities with a value of 0.3″ no extra emission is recovered, while the angular resolution is degraded.

Fig. A.1. [CII] line emission maps for W2246–0526 obtained using different CASA cleaning weightings (top, r = 0.5; bottom, r = 2.0 with a uvtaper = 0.3″).

Regarding the [OI]_145μm line in W2246–0526, a telluric atmospheric feature can potentially affect the observations. To ensure that the detected line profile was not affected by an increase in the noise toward the edge of the spectral window where the atmospheric line lies, we investigated the transmission at the Atacama Pathfinder Experiment (APEX) site¹ and confirmed that the center of the telluric feature is offset by at least 600 km s⁻¹ with respect to the redshifted center of the [OI]_145μm line. Hence, most of the line profile is not affected by the atmospheric feature.

Appendix B: CLOUDY modeling

B.1. Best-fit model ${\chi }_{\nu }^2$

To estimate the χ² of the different CLOUDY models we compare all the observed and modeled independent line ratios for W2246–0526 using [CII]_158μm as a reference. In Fig. B.1 we present the reduced χ² of the best-fit model presented in Sect. 3.4. The observed upper limits for [OIII]_88μm and [CI]_609μm are consistent with the best-fit model, and we do not observe any systematics in the residuals.

Fig. B.1. Comparison between the W2246–0526 observed and best-fit model line ratios used for the minimization. The reduced χ2 is shown in the top left. The bottom panel shows the differences between the two in units of standard deviations.

B.2. Fixed α_ox

Following the discussion in 4.2.1, we show in Fig. B.2 the posterior probability distributions for each of the CLOUDY parameters with a fixed α_ox = −1.4. The ionization parameter is unconstrained and extremely high.

Fig. B.2. Marginalized 1D posterior probability distributions for each parameter of the models without varying αox. The lines and shaded regions are as defined in Fig. 8.

Additionally, Fig. B.3 shows the same line ratios as in Fig. 3 for grids of CLOUDY models with α_ox = −1.4 and A_V = 10 mag. This showcases how lower X-ray to UV ratios and lower extinctions are not able to reproduce our observed line ratios.

Fig. B.3. Line ratio diagrams from CLOUDY AGN models for solar metallicity, αox = −1.4 with AV = 101 mag (labeled as A, grid with dashed lines), and αox = −0.8 with AV = 102.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.

B.3. Tests of additional parameters

Figure B.4 shows the same line ratios as in Fig. 3 for grids of CLOUDY models testing additional physical parameters. For all of 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. 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.

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.

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.

All Tables

All Figures

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 (H2). 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.

In the text

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.

In the text

	Fig. 2. continued.
In the text

Fig. 3. Line ratio diagrams from the CLOUDY models for the metallicity, αox and AV of the model with the lowest ${\chi }_{\nu }^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 ${\chi }_{\nu }^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.

In the text

Fig. 4. Line-to-FIR[42 − 122 μm] luminosity ratio as a function of FIR[42 − 122 μm] luminosity for [OI]63 μm, [OIII]88 μm, [NII]122 μm (top row), [OI]145 μm, [CII]158 μm, [NII]205 μm (middle row), [CI]370 μm, and [CI]609 μm (bottom row). W2246–0526 is indicated with a red square. The high-z quasars are taken from Weiß et al. (2003), Uzgil et al. (2016), Lu et al. (2017), Venemans et al. (2017), Walter et al. (2018), Lee et al. (2019), Novak et al. (2019), Li et al. (2020), Pensabene et al. (2021), Meyer et al. (2022). Upper limits for high-z quasars and W2246–0526 are shown with downward arrows. We note that the only high-z quasar shown for [OI]63 μm and [CI]609 μm corresponds to the Cloverleaf quasar (Weiß et al. 2003; Uzgil et al. 2016).

In the text

	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 Sect. 3.4 with log(U) = − 0.5 and αox = −0.8.
In the text

	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 1 Z⊙, the AV to 10 mag, and the hydrogen density to 103 cm−3. The gray dashed line (denoted by * in the legend) represents αox = −0.4 but with an AV = 100 mag.
In the text

	Fig. 7. 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 1024 cm−2.
In the text

Fig. 8. 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.

In the text

	Fig. 9. 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 fine-structure 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.
In the text

	Fig. A.1. [CII] line emission maps for W2246–0526 obtained using different CASA cleaning weightings (top, r = 0.5; bottom, r = 2.0 with a uvtaper = 0.3″).
In the text

	Fig. B.1. Comparison between the W2246–0526 observed and best-fit model line ratios used for the minimization. The reduced χ2 is shown in the top left. The bottom panel shows the differences between the two in units of standard deviations.
In the text

	Fig. B.2. Marginalized 1D posterior probability distributions for each parameter of the models without varying αox. The lines and shaded regions are as defined in Fig. 8.
In the text

Fig. B.3. Line ratio diagrams from CLOUDY AGN models for solar metallicity, αox = −1.4 with AV = 101 mag (labeled as A, grid with dashed lines), and αox = −0.8 with AV = 102.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.

In the text

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.

In the text