X-ray view of a massive node of the Cosmic Web at z=3 II. Discovery of extended X-ray emission around a hyperluminous QSO

TL;DR

This study reports the discovery of extended X-ray emission around a hyperluminous QSO at z=3, revealing hot gas in the circumgalactic medium with properties similar to those in galaxy clusters, providing insights into early hot halo formation.

Contribution

First detection of extended hot X-ray emitting gas around a high-redshift QSO, characterizing its properties and implications for galaxy formation models.

Findings

Hot gas temperature ~1.8 keV and mass ~2.6e12 solar masses.

X-ray luminosity exceeds typical values for similar systems.

Hot gas constitutes about 8.3% of the halo's virial mass.

Abstract

While the warm, ionized gas in the CGM at z>3 is now routinely observed around bright QSOs in Lya emission, little is known about the CGM hot phase due to its expected faintness in the X-ray band, often referred to as the ICM. Here, we report the analysis of 634 ks of Chandra X-ray observations in the MQN01 Cosmic Node, a region containing one of the brightest Lya nebulae and the largest galaxy overdensity discovered so far at z>3. We detect 66 net counts of X-ray emission in the 0.5-2 keV band extending to at least 30 kpc from the brightest QSO in MQN01. The morphology and spectrum are consistent with thermal emission from hot plasma in CIE. Photoionization is negligible, and IC is disfavored. A joint spatial and spectral MCMC analysis provides consistency with a beta-model with a steep density profile and a gas temperature kT~1.8 keV and virial halo mass Mvir~3e13 Mo. The inferred hot gas mass is Mhot(<Rvir)~2.6e12 Mo, which is ~8.3% of Mvir, or ~56% of the theoretical cosmological baryon budget of the halo. The hot gas also emits an exceptionally high Lx, with a measured L2-10~2.3e45 erg/s within the central 30 kpc. This system is a clear outlier in the Lx-Tx plane, indicating a thermodynamic state distinct from that of evolved lower-redshift hot halos. The cooling time in the inner 15-30 kpc is comparable to the local dynamical time, suggesting that the gas could become locally unstable in the absence of heating or feedback. Moreover, the thermal pressure associated with the detected CGM hot phase is large enough to confine the cold and dense clumps, which are required to reproduce the high Lya emission associated with the inner regions of the MQN01 structure. Although limited to a single system, our results provide unique information on the multi-phase properties of the CGM and a view of the nascent thermal hot gas phase observed in local galaxy clusters.

Tables

Table 1. Table 1: Median, 16th, and 84th percentiles of the posterior distributions from the MCMC analysis, along with the derived physical properties of the hot gas halo.

Parameter	Units	Value
$kT$	[keV]	${1.81}_{-0.37}^{+0.38}$
$\mathrm{norm}$	[${10}^{-4}{\text{cm}}^{-5}$]	${2.13}_{-0.82}^{+1.75}$
$Z^{†}$	$\left[Z_{\odot }\right]$	${0.014}_{-0.012}^{+0.083}$
$\beta$		${2.04}_{-0.75}^{+1.35}$
$r_{\text{core}}$	[kpc]	${36}_{-13}^{+16}$
$nor{m}_{\mathrm{pow}}$	[${10}^{-5}$]	${2.07}_{-0.18}^{+0.15}$
$\Gamma$		${2.09}_{-0.07}^{+0.06}$
$T$	[${10}^7$ K]	${2.1}_{-0.4}^{+0.4}$
$M_{\mathrm{vir}}$	[${10}^{13}{M}_{\odot }$]	$3\pm 1$
$R_{\mathrm{vir}}$	[kpc]	${190}_{-20}^{+19}$
$n_{\text{e},0}$	[${\text{cm}}^{-3}$]	${0.86}_{-0.19}^{+0.48}$
$M_{\mathrm{hot}\mathrm{gas}}\left(R_{\mathrm{vir}}\right)$	[${10}^{12}{M}_{\odot }$]	${2.6}_{-0.6}^{+1.7}$
$M_{\mathrm{hot}}/M_{\mathrm{vir}}$		${0.083}_{-0.030}^{+0.098}$
$f_{\mathrm{hot}}$		${0.56}_{-0.20}^{+0.65}$

Table 2. Table 2: Soft (0.5-2 keV) and hard (2-10 keV) absorbed fluxes (observed energy), and unabsorbed luminosities (rest-frame energy) of the thermal and power-law components in the central and intermediate annuli, including 1 σ \sigma uncertainties. These values are derived assuming the best-fit parameters reported in Table  1 .

Regions	Components	$F_{0.5-2\mathrm{keV}}$ / $F_{2-10\mathrm{keV}}$	$L_{0.5-2\mathrm{keV}}$ / $L_{2-10\mathrm{keV}}$
		[$\mathrm{erg}/{\mathrm{cm}}^2/s$]	[$\mathrm{erg}/s$]
		(${10}^{-15}$)	(${10}^{44}$)
0-2$″$	thermal	${5.38}_{-1.14}^{+1.64}$ / ${0.12}_{-0.05}^{+0.07}$	${12.22}_{-3.32}^{+4.68}$ / ${6.15}_{-1.29}^{+2.03}$
	power law	${44.11}_{-1.69}^{+1.74}$ / ${46.49}_{-1.21}^{+1.76}$	${51.63}_{-4.36}^{+4.61}$ / ${52.56}_{-1.97}^{+2.02}$
2$″$-3$″$	thermal	${2.59}_{-0.46}^{+0.53}$ / ${0.06}_{-0.02}^{+0.03}$	${5.79}_{-1.44}^{+2.11}$ / ${3.00}_{-0.55}^{+0.61}$
	power law	${0.79}_{-0.03}^{+0.03}$ / ${0.84}_{-0.02}^{+0.03}$	${0.93}_{-0.08}^{+0.08}$ / ${0.95}_{-0.04}^{+0.04}$
3$″$-5$″$	thermal	${1.82}_{-0.39}^{+0.40}$ / ${0.04}_{-0.01}^{+0.02}$	${4.03}_{-1.07}^{+1.64}$ / ${2.13}_{-0.45}^{+0.43}$
	power law	${0.31}_{-0.01}^{+0.01}$ / ${0.33}_{-0.01}^{+0.01}$	${0.36}_{-0.03}^{+0.03}$ / ${0.37}_{-0.01}^{+0.01}$

Equations

$$I_{\mathrm{anis}}=\frac{1}{N}\sum _{i=1}^N\left|\frac{C_i-\mu }{\mu }\right|,$$

$$\text{SB}(r)={\text{SB}}_0{\left[1+{\left(\frac{r}{r_{\text{core}}}\right)}^2\right]}^{-3\beta +1/2},$$

$$nor{m}_{r_1,r_2}=nor{m}_{2,3}\times \frac{{\int }_{r_1}^{r_2}\text{r SB(r)}dr}{{\int }_{2\text{arcsec}}^{3\text{arcsec}}\text{r SB(r)}dr}.$$

$$nor{m}_{r_{\mathrm{in}},r_{\mathrm{out}}}\text{leavevmode}[{\text{cm}}^{-5}]=\frac{10^{-14}{\int }_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}{n}_{\text{e}}(r)n_{\text{H}}(r)dV}{4\pi [D_A(1+z){]}^2},$$

$$n_{\text{e}}(r)=n_{\text{e},0}{\left[1+{\left(\frac{r}{r_{\text{core}}}\right)}^2\right]}^{-\frac{3\beta }{2}},$$

$$\int n_{\text{e}}^{2}(r)dV=\frac{4\pi r_{\text{core}}^{3}n_{\text{e},0}^{2}\,\tilde{\chi}(3\beta)}{\beta_{0}}\Biggl{[}\Biggl{(}1+\frac{r_{in}^{2}}{r_{\text{core}}^{2}}\Biggr{)}^{-\frac{\beta_{0}}{2}}-\Biggl{(}1+\frac{r_{out}^{2}}{r_{\text{core}}^{2}}\Biggr{)}^{-\frac{\beta_{0}}{2}}\Biggr{]}$$

$$\tilde{\chi }(a)={\int }_0^{+\infty }(1+x^2{)}^{-a}dx=\frac{\sqrt{\pi }}{2}\frac{\Gamma (a-\frac{1}{2})}{\Gamma (a)}.$$

$$n_{\text{e},0}=\sqrt{\frac{1.17\times 10^{14}[D_A(1+z){]}^{2}nor{m}_{r_{\mathrm{in}},r_{\mathrm{out}}}{\beta }_0}{r_{\text{core}}^3\tilde{\chi }(3\beta )[(1+r_{\mathrm{in}}^2/r_{\text{core}}^2{)}^{-\frac{{\beta }_0}{2}}-(1+r_{\mathrm{out}}^2/r_{\text{core}}^2{)}^{-\frac{{\beta }_0}{2}}]}}$$

$$M_{\text{hot gas}}=\frac{4\pi\,n_{\text{H},0}\,m_{p}}{X}\int_{0}^{R_{\text{vir}}}r^{2}\,\Biggl{[}1+\Biggl{(}\frac{r}{r_{\text{core}}}\Biggr{)}^{2}\Biggr{]}^{-3\beta/2}dr$$

$${\text{SB}}_{\text{Ann}}=ECF\times \frac{Net\text{leavevmode}Counts}{t_{\text{exp}}\times A_{\text{ann}}}\times \frac{ExpMa{p}_{\text{max}}}{ExpMa{p}_{\text{ann}}},$$

$$\tilde{\chi }(a)={\int }_0^{+\infty }(1+x^2{)}^{-a}dx=\frac{1}{2}{\int }_0^1{y}^{a-\frac{3}{2}}(1-y{)}^{-\frac{1}{2}}dy=\frac{1}{2}B\left(a-\frac{1}{2},\frac{1}{2}\right)=\frac{\sqrt{\pi }}{2}\frac{\Gamma \left(a-\frac{1}{2}\right)}{\Gamma (a)},$$

$$t_{\text{BV}}=\sqrt{\gamma \frac{h_S}{g}},$$

$$h_S(r)=\frac{r^2+r_c^2}{2\beta r}$$

$$t_{\text{BV}}(r)=\sqrt{\frac{5(r^2+r_c^2)}{6\beta rg(r)}}.$$

Full text

1 Introduction

Galaxy clusters host the largest reservoirs of hot gas known as the Intracluster Medium (ICM), which is prominently observed through X-ray emission. According to the hierarchical structure formation paradigm, these massive systems form through the gravitational collapse and merging of smaller structures (Peebles, 1980; Navarro et al., 1996). Gas infalling into the cluster potential well undergoes shock heating, reaching temperatures of $\sim 10^{7}-10^{8}\leavevmode\nobreak\ \rm K$ (Rosati et al., 2002; Kravtsov & Borgani, 2012). The dominant X-ray emission mechanism of the ICM is thermal bremsstrahlung, producing luminosities from $L_{X}\sim 10^{43}$ to more than $10^{45}\rm\leavevmode\nobreak\ erg\leavevmode\nobreak\ s^{-1}$.

At higher redshifts ($2<z<3$), massive halos are expected to host a nascent hot-Circumgalactic Medium (CGM) or proto-ICM, influenced by both gravitational processes and non-gravitational effects such as AGN and supernova feedback, although the formation scenario is still unclear (Shimakawa et al., 2018; Kooistra et al., 2022). Detecting extended X-ray emission from these structures is particularly challenging due to instrumental sensitivity limits and the need for high angular resolution to separate diffuse emission from point sources. This is crucial at $z>2$, where massive halos often host a central AGN that can dominate the X-ray signal. Models and simulations predict that the X-ray emission from hot gas at $z>2$ is faint and spatially compact, with an extent on the order of a few tens of arcseconds (i.e., $\sim$tens kpc; Saro et al., 2009). Therefore, observations require not only high sensitivity but also high angular resolution, a capability currently achieved only by the Chandra X-ray telescope (the mirror FWHM=0.2$\arcsec$ at 1.4 keV) and the ability to perform spatially resolved spectral analysis. Therefore, Chandra’s capabilities are unparalleled for studying extended X-ray emission on scales from kpc to a few Mpc, particularly in high-z systems. Constraining the thermal properties of the CGM/proto-ICM at these early epochs is essential for testing models of the redshift evolution of galaxy halos, such as the transition between cold and hot accretion modes (e.g., Dekel & Birnboim, 2006), and for understanding the role of feedback in shaping the thermodynamic state of diffuse gas in forming clusters.

The most distant potential X-ray detection of thermal emission from the ICM associated with a relaxed galaxy cluster has been reported in a system at $z\sim 2.5$ by Wang et al. (2016). A robust detection of thermal ICM emission was observed in the galaxy cluster XLSSC122 at $z\simeq 2$ by van Marrewijk et al. (2024). At redshifts $z>3$, most of the extended X-ray emission has been found around radio-loud galaxies (Yuan et al., 2003; Smail et al., 2012), and is thus associated with non-thermal emission from inverse Compton (IC) scattering of cosmic microwave background (CMB) photons by relativistic electrons in radio jets, as confirmed by the X-ray hard spectrum and the overlap with the radio jets. Recently, Tozzi et al. (2022) and Lepore et al. (2024) (hereafter L24) have reported direct evidence of proto-ICM at high redshift, presenting a detailed analysis of the extended thermal X-ray emission within 150 kpc of a radio galaxy in the Spiderweb protocluster at $z=2.156$. After accounting for non-thermal IC emission from the radio jet, they identified diffuse thermal emission from hot ($kT\sim 2$ keV) gas with a mass of $\sim 1.6\times 10^{12}\leavevmode\nobreak\ \rm M_{\odot}$, within a halo whose total mass is estimated to be $\sim 0.6-1.4\times 10^{13}\leavevmode\nobreak\ \rm M_{\odot}$. The combination of X-ray and Sunyaev-Zeldovich (SZ; Di Mascolo et al., 2023) detections revealed an entropy profile indicative of a cool-core system, with a cooling time of less than 100 Myr and a potential mass deposition rate of 250-1000 $\rm M_{\odot}\leavevmode\nobreak\ yr^{-1}$, consistent with infrared-based star formation rates of the central galaxy. These results suggest that AGN feedback and cooling flows may coexist in the early stages of protocluster formation.

This work is the second in a series of papers investigating deep (634 ks) Chandra X-ray observations of a protocluster at $z=3.25$, centered on the luminous quasi-stellar object (QSO) CTS G18.01 (hereafter ID1; see Travascio et al., 2024). This QSO lies at the center of the Multi Unit Spectroscopic Explorer (MUSE) Quasar Nebula 01 (MQN01), a giant ($>200\leavevmode\nobreak\ \rm kpc$) Ly$\alpha$  nebula in the sample of Borisova et al. (2016). Recent observations, based on a mosaic of eight VLT/MUSE AO-WFM pointings, have further extended its detection beyond two arcminutes (Cantalupo et al. in prep.). To investigate the connection between the properties of this cosmic web node and its galaxy population, multi-wavelength campaigns are actively characterizing the protocluster members (Pensabene et al., 2024; Galbiati et al., 2024). In our previous paper, we conducted a census of X-ray AGNs embedded in the protocluster, identifying six AGNs within an area of $\approx 16\rm\leavevmode\nobreak\ cMpc^{2}$ and a velocity range of $\pm 1000\leavevmode\nobreak\ \rm km\leavevmode\nobreak\ s^{-1}$ from the redshift of the central ID1. This corresponds to a signifcant overdensity of $log(L_{2-10\leavevmode\nobreak\ \rm keV/\rm erg\leavevmode\nobreak\ s^{-1}})>43$ X-ray AGN relative to the X-ray AGN space density in the field as constrained by Gilli et al. (2007).

In this second paper, based on the same Chandra X-ray dataset, we focus on the extended X-ray emission around the luminous QSO ID1. In Section 2, we describe the Chandra data reduction and astrometric correction, respectively. In Section 3.1, we report a significant detection of spatially resolved X-ray emission at observed energies below 2 keV, extending out to at least 30 kpc from the QSO. Section 3.2 examines the morphology of this emission in comparison with the Ly$\alpha$  nebula previously identified around the same QSO by Borisova et al. (2016), while Section 3.3 presents an analysis of its spectral properties. In Section 4, we investigate the physical properties of the hot gas under the assumption of a thermal origin. Discussions about the potential influence of QSO photoionization on the thermal emission and alternative emission mechanisms are presented in Sections 5.1 and 5.2. Finally, in Sections 5.4 and 5.6, we compare the properties of this system with those of the hot halo in the Spiderweb protocluster, with predictions from numerical simulations, and with observations of local galaxy clusters and groups.

Throughout this paper, all energies are reported in the observed frame unless stated otherwise, with fluxes consistently referring to observed energy bands. In contrast, luminosities are given in rest-frame energy bands. We adopt the same cosmological model as Travascio et al. (2024)111$H_{0}=67.4\rm\leavevmode\nobreak\ km\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ Mpc^{-1}$, $\Omega_{\Lambda}=0.714$, and $\Omega_{m}=0.286$, in which 1$\arcsec$ corresponds to 7.663 physical kpc. In Section 4, we assume a cosmic baryonic fraction of $f_{b}=\Omega_{b}/\Omega_{m}=0.15$ to evaluate the reliability of our best-fit model. Unless otherwise noted, all uncertainties shown in the plots correspond to the 1$\sigma$ (68$\%$) confidence level.

2 Methods

2.1 Chandra observations and data reduction

For this study, we used X-ray observations taken with the Chandra Advanced CCD Imaging Spectrometer (ACIS-I) during Cycle 23 (2022/2023), with a total exposure time of 634 ks (PI: S. Cantalupo; see Travascio et al., 2024, for further details). These observations were conducted in Very Faint (VFAINT) mode, which optimizes the separation of good and bad X-ray events by using the $5\times 5$ pixel event island for improved event classification. We performed the data reduction using CIAO 4.17 and the Chandra Calibration Database (CALDB 4.12.2) installed on Python 3.10. We removed flares identified in the light curves and corrected the astrometry of each ObsID according to the coordinates of the brightest QSO extracted in the GAIA catalog.

2.2 Refined astrometric alignment of ObsIDs

In Travascio et al. (2024), the different ObsIDs were aligned by cross-matching the positions of bright, point-like sources and estimating the least-squares minimization to determine the optimal transformation relative to the deepest ObsID. This is the best method to produce a catalog of X-ray sources in a wide field. However, in this paper, our main purpose is to explore the presence of extended X-ray emission around the AGNs identified in Travascio et al. (2024). This aim requires a more refined astrometric alignment, made ad hoc for the individual X-ray sources, minimizing positional uncertainties and ensuring that any observed extended emission is not an artifact of residual misalignment. However, this approach is feasible only for QSO ID1, which has enough photon counts in each ObsID to reliably determine the centroid position. For the other AGNs, multi-source alignment remains the most suitable option. More specifically, for QSO ID1, the 0.5-2 keV image of each ObsID was re-binned to one-fourth of the native pixel size and smoothed with a Gaussian kernel of 3 sub-pixels. The position of the QSO ID1 in each image was determined by computing the centroid’s coordinates within a 2$\arcsec$-radius circular region, centered on the approximate location of the QSO. Using the wcs_match tool, we generated a translation matrix for all ObsIDs, with the deepest ObsID used as the reference. We then applied this transformation to reproject the event and aspect files using the wcs_update tool. Figure 12 in Appendix A shows the results of this alignment, displaying the 0.5-2 keV images of the ObsIDs after re-binning and Gaussian smoothing. We then generated a new merged event file by combining the aligned ObsIDs.

3 Results

3.1 Detection and validation of extended X-ray emission around the hyperluminous QSO ID1

We investigated the presence of extended X-ray emission around the six X-ray-detected AGN in the MQN01 protocluster (Travascio et al., 2024). For this purpose, we compared the observed emission with the expected point spread function (PSF), considering different energy bands. At the position of each AGN and in each selected energy band, we simulated 500 individual PSFs predicted in each ObsID in our dataset. To generate these 500 simulated PSFs, we used the simulate_psf tool222To model the PSF, we tested blur parameter values ranging from 0 to the default value of 0.25. The blur parameter is meant to account for detector effects and uncertainties in aspect reconstruction that affect the extraction region. We found that variations in the blur parameter did not significantly impact our results, and we adopted a value of 0 for the final PSF modeling. in CIAO, which accounts for the instrument’s response and observational conditions associated with each ObsID. This method requires a spectral model for the energy band of interest. Therefore, we extracted a spectrum within a 2$\arcsec$ radius for each AGN, as this radius contains more than 98$\%$ of the PSF photons, independent of the spectral model. For each source and energy band, the composite PSF was constructed by weighting each simulated PSF by the exposure time of the corresponding ObsID and summing them.

After generating the PSF images, we compared their radial surface count profiles with those extracted from the observations. We estimated the radial surface count profile of the central AGN and subtracted the background level, determined from the outermost radial bins. For a consistent comparison with the PSF profile, we normalized the simulated PSF by scaling its total counts within a 2$\arcsec$ radius to match the corresponding counts in the observed, background-subtracted count profile centered on the QSO. To mitigate the uncertainties about the PSF centering, we aligned the PSF peak with the centroid of the observed X-ray emission and re-projected the PSF onto the same pixel grid using first-order interpolation. Finally, we quantified any residual excess emission in terms of $\sigma$ by subtracting the expected (PSF+background) counts from the observed counts in each radial bin and dividing by the combined uncertainty, which includes both statistical (Poisson) errors and background estimation uncertainties.

Our analysis reveals a significant extended X-ray emission associated with only one AGN in MQN01, the QSO ID1, in the soft X-ray regime (0.5-2 keV observed, corresponding to $\sim$2-8.5 keV rest-frame at $z=3.2502$). Figure 1 displays the radial profiles and residuals for ID1 in the 0.5-2.0 keV (left panel) and 2.0-10.0 keV (right panel) energy bands, extending out to $\sim$12$\arcsec$ ($\approx 90\leavevmode\nobreak\ \rm kpc$), by masking the nearest X-ray neighbors of the QSO, which are ID3 and ID4 in Travascio et al. (2024). The radial bins were chosen to ensure a minimum S/N of 3 in the data, except for the first bin, which has a fixed radius of 2$\arcsec$. Significant extended X-ray emission below 2 keV is observed at $\approx$15-30 kpc (2$\arcsec$-4$\arcsec$) from QSO ID1, with a significance of $\sim$2.5$\sigma$ per radial bin. For comparison, Figure 13 in Appendix B presents radial profiles and residuals for AGNs ID2, ID5, and ID6 in the 0.5-2.0 keV and 2.0-10.0 keV bands, where residuals remain below 2$\sigma$, consistent with PSF predictions. AGN ID3 and ID4 were excluded from this analysis due to blending, which prevents a reliable assessment of the presence of extended emission.

Figure 2 shows the 0.5-2 keV X-ray images of QSO ID1: the observed data (left), the simulated and rescaled PSF (middle), and the PSF-subtracted map (right). Red dashed circles indicate the radial bins used for the extraction of the count profile in Figure 1. The PSF was scaled such that its integrated counts within the central 2$\arcsec$ matched those of the observed data. The rightmost panel reveals the spatial distribution of residual counts beyond 2$\arcsec$. A similar image for AGN ID2 is shown in Figure 14 (Appendix B).

The detection of extended X-ray emission around the QSO ID1 was confirmed using both the standard merged event file, aligned with multiple reference sources (Travascio et al. 2024), and the newly merged event file obtained via the refined alignment procedure described in Section 2.2 and Appendix A. Both methods yielded consistent results, though the latter provided a more conservative detection, reducing excess counts within the 2$\arcsec$-4$\arcsec$ annulus by $\lesssim 10\%$. To assess potential PSF estimation errors due to spatial dependencies of the PSF and possible misalignments, we repeated the calculation assuming PSF models estimated from different locations in the field of view. We found that only PSF models generated at positions offset by more than $\sim$12$\arcsec$ from the QSO location produced a noticeable impact on the radial profiles. This is far larger than any realistic astrometric misalignment in our data (see, e.g., Figure 12).

Additionally, we tested an alternative approach in which the PSF was directly subtracted from each ObsID’s 0.5-2 keV image before merging, eliminating any potential biases from exposure-weighted PSF stacking. The final PSF-subtracted image remained unchanged, further confirming the robustness of our detection and PSF modeling.

After subtracting the background and the PSF from the QSO image in the 0.5-2 keV band $\sim 66\pm 8$ counts remain, indicating a detection with a significance of $\sim 8\sigma$. This estimate should be considered a lower limit on the extended emission, as it assumes no diffuse component within the central 2$\arcsec$ (i.e., $\sim$15 kpc).

3.2 Morphology of the extended X-ray emission

Around the QSO ID1, there is also evidence of warm ($T\sim 5\times 10^{4}\leavevmode\nobreak\ \rm K$) line-emitting gas in the CGM, as indicated by the Ly$\alpha$  nebula (i.e., MQN01) identified by Borisova et al. (2016), which extends over $200\leavevmode\nobreak\ \rm kpc$. If the extended X-ray emission traces hot ($T>10^{6}\leavevmode\nobreak\ \rm K$) halo gas, this system provides a unique high-redshift case for multi-phase CGM studies, particularly noteworthy given the absence of a detected radio jet, which rules out jet-driven excitation of the Ly$\alpha$  emission. This provides a remarkable example of the potential coexistence of warm and hot gas phases.In this section, we focus on comparing the morphology of the extended Ly$\alpha$  and X-ray emission. Future work will involve a detailed analysis of their physical properties, exploring the implications for gas mass, phase mass fractions, and the underlying processes involved.

Figure 3 shows the smoothed map of the 0.5-2 keV counts, obtained after subtracting the QSO’s PSF contribution as outlined above. The PSF was rescaled to match the sum of the counts within 2$\arcsec$ from the QSO’s center, which is marked with a black transparent area, while the black filled dot is centered on the QSO with a 1$\arcsec$ radius. The map was smoothed using the aconvolve tool in CIAO-4.15 with a Gaussian kernel of size $2\times 2$ pixels, a normalization factor of 1, and a standard deviation of 1 pixel. The magenta contours trace the Ly$\alpha$  nebula at surface brightness (SB) levels of 2.5, 4, 6, 8, 10 and $12\times 10^{-18}\rm\leavevmode\nobreak\ erg\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ cm^{-2}\leavevmode\nobreak\ arcsec^{-2}$, providing a direct visualization of the nebular extent as related to the X-ray extended emission morphology.

We observe that the majority of the soft X-ray counts are confined within the Ly$\alpha$  SB level of $\approx 6\times 10^{-18}\rm\leavevmode\nobreak\ erg\leavevmode\nobreak\ s^{-1}\leavevmode\nobreak\ cm^{-2}\leavevmode\nobreak\ arcsec^{-2}$, and seem to be distributed isotropically, with no preferential directions. The extended X-ray morphology appears to follow the south tail of the Ly$\alpha$  contours, while the east Ly$\alpha$  extension is not connected to any specific X-ray feature. The companion of the QSO ObjB, reported in Pensabene et al. (2024), is located approximately $1\arcsec$ from the QSO (red cross in Figure 3), well below the extension of the X-ray emission. Moreover, the observations do not indicate the presence of blended X-ray sources.

To investigate the spatial correspondence between the extended X-ray emission and the Ly$\alpha$  nebula, within a radial range of 2$\arcsec$-5$\arcsec$, we divided the region into eight sectors and calculated the mean flux in each azimuthal sector, effectively tracing the angular distribution of both components across the selected annulus. We adopt eight sectors to balance spatial resolution and statistical robustness, allowing us to detect potential azimuthal variations in both emissions. The resulting azimuthal profile, shown in Figure 4, shows a comparison of the trends in the normalized fluxes of the X-ray (black) and Ly$\alpha$  (red) emissions. While the X-ray emission does not show statistically significant azimuthal variations, we note that both X-ray and Ly$\alpha$  profiles exhibit a mild flux dip near an azimuthal angle of $\sim 360^{\circ}$, hinting at a possible, but not statistically robust, spatial alignment. This alignment hints at partial morphological coupling between the hot and cold CGM phases, although the X-ray emission remains largely isotropic. This relative isotropy of the hot gas component is consistent with expectations for a virialized or quasi-virialized halo, whereas the more anisotropic morphology of the Ly$\alpha$  emission likely reflects the clumpy and filamentary structure of the cooler, denser gas phases. To quantitatively assess the isotropy of the extended soft X-ray emission, we estimated an anisotropy index ($I_{\mathrm{anis}}$) in $N$ angular sectors. This index quantifies the average fractional deviation of the photon counts in each sector from a perfectly uniform distribution and is defined as:

[TABLE]

where $C_{i}$ is the number of counts in the $i$-th sector, $\mu$ is the mean count rate across all sectors, and $N$ is the total number of sectors. To assess deviations from uniformity in the extended emission, we applied both a $\chi^{2}$ test and a Kolmogorov-Smirnov (KS) test to the photon counts distributed in angular sectors. In all cases (including different energy bands and sector choices), we found low isotropy indices ($I_{\mathrm{anis}}<0.4$) and high $p$-values ($p>0.2$), with the KS test yielding a maximum $D$-statistic of 0.093 and $p$-value$\leavevmode\nobreak\ =1$. These results consistently indicate there is no statistically significant deviation from isotropy, suggesting that the emission is consistent with being isotropic within current statistical uncertainties.

3.3 Spectral Analysis of the Extended Emission

As a complementary way to investigate and characterize the presence of extended soft X-ray emission, we compare in Figure 5 two spectra: one extracted from the inner 2$\arcsec$ region, where the emission is dominated by the QSO, and one from the 2$\arcsec$-4$\arcsec$ annular region (corresponding to $\sim$15-30 kpc), where an excess of soft X-ray emission was detected above $2\sigma$ in Figure 1. The nuclear spectrum is of high quality, with a count rate of $4.15\times 10^{-3}\leavevmode\nobreak\ \rm cts\leavevmode\nobreak\ s^{-1}$ and a total of 2633 counts, allowing robust spectral fitting. The annular region contains a total of 172 counts in the 0.5-10 keV band and 96 counts in the 0.5-2.0 keV band.

We initially fitted the nuclear spectrum, binning it to a minimum of 50 counts per bin, using the C-statistic (C-stat) and a simple power-law model with Galactic absorption ($N_{\rm H,Gal}=1.15\times 10^{20}\leavevmode\nobreak\ \rm cm^{-2}$; Travascio et al., 2024). The analysis was performed with Sherpa (Freeman et al., 2001; Siemiginowska et al., 2024) and cross-checked with XSPEC (Arnaud, 1996). The best-fit model without intrinsic absorption yields a C-statistic of 41.19 for 41 degrees of freedom, with a photon index $\Gamma=2.3\pm 0.1$ and a normalization of $(2.74\pm 0.36)\times 10^{-5}$. The corresponding $\chi^{2}$ is 38.50 for 43 bins (null hypothesis probability $p=0.582$). We also tested the inclusion of intrinsic absorption using the xszphabs model. The fit improves slightly to a C-statistic of 39.75 for 40 degrees of freedom, with a best-fit intrinsic column density of $N_{\rm H,int}<6.25\times 10^{22}\leavevmode\nobreak\ \rm cm^{-2}$. The lower bound of $N{\rm H,int}$ is consistent with zero, and the parameter is fixed at the hard lower limit, indicating that the fit does not require intrinsic absorption. The photon index and normalization remain well determined. The small $\Delta$C-stat of 1.44 for 1 additional degree of freedom confirms that the inclusion of intrinsic absorption does not significantly improve the fit. Therefore, we adopt the simpler power-law model with Galactic absorption only. No evidence of a 6.4 keV iron emission line was found, consistent with at most a negligible contribution from X-ray reflection off cold material in the accretion disk or torus.

Then, we estimated the spillover of nuclear emission into the 2$\arcsec$-4$\arcsec$ annulus. To do it, the best-fit nuclear model was rescaled using energy-dependent correction factors derived from simulated PSFs. These factors, shown in the top panel of Figure 5, represent the ratio of expected PSF counts in the annular region and the central 2$\arcsec$ aperture. The bottom panel of Figure 5 shows the nuclear spectrum (black points), its best-fit model (solid line), the annular spectrum (red points), and the rescaled nuclear spillover component (dashed line). The annular spectrum was binned to a minimum of 10 counts per bin. The comparison shows significant residuals in the annular spectrum after accounting for PSF contribution, particularly below 2 keV, with deviations reaching approximately the 2-4$\sigma$ level in individual bins. This setup (the AGN model in particular) allows us to perform the spectral analysis of the soft emission.

4 Thermal model of the extended emission

4.1 Spectral and spatial analysis of the extended thermal emission

In this section, we investigate the physical origin of the spectral excess observed in the X-ray emission (see Section 3.3) and its implications for the hot halo in MQN01. The relatively soft and isotropic spatially extended nature of the emission suggests that it may arise from thermal gas heated to X-ray-emitting temperatures by gravitational shocks (or feedback processes) within the QSO halo. Under this assumption, we model the spectral excess using the xsmekal model in Sherpa, which describes the emission from hot, optically thin plasma in collisional ionization equilibrium (CIE). This model incorporates atomic data for thermal bremsstrahlung, recombination, and line emission from highly ionized species (hereafter referred to as thermal CIE emission; see Mewe et al., 1985, 1986; Liedahl et al., 1995). In Section 5.2, we also show that relaxing the assumption of CIE by including the effect of quasar photo-ionization does not affect our results in any way.

We performed a joint spectral and spatial analysis of the extended X-ray emission using Markov Chain Monte Carlo (MCMC) methods (e.g., Ruppin et al., 2021), implemented via the Python package emcee (Foreman-Mackey et al., 2013). For the comparison between the models and the data, we used four spectra extracted from concentric regions: the central $2\arcsec$, and the annuli spanning $2\arcsec$-$3\arcsec$, $3\arcsec$-$5\arcsec$, and $5\arcsec$-$8\arcsec$. The central spectrum was binned to a minimum of 50 counts per bin, while the others were binned to a minimum of 5 counts per bin. These were modeled using a xsmekal model for the extended emission and a power-law component for the nuclear one (i.e., the PSF).

As input parameters of the MCMC fitting, we allow the plasma temperature ($kT$), normalization ($norm$), and metallicity ($Z/Z_{\odot}$) to vary freely, adopting solar elemental abundances from the Abundanc table in Asplund et al. (2009). We also vary the normalization ($\texttt{norm}_{\rm pow}$) and photon index ($\Gamma$) of the power-law component. We did not include intrinsic absorption for this nuclear component, based on the results presented in Section 3.3, where the column density $N_{H}$ is consistent with zero within $1.5\sigma$ when fitting the AGN spectrum alone. Nevertheless, even when we include an extended thermal component and allow $N_{H}$ to be a free parameter, the resulting best‑fit values remain consistent with the current ones within 1$\sigma$, although $N_{H}$ itself remains close to the lowest detectable limit in our data.

The temperature, metallicity, and normalization ($norm_{2,3}$) of the xsmekal model are those used to fit the $2\arcsec$-$3\arcsec$ annular spectrum (red spectrum in Figure 5). While we assume that the metallicity and temperature are constant across all the annuli, the parameter $norm_{2,3}$ is proportional to the thermal flux from the $2\arcsec$-$3\arcsec$ annulus. It is physically related to the average squared gas density integrated along the line of sight (see Equation 4 below). To model the spatial distribution of the thermal gas, we assumed a classical $\beta$-model profile (Cavaliere & Fusco-Femiano, 1976), commonly used to describe the ICM in low-redshift galaxy clusters (Mohr et al., 1999; Dong et al., 2010; Conte et al., 2011; Paggi et al., 2021). This model is defined as:

[TABLE]

where $\beta$ and the core radius $r_{\text{core}}$ are free parameters. The thermal normalizations for each spectrum were computed by scaling the initial $norm_{2,3}$ parameter as follows:

[TABLE]

Furthermore, we assumed a constant temperature and metallicity at each radius. The $norm_{\rm pow}$ parameter was scaled using energy-dependent simulated PSF profiles (upper panel of Figure 5).

We adopted an exponential prior for $\beta$ ($P(\beta)\propto e^{-\beta/\beta_{*}}$, with $\beta_{*}=1$), motivated by results from local ICM SB profiles, which generally favor values of $\beta$ below or close to unity, with rare exceptions as large as $\sim$4 (see Mohr et al., 1999; Dong et al., 2010; Conte et al., 2011; Paggi et al., 2021; Xue & Wu, 2000; Wise et al., 2004; Mirakhor et al., 2022). A log-uniform (i.e., scale-invariant) prior was assumed for $r_{\text{core}}$, $Z$, and for the normalization parameters ($norm_{2,3}$ and $norm_{\rm pow}$) to allow for efficient exploration of the parameter space.

To compare the model predictions with the observed counts, we used the Poisson likelihood (aka C-statistics). To fit each spectrum, we include the background, consisting of a power-law continuum and five Gaussian emission lines, obtained by modeling a spectrum extracted from a large off-source region near the QSO. We used 100 walkers and 10000 steps per walker, resulting in a total of $10^{6}$ samples for the posterior distribution.

Figure 6 shows the marginalized 1D and 2D posterior probability distributions for the model parameters derived from the MCMC. The contours represent the 68.3%, 95.5%, and 99.73% confidence levels. In the histograms, the red solid line indicates the median value of each parameter distribution, while the red dashed lines mark the 16th and 84th percentiles. The median values are also shown in the 2D projections with a red square. As expected, we found a partial degeneracy between the parameters of the xsmekal model ($kT$ and $norm_{2,3}$), between the parameters of the nuclear power-law ($\Gamma$ and $norm_{\rm pow}$) and between those of the beta model ($\beta$, $r_{\text{core}}$). The latter is the most marked, reflecting the absence of deep data at large radii, which would be needed to better constrain the spatial distribution of the gas. The degeneracy between $kT$ and $norm_{2,3}$ of the xsmekal model is mostly due to the lack of data at lower energies (at this redshift, we can see with Chandra only the exponential tail of the Bremsstrahlung). Despite these limitations, the marginalized 1D posterior probability distributions are all well-behaved and unimodal for all parameters, allowing us to put clear quantitative constraints on the properties of the hot gas. In particular, we found as median temperature and normalization of the thermal component: $kT=(1.8\pm 0.4)\leavevmode\nobreak\ \rm keV$ (corresponding to $T=(2.1\pm 0.4)\times 10^{7}\leavevmode\nobreak\ \rm K$), and $norm_{2,3}=2.13_{-0.82}^{+1.75}\times 10^{-4}\leavevmode\nobreak\ \rm cm^{-5}$. The nuclear power-law parameters are well constrained, with $norm_{\text{pow}}=(2.1\pm 0.2)\times\leavevmode\nobreak\ 10^{-5}$ and $\Gamma=2.1\pm 0.1$. Note that the best fit slope of the nuclear component is slightly shallower than (albeit within $2\sigma$ from) that inferred from the analysis of the nuclear spectrum alone (2.3; Section 3.3), reflecting a small but not entirely negligible contribution of the thermal emission also at small radii (see also Table 2 below). The $\beta$-model parameters describing the spatial distribution of the hot gas are $\beta=2.04_{-0.75}^{+1.35}$ and $r_{\text{core}}=36_{-13}^{+16}\leavevmode\nobreak\ \rm kpc$. A complete list of the best-fit parameters is provided at the top of Table 1. The only parameter that is loosely constrained by our fit is the metallicity. The posterior probability distribution for this parameter is very broad, reflecting an insufficient S/N for a precise measurement. The posterior probability nonetheless shows a clear preference for relatively low metallicities ($Z\lesssim 0.1Z_{\odot}$), mainly as a consequence of the non-detection of the Fe K$\alpha$ complex from He-like and H-like ions at rest-frame 6.7-6.9 keV (observed $\approx$1.6 keV). We stress, however, that metallicity is included in our MCMC fit only for marginalization purposes, and the formal best-fit value reported in Table 1 should be taken with caution.

Figure 7 presents the best-fit models (red lines) for the four spectra extracted from the regions defined previously: the central 2$\arcsec$ aperture (top left), and the three concentric annuli 2$\arcsec$-3$\arcsec$, 3$\arcsec$-5$\arcsec$, and 5$\arcsec$-8$\arcsec$, shown in the top-right, bottom-left, and bottom-right panels, respectively. The individual spectral components are plotted separately: the nuclear power-law (blue), background (purple), and thermal (green) components. The residuals in the bottom sub-panels of each spectrum are expressed in terms of $\sigma$. The best-fit models reproduce the observed spectra within $\sim 2\sigma$ across all energy bins. However, in the 2$\arcsec$-3$\arcsec$ and 3$\arcsec$-5$\arcsec$ spectra, the models tend to underpredict the emission in the lowest-energy bin (below 0.7 keV). However, emission below 1 keV may be subject to systematic uncertainties related to calibration accuracy, as discussed in Section 4.3.

4.2 Physical properties of the hot halo gas

We used the median posterior values to estimate several physical properties of the hot gas. Assuming the system is virialized, we derived a virial mass (equivalent to $M_{200}$) of $M_{\text{vir}}\approx 3\times 10^{13}\leavevmode\nobreak\ \rm M_{\odot}$ from the gas temperature and the redshift of the system, following Equation (A10) from Dekel & Birnboim (2006). This estimate is consistent with scaling relations reported in the literature (e.g., Sato et al., 2000; Wang & Abel, 2008; Ehlert & Ulmer, 2009). The corresponding virial radius (i.e., $R_{200}$), calculated assuming the same overdensity criterion, is $R_{\text{vir}}\simeq 190\leavevmode\nobreak\ \rm kpc$, based on Equation (2) from Ehlert & Ulmer (2009). This yields a ratio $r_{\text{core}}/R_{\text{vir}}\approx 0.2$, which falls within the typical range of $10^{-2}$ to $0.6$ observed for the ICM in local galaxy clusters and groups (e.g., Vikhlinin et al., 2005).

The electron and hydrogen density distribution in the plasma is related to the best-fit $norm_{2,3}$ parameter through the expression333See https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/node198.html:

[TABLE]

where $D_{A}$ is the angular diameter distance to the source, and the volume integral includes only the region along the line of sight that contributes to the $r_{\rm in}-r_{\rm out}$ annulus in projection. Accounting for the contribution of electrons from helium (0.17 per hydrogen ion, assuming full ionization and a helium mass fraction $Y=0.25$), we adopt $n_{\text{e}}(r)\simeq 1.17\,n_{\text{H}}(r)$, making the $norm$ parameter proportional to $\int n_{\text{e}}^{2}(r)\,dV$. This expression depends on the electron density profile, which is inferred by deprojecting the observed SB profile under the assumption that the X-ray emissivity scales with $n_{\text{e}}^{2}(r)$. For the adopted $\beta$-model (equation 2), the deprojected electron density distribution is

[TABLE]

where $n_{\text{e},0}$ is the central electron density. To relate this 3D density profile to the observed 2D spectral normalization in a given annulus, we perform a proper deprojection of the spherical volume intersected by the line of sight through the annulus. Specifically, we compute the volume integral in Equation (4) using cylindrical coordinates ($R$,$\phi$,$z$), where the volume element is $dV=RdR\,d\phi\,dz$ and $z=\sqrt{r^{2}-R^{2}}$. The azimuthal integration yields a factor $2\pi$, using the analytical form of the line-of-sight projection for spherically symmetric functions, commonly known from the Abel transform framework, we obtained:

[TABLE]

where $\beta_{0}=6\beta-3$, while the function $\tilde{\chi}$ is defined as:

[TABLE]

A proof of the last equality in equation (7) is provided in Appendix C. Note that $\tilde{\chi}$ is related to the function $\chi$ defined in Pezzulli et al. (2017) by the relation $\tilde{\chi}(a)=\chi(2a)$.

This formulation (6) accurately captures the projection of the spherical density distribution into the annular region observed in the sky. By inserting this expression into Equation (4), we derived a direct relationship between the central density and the observed spectral normalization in the $n_{\text{e},0}$ corresponding $r_{\rm in}$-$r_{\rm out}$ annulus:

[TABLE]

where $D_{A}$ and $r_{\text{core}}$ are expressed in cm, $norm_{r_{\rm in},r_{\rm out}}$ in $\rm cm^{-5}$ and $n_{\text{e},0}$ in $\rm cm^{-3}$. Applying this method, we estimated a central electron density of $n_{\text{e},0}=0.9_{-0.2}^{+0.4}\leavevmode\nobreak\ \rm cm^{-3}$. We then computed the total mass of hot gas within the virial radius using:

[TABLE]

where $n_{\text{H},0}=n_{\text{e},0}/1.17$, $m_{p}=8.4\times 10^{-58}\leavevmode\nobreak\ M_{\odot}$, and $X=0.75$ the hydrogen mass fraction. This yields a hot gas mass of $M_{\text{hot gas}}(<R_{\text{vir}})=2.6_{-0.6}^{+1.7}\times 10^{12}\leavevmode\nobreak\ \rm M_{\odot}$ corresponding to a hot baryon fraction of $M_{\rm hot}/M_{\rm vir}\approx 0.083^{+0.098}_{-0.030}$. Assuming a cosmic baryon fraction $f_{b}=\Omega_{b}/\Omega_{m}=0.15$ (e.g. Mohr et al., 1999; Gonzalez et al., 2007; Ettori et al., 2009; Hinshaw et al., 2013; Ge et al., 2018) this implies that the hot gas accounts for $f_{\rm hot\leavevmode\nobreak\ gas}=(M_{\rm hot}/M_{\rm vir})/0.15=0.56_{-0.20}^{+0.65}$ of the baryons primordially associated with the halo. This fraction highlights that a substantial part of the halo’s theoretical baryon budget is found in the hot CGM phase. Here, by baryon budget, we refer to the total baryonic mass expected from the cosmic baryon fraction to be associated with the halo’s dark matter, independent of whether these baryons are currently located inside or outside the virial radius. This fraction is consistent with observations of low-mass groups and intermediate-mass clusters in the local Universe, where the hot gas typically accounts for about 30-85$\%$ of the halo’s baryonic content (e.g., Eckert et al., 2013; Morandi et al., 2015; Pratt et al., 2023). The fact that a similarly large mass of hot CGM is found already at $z\sim 3$, co-existent with a luminous Giant Ly$\alpha$ nebula, is in line with the predictions of Pezzulli & Cantalupo (2019), who explored the role of hot virialized gas in boosting the Ly$\alpha$ emissivity in MUSE QSO nebulae (including MQN01) through the compression of colder (Ly$\alpha$-emitting) gas. This is further discussed in Section 5.5.

Table 2 summarizes the soft (0.5-2 keV) and hard (2-10 keV) absorbed fluxes and unabsorbed luminosities of the thermal and nuclear components for each spectrum, assuming best fit (median) parameters, excluding the outermost annulus, where the emission lies below the background level. The results indicate that, within the central 2$\arcsec$ region, the observed thermal soft X-ray flux may account for approximately $12\pm 4\,\%$ of the total emission in the central $2\arcsec$. However, this estimate depends sensitively on the $\beta$-model parameters.

4.3 Systematic uncertainties on the X-ray luminosities

An important aspect to consider is the presence of systematic uncertainties affecting the measurement of the soft X-ray luminosities (Table 2). The detection of extended thermal emission in the $0.5$-$2$ keV rest-frame band may be susceptible to the shape of the spectral continuum at low energies, where the effective area of Chandra drops significantly. Since the thermal component is mainly constrained by photons below $\sim 2$ keV, any small number of counts or calibration inaccuracy in this range, although unlikely to fully mimic the observed signal, can steepen the spectral fit and bias the extrapolated luminosity estimate. However, it is worth noting that our main MCMC model already underpredicts the low-energy counts (see residual for the lowest energy bin in Figure 7), suggesting that any such systematic would only bring out the model in better agreement with the data. More generally, the mere detection of soft X-ray emission at such high redshift, where the observed energy range samples only the high-energy tail of the thermal bremsstrahlung spectrum, already implies that the emission must be intrinsically powerful.

Nonetheless, to mitigate any unmodelled systematics related to the lowest energy bin, we repeated the MCMC analysis using only the spectral data above $1$ keV (observed frame), where the Chandra effective area is more stable and remains above $60\leavevmode\nobreak\ \mathrm{cm}^{2}$. In this alternative fit, we applied a Gaussian prior on the temperature, centered at $kT=2.5$ keV (i.e., 64% higher than the previous best-fit value) with a standard deviation of $0.5$ keV. Imposing a prior with a higher temperature than our fiducial model is a conservative choice, as a higher temperature implies a shallower spectral slope, thereby reducing the extrapolated soft-band luminosity. Under these assumptions, the posterior distribution for $kT$ peaks around $\approx 1.9$ keV, with the corresponding total luminosity, within 30 kpc, reaching a minimum of $\approx 10^{45}\leavevmode\nobreak\ \mathrm{erg\,s^{-1}}$ within $1\sigma$. This test shows that, even when conservatively accounting for potential calibration systematics and introducing an explicit bias in favour of high temperatures, relatively low temperatures are still preferred and the inferred intrinsic soft-band luminosity remains significantly high.

5 Discussion

5.1 Impact of the QSO radiation on the thermal emission

Throughout our X-ray spectral analysis, we have implicitly assumed that the hot gas is in CIE, meaning that its ionization state and emitted spectrum are determined just by collisions with thermal electrons. However, this assumption may break down due to the presence of intense photoionizing radiation from the central QSO. In our case, the QSO hosted in the center of MQN01 protocluster is hyper-luminous, emitting copiously in both the UV and X-ray bands. This extreme radiation field has the potential to alter the ionization balance and emissivity of the surrounding gas, particularly within the central 30 kpc. It is therefore important to test whether AGN photoionization can significantly affect the thermal X-ray emission, and hence bias our interpretation of the physical properties derived under the CIE assumption.

We used Cloudy (Ferland et al., 1998) to model the emission from hot CGM gas, both with and without the influence of AGN photoionization. Specifically, we simulated a coronal gas as in our fiducial model. For this test, we focused on a shell with radii 15-23 kpc, which corresponds in projection to the $2\arcsec$-$3\arcsec$ annulus where the signal of the extended emission is the strongest (see Figure 7), and we assume a temperature and metallicity of $kT=1.8\leavevmode\nobreak\ \rm keV$ and $Z\simeq 0.014\,Z_{\odot}$, and a gas density $\langle n_{\text{H}}\rangle\simeq 0.4\leavevmode\nobreak\ \rm cm^{-3}$, which is the average density value in such shell for our fiducial beta model (see Section 5.4). The gas was placed at this distance in a shell geometry within Cloudy.

The energy injected by the central AGN was described using a spectral energy distribution (SED; panel (a) of Figure 8) tabulated from observations (”AGN T=3e5 K a(ox)=-1 a(uv)=-2 a(x)=-2.2”), constrained to match UV ($\nu L_{\nu}(400\,\text{\AA })=5\times 10^{46}\leavevmode\nobreak\ \rm erg\,s^{-1}$; Borisova et al., 2016) and X-ray ($\nu L_{\nu}(2\,\text{keV})=3.3\times 10^{45}\leavevmode\nobreak\ \rm erg\,s^{-1}$; Travascio et al., 2024) photometric measurements, adopting a X-ray power law photon-index of $\Gamma\approx 2.1$. This setup allows us to test the potential impact of AGN photoionization on the thermal emission from the surrounding gas. We also ran an additional Cloudy simulation where the QSO SED was suppressed by 6 orders of magnitude, to simulate a condition where photoionization is negligible and CIE applies. This was done to allow for a direct comparison with the xsmekal model (which also assumes CIE) and test the dependence of the results on different codes at fixed physical assumptions.

We then reproduced the corresponding thermal emission using the xsmekal function in sherpa, adopting the same parameters as before: $kT=1.8\leavevmode\nobreak\ \rm keV$, $Z=0.014\leavevmode\nobreak\ Z_{\odot}$, and a normalization, $norm_{2,3}$, related to the gas density $n_{\text{H}}=0.4\leavevmode\nobreak\ \rm cm^{-3}$ via Equation (4), and adopting, in this case, a volume equal to the volume of the shell. The resulting model was rescaled to express the emission in units of $\nu L_{\nu}$.

The panel (b) in Figure 8 presents the spectra obtained for the three models described above. The Cloudy thermal models with and without QSO photoionization are shown as blue and cyan lines, respectively, while the red line corresponds to the thermal emission from the xsmekal model. The gray transparent band highlights the energy range (0.5-2 keV observed) where extended X-ray emission is detected in the Chandra data, and where we aim to identify differences between the models. The panel (c) shows the ratio between the Cloudy and xsmekal thermal models without QSO photoionization, while the bottom panel displays the ratio between the Cloudy thermal models without and with the QSO radiation field. As shown by the similarity between the blue and cyan curves, whose ratio is shown in panel (d), the impact of QSO photoionization on the thermal emission of the hot gas is negligible in this case. Additionally, the comparison between the cyan and red curves, whose ratio is shown in panel (c), shows that, even under identical physical assumptions, the predicted emission varies only marginally between different modeling codes. These small differences are negligible for our purposes and do not affect our conclusion that QSO photoionization has no significant impact on the thermal emission of the hot gas.

We also tested how the presence of the QSO affects the gas heating-cooling balance, and thus its cooling time, using time-dependent Cloudy simulations. We enabled non-equilibrium cooling (set dynamics relax 3) and tracked the thermal evolution over 300 iterations. A series of time steps, ranging from 10 to 100 Myr, was defined (using time commands), starting with an initial timestep of $10^{4}$ years. The simulation stops once the gas cools below $10^{5}$ K, allowing us to monitor and compare the cooling progression with and without QSO irradiation. We found that the cooling time to reach temperatures below $10^{5}\leavevmode\nobreak\ \rm K$ was $\sim 80$ Myr when the effect of QSO radiation on the thermal balance was neglected, and increased to $\sim 90$ Myr when this contribution was included. This indicated that QSO photoionization could slow down the cooling rate (e.g., Cantalupo, 2010; Gnedin & Hollon, 2012), although at temperatures $kT\lesssim 2\leavevmode\nobreak\ \rm keV$ the effect remained negligible, provided that the metallicity was low ($Z\lesssim 0.1\leavevmode\nobreak\ Z_{\odot}$). We emphasize that these results do not imply that the gas will necessarily cool within 90 Myr, such as mechanical feedback or gravitational shocks, that are not included in the Cloudy calculations. The implications of the inferred cooling times are discussed further in Section 5.7.

5.2 Alternative scenarios for the extended soft X-ray emission

We developed analogous models to those of the previous section in which the hot thermal gas is replaced by cold/warm clouds photoionized by the QSO, exploring a grid of gas densities and metallicities. As expected, the emission from these photoionized clouds fails to account for the observed extended X-ray excess, even under the extreme assumption that the warm/cold phase occupies all the volume with filling factor $f_{V}=1$. A more detailed analysis of these models, including constraints derived from the Ly$\alpha$  nebula, will be presented in a future paper, as this lies beyond the scope of the present work.

We tested a non-thermal IC scenario, fitting the X-ray excess with a power law, which results in unphysically steep photon indices of $\Gamma\sim 6$ (see Worrall et al., 2016), and found no evidence of powerful radio jets associated with this QSO (e.g., Sydney University Molonglo Sky Survey at 843 MHz; Mauch et al. 2003, and 0.8 mJy/beam at 1.367 GHz and 0.78 mJy/beam at 887.5 MHz in the Rapid ASKAP Continuum Survey).

We also explored a scenario in which the extended X-ray emission originates from thermal Compton up-scattering of AGN seed photons by a hot electron population produced by an extended AGN wind, modeled using compTT. To reproduce the spectral shape of the extended X-ray emission, we assumed, as a conservative scenario, seed photons representative of the standard soft X-ray excess commonly observed in local AGN, typically peaking around $\sim 0.1\leavevmode\nobreak\ \rm keV$ rest-frame (Crummy et al., 2006; Done et al., 2012). We adopted an optical depth of $\tau=0.01$, which is the upper limit derivable from our constraint on the electron column density in the X-ray quasar spectrum, corresponding to an electron density of $n_{e}=\tau/(l\,\sigma_{\rm T})\sim 0.3\leavevmode\nobreak\ \rm cm^{-3}$, where $l\sim 15\leavevmode\nobreak\ \rm kpc$ and $\sigma_{\rm T}=6.65\times 10^{-25}\leavevmode\nobreak\ \rm cm^{2}$ is the Thomson cross-section. From spectral fitting, we obtained $\log(\texttt{norm}_{\texttt{compTT}})=-2.77^{+0.71}_{-0.88}$ and $kT_{e}=47^{+17}_{-11}\leavevmode\nobreak\ \rm keV$. Although this model can, in principle, reproduce the observed extended emission with a required $L_{\text{seed}}\approx 2\times 10^{46}\leavevmode\nobreak\ \rm erg\leavevmode\nobreak\ s^{-1}$, about a factor of two lower than the UV luminosity at 400$\,\AA$, we regard it as energetically implausible, as it requires the presence of a hot ($kT_{e}\sim 47\leavevmode\nobreak\ \rm keV$) and dense AGN wind that is capable of reaching distances of several tens of kiloparsecs from the nucleus. Although such AGN winds have been invoked by observations and numerical simulations (Faucher-Giguère & Quataert, 2012; Zubovas & King, 2012; Weinberger et al., 2017; Powell et al., 2018), these are typically confined within the galactic nucleus or the galaxy ISM, given their short cooling times. In addition, we do not see any signs of such fast and extended wind in the warm gas kinematics, as the Ly$\alpha$  nebula appears kinematically quiet (Borisova et al., 2016) without detectable CIV emission. Similarly, we do not see any Broad Absorption Line in the quasar spectra. For all these reasons, we consider a Compton up-scattering origin of the detected X-ray emission as unlikely.

5.3 Constraints from SZ non-detection

The detection of extended X-ray emission suggests the presence of hot gas within a $\sim 3\times 10^{13}\leavevmode\nobreak\ \rm M_{\odot}$ halo, potentially producing a thermal Sunyaev-Zeldovich (SZ) effect (Sunyaev & Zeldovich, 1972). However, an analysis of ALMA Band 3 data (Pensabene et al., 2024) reveals no clear evidence of an SZ signature in the current observations. To estimate the expected SZ significance ($\sigma_{\text{SZ}}$), we generate the SZ signature for each set of $\beta$-model parameters, $kT$, and $n_{\text{e}}$ from the previous analyses, and inject the resulting models into jackknifed realizations of the available ALMA data (see Di Mascolo et al., 2023; van Marrewijk et al., 2025, for details on the jackknifing procedure). Our analysis indicates that the non-detection within the sensitivity limits of the current ALMA data is fully consistent with the hot gas properties as derived by the X-ray analysis. The posterior parameters from the MCMC analysis predict a marginal SZ signal with a significance of only $\simeq(2.6\pm 0.4)\,\sigma_{\text{SZ}}$. Additionally, continuum emission from companion quasars and deviations from thermal pressure support could further suppress the SZ signal, potentially explaining the current non-detection. New ALMA Cycle 12 observations (Project ID: 2025.1.00107.S, priority grade B) have been awarded to our team. These data, once executed, will provide the sensitivity and angular resolution needed to further probe the thermal and non-thermal components of the ICM and to test more robustly the presence of an SZ signal in this system.

5.4 Comparison with Spiderweb and simulated hot halos

Studying the nascent hot phase of the CGM (or proto-ICM) in $z>2$ halos is difficult because it requires long integration times with current X-ray telescopes. Previous to this work, the only reported direct imaging detection of extended thermal ICM emission at $z>2$ was around the Spiderweb galaxy in $\approx 700\leavevmode\nobreak\ \rm ks$ L24, who reported residual extended emission after the subtraction of the X-ray Inverse Compton emission due to radio jets. For both Spiderweb and MQN01, the extended thermal emission appears to originate from the halo of a proto-BCG at the center of a forming (or proto-)cluster.

Comparing luminosity, density, and morphology of these systems at $z=2.16$ and $z=3.25$ could provide some insights into the possible evolution of the hot CGM/proto-ICM with redshift and with their halo properties. Following the approach of L24, we computed the average SB of the extended X-ray emission in concentric annuli using:

[TABLE]

where $A_{\text{ann}}$ is the area of the annulus, $t_{\text{exp}}$ is the total exposure time, and the ratio $ExpMap_{\text{max}}/ExpMap_{\text{ann}}\approx 1$ accounts for small variations of the Chandra sensitivity within the relevant portion of the field of view. The energy conversion factor (ECF) depends on the temperature and metallicity parameters of the xsmekal model used to fit the spectra. The extended X-ray SB profile was obtained by subtracting the PSF-derived profile from the data. We estimated the SB profile in the 0.5-2 keV band (hereafter $\text{SB}_{X}(r)$) for the extended X-ray emission around QSO ID1 in the MQN01 protocluster. This used an ECF of $\sim 4.6\times 10^{-11}\leavevmode\nobreak\ \rm erg\leavevmode\nobreak\ cm^{-2}\leavevmode\nobreak\ cnts^{-1}$, based on the best-fit temperature ($kT\sim 1.8\leavevmode\nobreak\ \rm keV$), and $0.014\,Z_{\odot}$ from our MCMC fit (see Section 4). As indicated by the MCMC results, the PSF model was normalized to match the observed counts within the central 2$\arcsec$ (i.e., $\sim$15 kpc), under the assumption that $\sim 12\,\%$ of the total emission originates from thermal processes (i.e., $f_{\rm th}=0.12$; see Table 2).

Figure 9 shows the redshift-dimming corrected $\text{SB}_{X}(r)$ of the extended X-ray emission in MQN01 (red) and Spiderweb (blue). The red dots, indicating the $\text{SB}_{X}(r)$ of the MQN01 halo, consist of three data points measured between 15 and 30 kpc (i.e., $2\arcsec$-$4\arcsec$; where the extended residual X-ray emission is detected following the subtraction of the nuclear emission), one point within 15 kpc ($<2\arcsec$; shaded grey region), which is extrapolated from the best-fitting $\beta$ model (see Section 4), and upper limits at all radii above 30 kpc. The $\text{SB}_{X}$ profile of the X-ray extended halo in Spiderweb is taken from L24. To correct for different rest-frame bands due to redshift, we applied a 0.536 factor to the Spiderweb $\text{SB}_{X}(r)$ halo. This factor was derived assuming the best-fit $kT=2\leavevmode\nobreak\ \rm keV$ xsmekal model reported in Tozzi et al. (2022). The black lines show the $\text{SB}_{X}(r)$ profiles of simulated hot halos, where the solid, dashed, and dotted lines represent the median, 16th-84th percentile, and full range (minimum to maximum) of the $\text{SB}_{X}$ profiles, respectively. These profiles were computed from a sample of 12 halos at $z=3$ extracted from the DIANOGA cosmological hydrodynamical simulations of galaxy clusters (Esposito et al., 2025). This version of the DIANOGA simulations includes 14 target clusters at $z=0$, in the mass range $(0.2-3)\times 10^{15}\leavevmode\nobreak\ \rm M_{\odot}$, simulated with OpenGADGET-3 (e.g., Groth et al., 2023; Damiano et al., 2024). The 12 halos were selected in the mass range $(2-6)\times 10^{13}\leavevmode\nobreak\ \rm M_{\odot}$ to have a median mass matching the virial mass of the MQN01 halo ($M_{\rm vir}\sim 3\times 10^{13}\leavevmode\nobreak\ M_{\odot}$). The X-ray maps of the simulated halos were obtained with the post-processing tool SMAC (Dolag et al., 2005) in cubic regions of 500 kpc per side around the center of each halo, with 1 kpc resolution.

While the $\text{SB}_{X}$ profiles of MQN01 and Spiderweb differ significantly in normalization, by a factor of 3-10 within 30 kpc from the QSOs, MQN01’s profile is also noticeably steeper, consistent with a highly compact core. The $\text{SB}_{X}$ profile in MQN01 is consistent, for $r>20\leavevmode\nobreak\ \rm kpc$, with the upper envelope of the $\text{SB}_{X}$ profiles from DIANOGA simulations. Similarly, the $\text{SB}_{X}$ profile of the Spiderweb halo closely follows the 84th percentile of the distribution from the DIANOGA hot halos. Within the limited statistics, simulations consistently predict fainter and shallower $\text{SB}_{X}$ profiles than observed. Part of the reason for this discrepancy may lie in the different sub-grid prescriptions adopted in the simulations. In particular, the DIANOGA runs use a relatively low density threshold for star formation ($0.13\leavevmode\nobreak\ \rm cm^{-3}$), which inhibits the formation of high-density gas regions, especially near halo centers. To test the impact of this, we compared the hot gas density profiles from DIANOGA with those of the most massive halos at $z=3$ in the DaLya simulations (priv. comm., Lazeyras et al. in prep.), which adopt a star formation threshold 100 times higher. For a consistent comparison, we applied the same temperature cut ($kT>0.1\leavevmode\nobreak\ \rm keV$) as in the DIANOGA analysis to isolate the X-ray-emitting hot phase in both simulations. We found that the density of $kT\sim 0.1\leavevmode\nobreak\ \rm keV$ gas in DaLya is 5-10 times higher than in DIANOGA. While other factors, such as feedback modeling, numerical resolution, or selection effects in the observations, may also contribute to the discrepancy, the low star formation threshold in DIANOGA appears to suppress central gas densities and flatten the hot gas profiles (see Figure 10), thus significantly reducing the X-ray emissivity of the hot CGM.

Figure 10 compares electron density profiles, $n_{\text{e}}(r)$, of hot halos in MQN01 (red line and black dots), Spiderweb (blue), and DIANOGA simulations (gray). The red curve shows the best-fit $n_{\text{e}}(r)$ for the MQN01 halo, with the shaded area indicating the 68% confidence interval from Monte Carlo realizations on $n_{\text{e},0}$, $\beta$, and $r_{\text{core}}$. The $n_{\text{e}}(r)$ of the Spiderweb hot halo is from L24, while DIANOGA $n_{\text{e}}(r)$ are directly extracted from simulations.

Between 15 and 30 kpc, the electron density of the MQN01 hot halo exceeds that of the Spiderweb halo by a factor of $4.3\pm 1.0$ on average. At 15 kpc the ratio is $4.9\pm 3.1$, decreasing to $3.5\pm 4.5$ at 30 kpc (uncertainties propagated from the MQN01 $n_{\rm e}$ and the Spiderweb bounds).

These SBX and $n_{\rm e}(r)$ differences broadly match expectations from cosmological evolution: in a self-similar scenario at fixed halo mass $n\propto(1+z)^{3}$ (Kaiser, 1986), and since SB${}_{X}\propto\int n_{\rm e}n_{\rm H}\Lambda(T,Z)\,dl$ one expects SB${}_{X}\propto(1+z)^{6}$ for fixed $T$ and $Z$. For MQN01 ($z=3.25$) vs. Spiderweb ($z\sim 2.2$) this yields a factor $\sim 5.5$ in $\text{SB}_{X}$ and $\sim 2.34$ in $n_{\rm e}$. After accounting for this self-similar redshift scaling, the mean density ratio between 15-30 kpc range is reduced to $1.84\pm 0.42$, with corrected point estimates of $2.08\pm 1.32$ at 15 kpc and $1.48\pm 1.91$ at 30 kpc (errors adjusted accordingly). The ratio at 30 kpc remains poorly constrained due to large uncertainties. To account for unresolved gas inhomogeneities, we model the electron density at each radius as $n_{\text{e}}(r)=C^{1/2}\langle n_{\text{e}}\rangle$, where $C$ is the clumping factor (e.g., Ettori et al., 2013). Reconciling the residual differences at 15 kpc would thus require modest clumping in MQN01 relative to the Spiderweb, with $C\simeq 4_{-3}^{+8}$. This range is consistent with typical clumping factors found in local ICM ($\leq 5$; Nagai & Lau, 2011; Vazza et al., 2013; Eckert et al., 2015; Mirakhor & Walker, 2021). However, the notably high central densities in MQN01 might also reflect additional effects, such as localized gas condensation or AGN-driven compression, which could enhance central gas densities beyond what is expected from redshift evolution alone.

5.5 Pressure confinement of cold CGM clouds

The hot gas pressure at 15 kpc and 30 kpc, that are the minimum and maximum radii at which we directly observe the extended X-ray emission based on our analysis, is $0.92_{-0.63}^{+1.24}\leavevmode\nobreak\ \rm keV\leavevmode\nobreak\ cm^{-3}$ and $0.30_{-0.12}^{+0.53}\leavevmode\nobreak\ \rm keV\leavevmode\nobreak\ cm^{-3}$, respectively. These are about one to two orders of magnitude higher than typical values observed for the ICM in local galaxy clusters (Ettori et al., 2009), and about three/nine times higher than those obtained for the hot halo in Spiderweb i.e., $<0.1\leavevmode\nobreak\ \rm keV\leavevmode\nobreak\ cm^{-3}$ L24 from independent measurements using ALMA SZ observation (Di Mascolo et al., 2023). Such a high thermal pressure is sufficient to confine the cold, Ly$\alpha$ -emitting phase of the CGM. Assuming pressure equilibrium and a cold gas temperature of $T_{\rm{cold}}=5\times 10^{4}\rm\leavevmode\nobreak\ K$, the implied temperature ratio $T_{\rm{hot}}/T_{\rm{cold}}\sim 400$ leads to a density contrast of the same factor. Therefore, for a hot gas density of $n_{\rm{hot}}=0.1\leavevmode\nobreak\ \rm cm^{-3}$ at $\sim 50\leavevmode\nobreak\ \rm kpc$ (see Figure 10), the warm gas can reach densities up to $n_{\rm{cold}}\sim 40\leavevmode\nobreak\ \rm cm^{-3}$, supporting a scenario in which the Ly$\alpha$  nebula in MQN01 arises from a distribution of dense clumps pressure-confined by the surrounding hot halo (Pezzulli & Cantalupo, 2019; Cantalupo et al., 2019), at least in the CGM inner regions. We stress, however, that our results (in particular the exceptionally high densities) may not be representative of the average Ly$\alpha$  nebulae detected so far around luminous $z\sim 3$ quasars, as the latter do not generally reside in overdensities as large as that of MQN01 and their halo masses are expected to be, on average, smaller than what is derived for MQN01 (Pezzulli & Cantalupo, 2019; de Beer et al., 2023).

5.6 Comparison with local galaxy groups and clusters

Figure 11 shows the relation between $L_{0.5-2\leavevmode\nobreak\ \rm keV}$ and $kT$ for a variety of massive structures, including low- and high-redshift systems ranging from galaxy groups to rich clusters. Notably, the red square represents the average temperature and X-ray luminosity of the extended X-ray emission in the Spiderweb protocluster ($kT=2.0_{-0.4}^{+0.7}\leavevmode\nobreak\ \rm keV$, $L_{0.5-2\leavevmode\nobreak\ \rm keV}=(2.0\pm 0.5)\times 10^{44}\leavevmode\nobreak\ \rm erg\leavevmode\nobreak\ s^{-1}$). The Spiderweb point falls within the scatter of the group/cluster population, although the temperature is likely overestimated compared to the SZ-based value ($kT_{SZ}=(0.7\pm 0.3)\leavevmode\nobreak\ \rm keV$ Di Mascolo et al., 2023). Adopting the lower SZ-based temperature would shift the Spiderweb point even further above the typical relation. The extended X-ray emission analyzed here at $z=3.25$ (red circle) shows a notably higher 0.5-2 keV luminosity than other structures with comparable temperatures, even after accounting for uncertainties on the $f_{th}$ value, placing MQN01 well outside the typical distribution of hot gas in groups and clusters.

One possible explanation for the observed high luminosity could be related to redshift evolution. In particular, the self-similar evolution model (see Maughan et al., 2012, where $L_{X}\propto E(z)$ as first approximation) predicts luminosities of $\sim 6.3\times 10^{43}\leavevmode\nobreak\ \rm erg\leavevmode\nobreak\ s^{-1}$ for Spiderweb and $\sim 4.1\times 10^{44}\leavevmode\nobreak\ \rm erg\leavevmode\nobreak\ s^{-1}$ for MQN01. This provides a baseline to compare with local groups and clusters, under the assumption that the self-similar model applies. However, we note that this evolution is formally derived for bolometric luminosities, whereas our comparison is carried out in the soft band. Since the fraction of the bolometric luminosity falling in the soft band depends on the plasma temperature, applying the self-similar scaling directly to soft-band luminosities can overestimate the true redshift evolution, particularly for low-temperature systems. Even accounting for this, the MQN01 point remains significantly above the local $L_{X}$-$T$ population, while the Spiderweb point appears broadly consistent with it.

Within the framework of thermal bremsstrahlung and line emission, this is not surprising, since we are probably observing a distinct and earlier phase of the ICM, which is characterized by different physical conditions. This pronounced deviation suggests that the enhanced luminosity of MQN01’s diffuse X-ray gas may be due to higher gas densities at fixed mass (i.e., greater compactness), boosting emissivity. This compactness likely results from its higher redshift combined with a relatively steep gas density profile, supported by the $\beta$ value from our MCMC analysis (Table 1), which exceeds typical local ICM values of 0.1-1 (Mohr et al., 1999; Dong et al., 2010; Conte et al., 2011; Paggi et al., 2021; Xue & Wu, 2000; Wise et al., 2004; Mirakhor et al., 2022).

5.7 Thermal instability and hot halo condensation

In this work, we have presented the first detection of extended hot ($T\gtrsim 10^{7}K$) gas in a massive halo at $z>3$. Cosmological theory (e.g., Dekel & Birnboim, 2006) predicts that such structures should form as a consequence of the virialization of gas infalling onto dark matter haloes from the IGM and being shock-heated as a result of the collision of multiple infalling gas streams. The early cosmic time of our observation, along with with the peculiarly compact structure (see Section 5.4), leads us to speculate that we may be witnessing the initial stages of the formation and virialization of a hot halo, possibly leading, through further gas accretion and the expansion of the shock front to larger radii, to the more extended hot gas structures typically observed in the nearby Universe.

At the same time, the extreme luminosity of MQN01 (see Figure 11), indicative of large radiative losses, also raises the question of whether some of the recently virialized gas may already be undergoing substantial cooling, possibly contributing to the fueling of star formation and AGN activity in the central galaxy. A similar interpretation has been discussed by L24 in the context of the hot halo of Spiderweb at $z\simeq 2.16$, based on the inferred low temperatures ($T\sim 0.3$-$0.9\leavevmode\nobreak\ \rm keV$) and short cooling times ($<100\leavevmode\nobreak\ \rm Myr$)444Shorter than those typically observed in cluster cores at low redshift ($\approx 10^{3}\leavevmode\nobreak\ \rm Myr$; Edge et al., 1992; Hudson et al., 2010). In addition, L24 report a nominal isobaric mass deposition rate of $\dot{M}_{\text{cool}}\approx 0.25$-$1\times 10^{3}\leavevmode\nobreak\ \rm M_{\odot}\leavevmode\nobreak\ yr^{-1}$, which exceeds typical cooling rates inferred for low-redshift ICMs (McDonald et al., 2018), although it does not account for heating or feedback processes.

Similarly, applying the same steady-state, isobaric cooling-flow formalism to MQN01 (Equation (2) Fabian, 1994), and using $L_{\text{cool}}\approx L_{\rm 0.5-2\,keV}^{15-23\rm\leavevmode\nobreak\ kpc}=5.8_{-1.4}^{+2.1}\times 10^{44}\leavevmode\nobreak\ \rm erg\leavevmode\nobreak\ s^{-1}$ (as in Jones & Forman, 1984), would result in a mass deposition rate of $\dot{M}_{\text{cool}}\simeq 1.3_{-1.1}^{+1.7}\times 10^{3}\leavevmode\nobreak\ \rm M_{\odot}\leavevmode\nobreak\ yr^{-1}$ in this central region. This estimate assumes that the entire X-ray-emitting gas cools efficiently and steadily, and should be interpreted as an upper limit. We emphasise that these estimates ignore any source of heating or mechanical energy injection (e.g., AGN feedback, turbulence), which are thought to largely compensate for the cooling, so that the nominal mass deposition rates can be overestimated by factors 10-100 (e.g., Fabian, 1994; Gaspari et al., 2013). Moreover, a sustained mass accretion rate of $\simeq 10^{3}\leavevmode\nobreak\ \rm M_{\odot}\leavevmode\nobreak\ yr^{-1}$ could at most occur episodically. Indeed, if maintained for 1 Gyr, it would lead to an accumulation of $10^{12}\leavevmode\nobreak\ \rm M_{\odot}$ of cold gas, exceeding that observed in any known galaxy.

To assess whether a less extreme scenario of localized thermal condensation is likely in our hot halo in MQN01, we compared the radiative cooling time, $t_{\rm cool}$, to the dynamical free-fall time, $t_{\rm ff}$, at $r=$15 and 30 kpc. The ratio $t_{\rm cool}/t_{\rm ff}$ is a widely used diagnostic for precipitation/condensation in hot halos (e.g. Voit et al., 2015; Choudhury & Sharma, 2016; Voit et al., 2017; Stern et al., 2021; Donahue & Voit, 2022). Using the ’classical’ formula (as shown by Donahue & Voit 2022) and a representative cooling function $\Lambda(T,Z)\simeq 10^{-23}\rm erg\leavevmode\nobreak\ cm^{3}\leavevmode\nobreak\ s^{-1}$, we found $t_{\rm cool}(15\,{\rm kpc})=55_{-14}^{+24}$ Myr and $t_{\rm cool}(30\,{\rm kpc})=181_{-83}^{+268}$ Myr. Estimating $t_{\rm ff}$ from the local gravitational field (e.g., McCourt et al., 2012; Wibking et al., 2025), which we computed assuming a NFW potential with a concencentration $c\simeq 3.5\pm 0.5$ (as suggested by Correa et al., 2015, for a $10^{13}\leavevmode\nobreak\ \rm M_{\odot}$ halo at $z\approx 3$), we obtained $t_{\rm ff}(15\,{\rm kpc})=36_{-4}^{+5}$ Myr and $t_{\rm ff}(30\,{\rm kpc})=58_{-5}^{+8}$ Myr. These yield ratios $t_{\rm cool}/t_{\rm ff}(15\,{\rm kpc})=1.9_{-0.9}^{+1.9}$ and $t_{\rm cool}/t_{\rm ff}(30\,{\rm kpc})=3.1_{-1.4}^{+4.6}$. All fall within the commonly cited precipitation threshold, $1<t_{\rm cool}/t_{\rm ff}<10$, under which local thermal instabilities can seed multiphase condensation in roughly hydrostatic, feedback-regulated atmospheres (Donahue & Voit, 2022).

We note that, despite its common usage in the literature, the free-fall time $t_{\rm ff}$ is not strictly the most relevant timescale for local thermal instabilities. A more direct indicator of potential condensation is the Brunt-Väisälä time $t_{\text{BV}}$, which is the timescale for the reaction of buoyancy against thermal instability. While $t_{\text{ff}}$ and $t_{\text{BV}}$ are often comparable in compact halo cores, $t_{\text{BV}}$ remains a more precise diagnostic because it is defined locally through perturbation analysis (e.g., Nipoti & Posti, 2014; Wibking et al., 2025). To compute $t_{\text{BV}}$, we used Equation (D3) in Appendix D, again assuming the gravitational field $g$ associated with the NFW potential. We found $t_{\rm cool}/t_{\text{BV}}(15\,{\rm kpc})=1.3_{-0.5}^{+1.0}$,and $t_{\rm cool}/t_{\text{BV}}(30\,{\rm kpc})=4.5_{-2.4}^{+8.7}$, which are still in the range $1<t_{\rm cool}/t_{\rm BV}<10$.

This does not imply a global cooling catastrophe, but it is consistent with intermittent condensation in portions of the hot halo that can feed cold inflows and/or SMBH fueling. Note, however, that classic analyses demonstrate that buoyancy in atmospheres with positive entropy gradients can suppress linear thermal instability (Binney et al., 2009). Recent simulations (see Wibking et al., 2025) further suggest that whether condensation occurs depends on turbulence, mixing, and the geometry of heating. Some systems with $t_{\rm cool}/t_{\rm ff}\lesssim 10$ can remain stable if turbulent support is weak or heating is centrally concentrated (see discussions in Donahue & Voit, 2022). Consequently, our results should be interpreted as strong plausibility for localized condensation, contingent on the local balance of turbulence, mixing, and feedback.

To further probe the dynamical state of the hot gas halo, we examined the balance between the pressure gradient force and the gravitational acceleration derived from the NFW dark matter profile. Using the MCMC posterior parameters of $kT$, $r_{\text{core}}$, and $M_{\text{vir}}$ with their associated uncertainties, we computed the ratio of the hydrostatic acceleration from the pressure gradient, $1/\rho(dP/dr)$, to the gravitational acceleration, $g(r)$, at radii of 15 and 30 kpc. Our Monte Carlo analysis, sampling the uncertainties, yields ratios of $\simeq 0.7_{-0.5}^{+1.4}$ at 15 kpc and $\simeq 1.3_{-0.8}^{+1.7}$ at 30 kpc, consistent with hydrostatic equilibrium, although significant deviations are allowed within the credible intervals. These results suggest that, while the hot gas is roughly in hydrostatic balance, there may be local departures indicative of dynamical processes such as inflows, outflows, or turbulence. Such deviations can enhance thermal instability and promote condensation, consistent with the $1<t_{\text{cool}}/t_{\text{ff}}<10$ found above. Therefore, the combined diagnostics support a scenario in which the hot halo is marginally stable, possibly allowing localized cold gas condensation and potentially fueling galaxy growth and AGN activity.

6 Summary and Conclusions

In this paper, we report a significant (at least $8\sigma$ between 15-30 kpc) detection of 0.5-2 keV ($\sim$2-8 keV rest-frame) X-ray emission, extended out to $\sim 30\leavevmode\nobreak\ \rm kpc$, around the brightest QSO (ID1) within the MQN01 protocluster at $z=3.25$ (Section 3). This QSO is also surrounded by a giant ($>200\leavevmode\nobreak\ \rm kpc$) Ly$\alpha$  nebula (#1 in Borisova et al., 2016). The extended X-ray emission exhibits an isotropic morphology and appears spatially correlated with the Ly$\alpha$  nebula (Section 3.2).

The steepness of the spectrum and the isotropic morphology of the extended X-ray emission suggest a thermal origin, with the emission powered by bremsstrahlung and collisional excitation mechanisms. Alternative emission mechanisms, such as photoionization of cold or hot gas clouds by the AGN, inverse Compton scattering, or Compton up-scattering, fail to fully or even partially reproduce the observed extended X-ray excess, as discussed in Section 5.2. Therefore, this likely represents the first evidence of thermal emission from proto-ICM (or hot CGM) at $z>3$.

Assuming the QSO halo is virialized and its X-ray emission is characterized by a $\beta$-model, we performed a joint spatial and spectral MCMC analysis. The median posterior values imply the following (Section 4):

• •

A temperature of $\approx 1.8\leavevmode\nobreak\ \rm keV$ corresponding to a virial halo mass of $M_{\mathrm{vir}}\simeq(3\pm 1)\times 10^{13}\leavevmode\nobreak\ M_{\odot}$.

• •

Hot gas within 30 kpc exhibits electron densities ranging from $0.9$ to $0.2\leavevmode\nobreak\ \rm cm^{-3}$ and a 0.5-2 keV X-ray luminosity of $L_{\text{0.5-2 keV}}=2.25_{-1.38}^{+0.77}\times 10^{45}\leavevmode\nobreak\ \rm erg\leavevmode\nobreak\ s^{-1}$.

• •

The spatial fit yields $\beta$-model parameters of core radius $r_{\text{core}}\simeq 36_{-13}^{+16}\leavevmode\nobreak\ \rm kpc$ and $\beta\simeq 2.0_{-0.8}^{+1.4}$. The latter is significantly higher than the typical $\beta=0.1$-$1$ observed in local galaxy groups and clusters, indicating a steeper gas density profile.

• •

The analysis implies that the hot gas in the halo of MQN01 contains a substantial fraction $f_{\textrm{hot}}=56_{-20}^{+65}\%$ of the baryons that are theoretically associated to the halo (i.e. the baryonic mass that one would expect assuming a cosmic baryon fraction $\Omega_{b}/\Omega_{m}\simeq 0.15)$.

Despite the large uncertainty, the last result, united to the fact that some baryonic mass must also be stored in the colder (Ly$\alpha$-emitting) phase of the CGM, suggests that the halo of MQN01 contains a substantial fraction of its theoretical baryon budget within the virial radius. We plan to better quantify the total baryon budget of MQN01 in a future combined analysis of both the X-ray and Ly$\alpha$ emission.

Based on our results, the currently available ALMA data are not expected to allow a clear detection of any thermal Sunyaev-Zeldovich (SZ) signal, as the expected significance is only $(2.6\pm 0.4)\,\sigma_{\text{SZ}}$, even assuming perfect subtraction of contaminating sources (Section 5.3). However, upcoming tailored ALMA observations will be able to confirm (or challenge) our findings and add further constraints on the spatial distribution and thermal state of the hot gas in MQN01.

MQN01’s hot halo is notably brighter and denser than that in the Spiderweb protocluster at $z=2.16$ (L24), with X-ray surface brightness and electron density profiles exceeding those of Spiderweb by factors of 3-10 and 3-5 times within 30 kpc. These differences likely reflect cosmic evolution in density with redshift and potentially increased gas clumpiness. Compared to cosmological simulations (DIANOGA), MQN01’s halo represents the luminous extreme of predicted hot halos, with higher densities possibly linked to a too low density threshold of star formation assumed in the simulations, which prevents dense gas from remaining in the X-ray emitting phase (Section 5.4).

On the other hand, the total X-ray luminosity within 30 kpc ($L_{0.5-2\leavevmode\nobreak\ \rm keV}\approx 2.3\times 10^{45}\rm\leavevmode\nobreak\ erg\leavevmode\nobreak\ s^{-1}$) of the MQN01 halo is unusually high for its temperature, placing it well above the $L_{X}$-$kT$ relation for local groups and clusters. Even after correcting for self-similar evolution to $z=0$, MQN01 remains an outlier, reflecting a uniquely compact, dense gas distribution during an early intracluster medium formation phase. This enhanced emissivity likely arises from a steep density profile and concentrated gas within its dark matter halo, setting MQN01 apart from typical systems across cosmic time (Section 5.6).

At roughly 15 kpc, MQN01’s hot gas exhibits low temperature ($<2\leavevmode\nobreak\ \rm keV$), short cooling times ($\approx 55\leavevmode\nobreak\ \rm Myr$), and relatively low entropies ($K\simeq 3\leavevmode\nobreak\ \rm keV\leavevmode\nobreak\ cm^{-2}$), and a relatively low cooling time to free-fall (or to Brunt-Väisälä time) time ratio $t_{\rm cool}/t_{\rm ff}\simeq 1.9_{-0.9}^{+1.9}$ ($t_{\rm cool}/t_{\rm BV}\simeq 1.3_{-0.5}^{+1.0}$), consistent with localized cold gas condensation according to some instability criteria proposed in the literature (Section 5.7). On the other hand, a classical, large-scale cooling flow is deemed unlikely, but we report that the estimated mass deposition rate, assuming steady-state isobaric cooling, would be $\dot{M}_{\rm cool}\simeq 1.3_{-1.1}^{+1.7}\times 10^{3}\leavevmode\nobreak\ \rm M_{\odot}\leavevmode\nobreak\ yr^{-1}$, an extreme upper limit given likely heating and non-radiative effects (Section 5.7).

The hot gas pressure at 15 kpc and 30 kpc ($\sim 0.30_{-0.12}^{+0.53}$ and $\sim 0.92_{-0.63}^{+1.24}\leavevmode\nobreak\ \rm keV\leavevmode\nobreak\ cm^{-3}$) are exceptionally high, one to two orders of magnitude above local cluster values and three to nine times higher than in the Spiderweb protocluster. Such pressure can efficiently pressure-confine denser Ly$\alpha$ -emitting clouds, supporting a scenario for which the extended Ly$\alpha$  nebula MQN01 could originate from colder clumps embedded within a hot, high-pressure medium.

Multi-wavelength follow-up observations with next-generation X-ray telescopes, such as the Advanced X-ray Imaging Satellite (AXIS) and NewAthena (Advanced Telescope for High ENergy Astrophysics), will be crucial to better characterize the physical properties of the extended gas distribution within the MQN01 Cosmic Web node. These will be enabled by their expected high spatial resolutions ($<1.5\arcsec$ and $<3\arcsec$) and spectral resolutions ($<70\leavevmode\nobreak\ \rm eV$ and $<1.5\leavevmode\nobreak\ \rm eV$ at 1 keV) in the soft X-ray band (0.3-10 keV and 0.2-12 keV) as reported in Reynolds et al. (2023) and Barret et al. (2023), respectively. ALMA observations to detect the expected SZ signal on small scales are scheduled in Cycle 12 (Proposal ID 2025.1.01488.S). Very Long Baseline Interferometry (VLBI) radio observations may also prove valuable for exploring potential jet activity. These investigations will be key to either confirming or challenging the thermal emission scenario proposed in this work, and to further constraining the physical properties of the nascent hot CGM/proto-ICM in this exceptional system.

Acknowledgements.

This project was supported by the European Research Council (ERC) Consolidator Grant 864361 (CosmicWeb). FF acknowledges support by HORIZON2020: AHEAD2020-Grant Agreement n. 871158. PT acknowledges support from the Next Generation European Union PRIN 2022 20225E4SY5 - ”From ProtoClusters to Clusters in one Gyr.”. FV acknowledges support from the ”INAF Ricerca Fondamentale 2023 Large GO” grant. LDM was supported by the French government, through the UCAJ.E.D.I. Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01. RM acknowledges financial support from the ASI-INAF agreement n. 2022-14-HH.0. The authors also thank Alessandro Lupi, Silvano Molendi, Piero Rosati, and Ákos Bogdán for useful discussions.

Appendix A Alignment of QSO ID1

Since the primary aim of this paper is to investigate potential extended X-ray emission around QSOs, with a particular focus on QSO ID1 due to its high photon count, we developed a refined alignment and merging procedure for the individual ObsIDs. This method minimizes misalignment in the PSF of the ID1 source across different observations, allowing for a more reliable search for extended X-ray emission in the merged event file. Section 2.2 contains a more detailed description of this method. Figure 12 presents the results of this alignment, showing the 0.5-2 keV images of the ObsIDs, centered on the QSO ID1, after re-projection and Gaussian smoothing. The green circles mark 2$\arcsec$-radius regions centered on a common reference coordinate corresponding to the QSO position. This alignment, along with additional tests discussed in the main text, confirms that the observed extended emission is not an artifact caused by ObsID misalignment or inaccuracies in PSF modeling.

Appendix B Lack of Extended X-ray Emission around the others AGN within MQN01

Figure 13 presents the radial profiles extracted from observational data (black points) compared to those derived from the simulated PSF (red points) for other X-ray AGNs reported in Travascio et al. (2024). AGNs ID3 and ID4, located near the brightest QSO (ID1), were excluded due to their proximity, which complicates the identification of extended emission. As reported in Section 3.1, no significant extended emission is detected around these AGNs, as the normalized PSF accurately reproduces the observed profiles within the uncertainties. Beyond 4$\arcsec$ (corresponding to 30 kpc at $z=3.25$), the uncertainties are primarily dominated by the local background. This result further validates the reliability of our analysis and the accuracy of the generated PSFs in reproducing the observational data.

To further test this, we estimated the ratio of net counts in the 0.5-1.5 keV energy range (where the extended emission is detected) within the 2$\arcsec$-3$\arcsec$ annulus and the 2$\arcsec$-radius region centered on ID1. This was then compared to the ratios for other bright X-ray point sources in the Chandra field of view (FoV), located as close as possible to ID1. Among these, the only source with a well-defined, non-elongated PSF, a high count rate, and a position near the aim point is the AGN ID2. To estimate the net counts, the background was derived from 100 large annuli around the sources, with their minimum and maximum radii chosen randomly within specific limits to minimize the influence of local background variations. By computing the median, along with the 25th and 75th percentiles, of the count ratio for ID1 and ID2, we found that the excess emission in ID1 compared to ID2 is detected at a significance level of $4.5\sigma$ (ranging between $4\sigma$ and $5\sigma$). This calculation assumes a similar PSF for both sources, which is reasonable given their proximity.

Figure 14 displays the 0.5-2 keV images for both the observational data (left) and the simulated PSF (right). This comparison highlights the excess counts beyond 2$\arcsec$ around QSO ID1, as also shown in the left panel of Figure 2.

Appendix C Proof of equation (7)

Equation (7) can be derived for instance from 3.241-4 in Gradshteyn & Ryzhik (2007). However, the original source, containing the mathematical proof, is of uneasy retrieval. Furthermore, it is unclear from Gradshteyn & Ryzhik (2007) whether the formula has been proven for non-integer values of $a$. We therefore provide here a simple, explicit proof of equation (7), valid for all $a>0.5$:

[TABLE]

where we used the substitution $y=1/(1+x^{2})$, then the definition of Euler’s beta function $B$ and finally the well known properties $B(z_{1},z_{2})=\Gamma(z_{1})\Gamma(z_{2})/\Gamma(z_{1}+z_{2})$ and $\Gamma(\frac{1}{2})=\sqrt{\pi}$. Note that for $a\leq 0.5$ the integral diverges and $\tilde{\chi}(a)$ is undefined, or formally infinite.

Appendix D Brunt-Väisälä time-scale for an isothermal $\beta$ model

Following Wibking et al. (2025), we can define the Brunt-Väisälä time-scale as:

[TABLE]

where $g$ in the gravitational field and $h_{S}$ is the local entropy scale-length defined as $h_{S}=(dS/dr)^{-1}$, with $S=\ln(Tn^{1-\gamma})$. For an isothermal $\beta$-model and assuming $\gamma=5/3$ (for a monoatomic gas), we find:

[TABLE]

and therefore:

[TABLE]