AGN radiation imprints on the circumgalactic medium of massive galaxies

Active Galactic Nuclei (AGN) in cosmological simulations generate explosive feedback that regulates star formation in massive galaxies, modifying the gas phase structure out to large distances. Here, we explore the direct effects that AGN radiation has on gas heating and cooling within one high-resolution $z=3$ dark matter halo as massive as a quasar host ($M_{\rm h}=$10$^{\rm 12.5}$M$_{\rm\odot}$), run without AGN feedback. We assume AGN radiation to impact the circumgalactic medium (CGM) anisotropically, within a bi-cone of angle $\alpha$. We find that even a relatively weak AGN (black hole mass $M_{\rm\bullet}=10^{\rm 8}$M$_{\rm\odot}$ with an Eddington ratio $\lambda=0.1$) can significantly lower the fraction of halo gas that is catastrophically cooling compared to the case of gas photoionized only by the ultraviolet background (UVB). Varying $M_{\rm\bullet}$, $\lambda$ and $\alpha$, we study their effects on observables. A 10$^{\rm 9}$M$_{\rm\odot}$ AGN with $\lambda=0.1$ and $\alpha\approxeq60^{^{\rm o}}$ reproduces the average surface brightness (SB) profiles of Ly$\alpha$, HeII and CIV, and results in a covering fraction of optically thick absorbers within observational estimates. The simulated SB$_{\rm CIV}$ profile is steeper than observed, indicating that not enough metals are pushed beyond the very inner CGM. For this combination of parameters, the CGM mass catastrophically cooling is reduced by half with respect to the UVB-only case, with roughly same mass out of hydrostatic equilibrium heating up and cooling down, hinting to the importance of self-regulation around AGNs. This study showcases how CGM observations can constrain not only the properties of the CGM itself, but also those of the AGN engine.


INTRODUCTION
Super massive black holes (SMBHs, ≥10 5 M ⊙ ) can inhabit all types of galaxies, from dwarfs (e.g.Reines et al. 2013) to bright cluster galaxies (e.g.Gear et al. 1985).Their strong gravitational field, very spatially localized, has only been directly measured in our own galaxy (e.g.Genzel et al. 1997), or in the very nearby universe (e.g.Event Horizon Telescope Collaboration et al. 2019).At large distances (redshifts) their presence is typically inferred from the strong radiation field emitted from their vicinity during periods of accretion, hence the name AGN (Lynden-Bell 1969).
AGNs are distinguished from other astrophysical radiation sources by: extremely high bolometric luminosities (up to  bol ≲10 49 erg s −1 ), a very strong evolution of their luminosity function (e.g.Shen et al. 2020), and SEDs that are detectable across all wavelengths, from radio to gamma rays (e.g.Bianchi et al. 2022).
★ E-mail: obreja@usm.lmu.deTheir extreme luminosity makes them valuable cosmic beacons to find and explore the properties of diffuse gas even in the highest redshift universe, both in emission and in absorption against their light (e.g.Farina et al. 2019;D'Odorico et al. 2023).The combination between the AGN's high power and the compactness of their emitting regions implies very high energy densities, bound to have an effect on their host galaxies (e.g.Bower et al. 2008).Current observational signposts of AGN effects on their host galaxies include: highly collimated jets on scales up to hundreds of kpc (e.g.Waggett et al. 1977;Biretta et al. 1983), outflows/(ultra)fast winds on scales of the order of kpc to tens of kpc (e.g.Pounds et al. 2003;Ganguly et al. 2007;Ajello et al. 2021), high energy bubbles and cavities in the circumgalactic medium anchored on the galaxy's center (e.g.Su et al. 2010;Predehl et al. 2020).On the other hand, scaling relation between SMBH masses and galaxy properties like bulge mass (Magorrian et al. 1998), central stellar velocity dispersion (Ferrarese & Merritt 2000;Gebhardt et al. 2000), and hot circumgalactic medium temperature and X-ray luminosity (Gaspari et al. 2019;Lakhchaura et al. 2019) provide further contextual evidence for such effects.All these observations taken together suggest that the energy and momentum generated by the accretion heats, ionizes and pushes outwards the gas surrounding the AGN, thus regulating the SMBH growth and star formation in the galaxy (e.g.Begelman & Nath 2005;Fabian 2012).
One of the strongest indirect evidences in favor of an active role of AGNs in galaxy evolution comes from semi-analytical models, which have shown that the dependency of stellar mass on dark matter halo mass is best described by a double power-law (e.g.Moster et al. 2018;Behroozi et al. 2019).This translates into a particular halo mass (≈10 12 M ⊙ ) for which the efficiency of star formation is maximum.In both semi-analytical studies and cosmological simulations of galaxy formation, this power-law dependency of galaxy stellar mass on dark matter halo mass can only be reproduced if the hydrodynamics equations take into account not only the gas cooling and star formation, but also the back-reaction of gas under energetic phenomena like supernovae at low halo masses (SNe, e.g.Stinson et al. 2006;Dalla Vecchia & Schaye 2008;Oppenheimer & Davé 2008;Keller et al. 2014) and AGN at high halo masses (e.g.Springel et al. 2005;Sĳacki et al. 2007;Booth & Schaye 2009;Debuhr et al. 2011;Costa et al. 2014;Steinborn et al. 2015;Weinberger et al. 2017).The impact of these energetic phenomena on the formation and evolution of galaxies are generally referred to as feedback processes.In this framework, simulations, and especially cosmological simulations that also trace the formation of galaxies within the cosmic web, are a powerful tool to follow and understand the complex astrophysics of structure build-up through cosmic time.
In cosmological simulations, AGN feedback is generally divided in a quasar (radiative) and a radio (jet) mode, none of which can be modeled from first principles due to resolution limitations.Quasar mode feedback (regime of high accretion rates) has been implemented in two distinct ways: as a thermal energy transfer from the accreting SMBH to the nearby gas (e.g.Springel et al. 2005;Di Matteo et al. 2005;Booth & Schaye 2009;Steinborn et al. 2015), and as a momentum (kinetic) transfer from the AGN photons to the nearby gas (e.g.radiation driven winds Costa et al. 2014Costa et al. , 2018)).In these radiative feedback schemes, the energy and/or momentum imparted to the neighboring gas is proportional to the SMBH's bolometric luminosity, which is assumed to scale linearly with the accretion rate.The radio mode feedback is observationally based on the relativistic jets associated with low accreting AGN (e.g.Blandford et al. 2019;Kondapally et al. 2023), and it is thought to be responsible for the large high energy bubbles observed around galaxies (e.g.Guo & Mathews 2012;Yang et al. 2022).The implementations of radio mode feedback in cosmological simulations is severely limited by resolution, as the modeling of the jets requires general relativistic magnetohydrodynamical simulations with very high resolutions (e.g.Ressler et al. 2017).For this reason, this feedback mode is implemented also as a momentum transfer in a very approximate manner, and often referred to as kinetic feedback (e.g.Sĳacki et al. 2007;Weinberger et al. 2017;Davé et al. 2019).A third, much less explored feedback mode, dubbed electromagnetic (EM) feedback by Vogelsberger et al. (2013) aims to quantify the direct photoionization/photoheating of gas by AGN radiation in simple toy models or isolated galaxy simulations (e.g.Ciotti & Ostriker 1997, 2007;Sazonov et al. 2005;Choi et al. 2012;Gnedin & Hollon 2012;Xie et al. 2017).In cosmological simulations, Kim et al. (2011) implemented EM feedback by solving the radiation transfer on-the-fly assuming a 2 keV monochromatic spectrum that scales with the AGN's bolometric luminosity, while Vogelsberger et al. (2013) used a universal, time independent AGN spectrum, to compute the photoionization and photoheating rates of hydrogen and helium, under the assumption that gas is optically thin to the AGN radiation.
While the various flavours of AGN feedback models reproduce fairly well galaxy properties, they produce diverging results for the amount of baryons locked in virialized structures.In particular, current AGN feedback models reduce the baryon fractions in massive halos (e.g.Wright et al. 2020;Sorini et al. 2022) in comparison with the cosmic mean (e.g.0.157, Planck Collaboration et al. 2014).This prediction might be in tension with X-ray observations at low redshifts (e.g.Chadayammuri et al. (2022), but see also Comparat et al. (2022) for a different conclusion).Mapping the hot CGM in galaxies is currently becoming a reality (e.g.Nicastro et al. 2023), and the development of new instruments will, hopefully, allow measuring the baryon fraction across galaxy populations (e.g.Schellenberger et al. 2023).Another observational estimate in disagreement with current simulations is the non-thermal to thermal pressure ratio in galaxy clusters, whose value is significantly smaller than what such simulations predict (e.g.Sayers et al. 2021).This indicates that current AGN implementations introduce too much turbulence in the interstellar/circumgalactic medium.All in all, CGM observations across the full electromagnetic spectrum are needed to constrain AGN feedback models, and galaxy formation models in general.
A powerful method to estimate the properties of the multi-phase CGM, is the study of absorption features in the spectra of bright background sources, i.e. quasars (e.g.see review by Tumlinson et al. 2017).In this method, one can statistically probe various impact parameters within the CGM regions of foreground galaxies.The subsequent modeling of the absorption features can put constraints on the nature of the absorber, and therefore on the gas properties in the CGM of the foreground system like, e.g., the ionization state and metallicity (e.g.Prochaska et al. 2013b;Werk et al. 2014;Lau et al. 2016).Studies of quasars' CGM in absorption at  = 2 − 3 have found evidence for large reservoirs (>10 10 M ⊙ ) of cool gas with relatively high metallicities (> 0.1  ⊙ ) (e.g.Prochaska et al. 2013b;Lau et al. 2016), and with kinematics mostly in agreement with those of the hosting dark matter halos (e.g.Prochaska et al. 2014;Lau et al. 2018).Importantly, the distribution of optically thick H i systems around quasars is found to be highly anisotropic, with the transverse direction preferred over the line of sight (Hennawi & Prochaska 2007).This fact suggests that quasars illuminate their surrounding anisotropically.
Besides absorption studies, a complementary and more direct method to constrain the effects of AGN radiation on the CGM is provided by line emission maps.Fortunately, the high redshift CGM is now accessible with sensitive integral field unit (IFU) spectrographs like the Multi-Unit Spectroscopic Explorer on the Very Large Telescope (MUSE, Bacon et al. 2010), or the Keck Cosmic Web Imager on Keck (KCWI, Morrissey et al. 2018), which can directly probe rest-frame UV emission lines from the high-z cool CGM.The usually strongest emission line is the hydrogen Lyman  (Ly) transition at 1215.67Å restframe.In particular, MUSE and KCWI have facilitated an explosion in the number of surveyed objects for CGM emission at high-, resulting in the detection of extended Ly glows around a few hundred quasars (e.g.Borisova et al. 2016;Arrigoni Battaia et al. 2019;Cai et al. 2019;Farina et al. 2019;Fossati et al. 2021).These Ly nebulae around quasars have been long ago predicted (e.g.Rees 1988), but the level of emission is so faint that only in the last decade they started to be routinely observed.The largest homogeneous sample of Ly glows around quasars is the QSOMUSEUM survey (Arrigoni Battaia et al. 2019) which targeted 61 quasar fields at  ∼ 3.This survey found a wide range of nebula morphologies and sizes, from an asymmetric 'enormous Lyman alpha nebula' (ELAN) spanning ∼300 kpc to small symmetric ones that barely reach the CGM region (Arrigoni Battaia et al. 2023).These observations are a unique probe into the CGM of quasars, as Ly traces cool gas ( ∼ 10 4 K), which can subsequently fuel star formation and AGN activity on galaxy scales.A detailed analysis of the ELAN in the QSOMUSEUM survey showed that the cool gas is likely inspiraling towards the central galaxy, and has velocity dispersions in agreement with pure gravitational motions within a ∼10 12.5 M ⊙ dark matter halo (Arrigoni Battaia et al. 2018).How this cool and likely dense (e.g.Hennawi et al. 2015;Arrigoni Battaia et al. 2015b) gas manages to survive in the hot (virial temperatures of the order of 10 6.5 K) and harsh environment of quasars is an active topic of research (e.g.Bennett & Sĳacki 2020;Costa et al. 2022;Gronke et al. 2022).
In this work, we quantify, in post-processing the intrinsic and observable effects induced by photoionization and photoheating from an accreating SMBH on the CGM of a high resolution cosmological  = 3 quasar host halo.The paper is structured as follows.In Sections 2 and 3 we describe the simulation and the spectra of the radiation sources considered.Section 4 explains how we construct the photoionization models, while Section 5 quantifies the intrinsic effects of an anisotropic AGN radiation field on CGM heating and cooling rates.Our comparisons with CGM observations in absorption and emission are given in Section 6 and 7, while in Section 8 we discuss the implications for on-going and future observational efforts to map the CGM of high redshift quasars.Finally, Section 9 summarizes our results.This study is a first step in the development of an updated AGN EM feedback model for cosmological simulations.

THE SIMULATED HALO
We use a zoom-in cosmological simulation of a quasar host halo at  ∼ 3 to study in post-processing the radiation effects on the CGM of a central supermassive black hole.This simulation has been run without AGN feedback, using the exact same galaxy model as NIHAO (all free parameters are set to the values in Wang et al. 2015), in- cluding the same cold dark matter cosmology (Planck Collaboration et al. 2014).The code used to simulate this galaxy, named g3.16e12, is the updated version (Gasoline2, Wadsley et al. 2017) of the Nbody smoothed particle hydrodynamics code Gasoline (Wadsley et al. 2004).Stellar particles, representing single stellar populations, are formed stochastically from dense ( H > 10.3 cm −3 ) and cold ( <15000 K) gas such that a Kennicutt-Schmidt type of relation is recovered.Gas is metal enriched by supernovae type I and type II (Raiteri et al. 1996), computed assuming the initial stellar mass function of Chabrier (2003), and the metal diffusion and mixing follows the prescriptions of Shen et al. (2010) and Wadsley et al. (2008).Non-equilibrium cooling of hydrogen and helium is computed onthe-fly (Shen et al. 2010), while for the metal line cooling the gas is assumed in photoionization equilibrium with the UVB of Haardt & Madau (2012).Stellar particles further impact their surrounding gas through two types of feedback: i) heating prior to SNe II phase also known as 'early stellar feedback' (Stinson et al. 2013), and ii) blast-waves from the SNe II events (Stinson et al. 2006).To prevent the rapid dissipation of the SNe II thermal energy transferred to the gas (Stinson et al. 2006), the cooling of gas particles affected by this feedback is delayed by 40 Myr.The zoom-in region of this simulation has been run with mass resolutions of 6.5×10 5 M ⊙ and 3.6×10 4 M ⊙ , and with gravitational softening lengths of 168 pc and 72 pc for the dark matter and gas particles, respectively.This resolution allows us to follow the gas distribution up to densities of ∼10 3 cm −3 .and hot ( ≥10 5 K, orange) gas, and all gas (purple).Center: number density of hydrogen as a function of radius for the cool and hot gas phases.Bottom: metallicity as a function of radius for the cool and hot gas phases.In the middle and bottom panels, the points and their errors give the mean and standard deviation, the solid lines give the medians, and the shaded regions enclose from the 16 th to the 84 th percentiles (darker shade), and from the 2 nd to the 98 th percentiles (lighter shade).We also show the collapsed density distribution functions of hydrogen number density and metallicity in the right panels.The vertical grey line marks the virial radius  200 .In all panels, we exclude the ISM region  < 20 kpc, where  is the 3D radius.
Figure 1 shows the dark matter (top left), stellar (top right), gas mass (bottom left), and gas-phase metal mass (bottom right) distributions in a 4 ×  200 -box centered on the quasar host halo at  = 3.The dashed white circle marks  200 = 134kpc, defined using the Amiga Halo Finder (Knollmann & Knebe 2009).The gas column density panel (bottom left) shows nicely the cosmic-web filaments feeding this halo, and the in-falling galaxy sitting roughly at  200 from the center of the main halo.The stellar light surface brightness (top right) more clearly displays the on-going formation process, while many of the dark matter substructures visible just outside of the virial radius are clearly devoid of stars (top left).At this redshift ( = 3), the dark matter halo mass is  h =2.84×10 12 M ⊙ , while the virial star and gas masses are  ★ =1.09×10 11 M ⊙ and  gas =3.25×10 11 M ⊙ , respectively.This implies an integrated baryon conversion efficiency  ≡  ★ /( ★ +  gas )=0.25, which is within the 1 scatter of the abundance matching predictions for this galaxy's dark matter halo mass and redshift: ( h ,  = 3) = 0.19 +0.08  −0.06 (Moster et al. 2018).Figure 2 shows the mass weighted phase space for the selected box (top), as well as for what is typically referred to as the CGM (center), and ISM (bottom) regions1 .The left and right panels show the phase space weighted by the hydrogen and median metallicity, respectively.Both the CGM and the ISM regions show two distinct blobs.In the bottom left panel, the lower blob is the cool/cold ISM, while the upper one is made by the ISM gas heated by SNe.The ISM gas heated by SNe also carries a significant amount of the metals, as shown in the bottom right panel.The two blobs in the central left panel represent the hot corona (upper blob) and the "cooling flows" (lower blob) in the CGM.The other important feature in phase-space is visible as a straight line ( ∝    with  > 0) in the region of low temperature and low density of the full box upper panels, and represents gas of the IGM.
Looking at the phase diagram of CGM gas, we see that  ∼10 5 K can separate relatively cleanly the cool from the hot phase.Therefore, in the two lower panels of Figure 3 we show how gas density and metallicity varies with the radius for the two phases.The hydrogen number density for the hot gas drops smoothly with  (orange curve and points in the central panel), with an almost log normal distribution in  H when considering all the region with 20 kpc ≈ 0.15 vir <  < 300 kpc.Contrary, the radial variation of  H for the cool gas also drops with , but it is significantly more wiggly as it also incorporates the ISM of nearby smaller galaxies.At fixed  the distribution of  H for the cool phase is significantly more wide than for the hot one, and has the median at larger  H values.On the other hand, the radial distributions of metallicities (bottom panel) for the two phases are very similar.In the upper panel of Figure 3 we also show how the mass of dark matter, hot and cool gas, and total gas is distributed with .
To summarize, g3.16e12 at  = 3 is a high resolution galaxy (resolved densities in the ISM of ∼10 3 cm −3 ), with a halo mass within the observed range for quasar hosts (e.g.Lau et al. 2016;Petter et al. 2023), and with a stellar mass compatible with the constraints from abundance matching (e.g.Moster et al. 2018).The simulation has been run without on-the-fly AGN feedback.In the next sections, we post-process the simulation to study the effects of AGN radiation (from a central SMBH with mass compatible with the halo mass) on its CGM's intrinsic properties, as well as on observational constraints like column densities of various ions and surface brightness in various emission lines.

THEORETICAL MODEL SPECTRA
In a galaxy there are a few different types of local ionizing sources (e.g.young massive stars, AGN, X-ray binaries, hot gas emitting Bremsstrahlung), and few cosmological galaxy simulation studies have tried to incorporate at least the heating/ionizing effect of some of them (e.g.Kannan et al. 2014Kannan et al. , 2016;;Obreja et al. 2019;Hopkins et al. 2020;Costa et al. 2022).While most works have focused on the ionizing continuum from stellar sources, few studies have looked at the effects of ionizing radiation from AGN (e.g.Kim et al. 2011;Vogelsberger et al. 2013), but none studied the EM feedback alone, making it impossible to disentangle its effects from those induced by the widely used thermal and kinetic AGN feedback models.
AGN have been classified in a myriad of sub-classes based on features observed in their spectra and/or their wavelength of first detection.In time, only three properties emerged as fundamental: the radiative efficiency, the presence/absence of a jet, and the galaxy host environment (Antonucci & Miller 1985;Antonucci 1993;Bianchi et al. 2022).The radiative efficiency measures how efficient the accretion is, and it is quantified by the ratio between the observed and Eddington luminosity  = / Edd , where  Edd = 1.3 × 10 46 (/10 8  ⊙ )erg s −1 is the maximum luminosity of a mass  accreting isotropically.The third property -galaxy host environment -was the key one for what it is now known as the Unification Model (Antonucci 1983).The idea behind this class of models is that all types of AGN have the same nuclear engine, and that an obscuring axisymmetric structure (a dusty torus) around the SMBH defines how the AGN spectrum appears to a distant observer.The presence of such an optically thick torus, on parsec scales and aligned with the symmetry axis, obscures the inner broad line region (typical in the spectra of type I AGN) when the observer line-of-sight (LOS) makes a sufficiently large angle with the symmetry axis, while leaving unaffected the outer narrow line region (typical in the spectra of type II AGN).As the multi-wavelength coverage of AGN spectra has grown, it became clear that this class of models is too simplistic (e.g.Combes 2021), and that different absorbers on different scales are needed to reproduce well high resolution high sensitivity AGN spectra (e.g.Gaspari et al. 2017).Notwithstanding these complications on small scales, the radiation from an AGN is currently though to illuminate anisotropically the surrounding CGM with its so-called ionization cones.

AGN spectra
For the spectral energy distribution (SED) of accreting SMBHs, we use the optxagnf 2 model (Done et al. 2012).This theoretical SED model is made of three components: • a pseudo-thermal accretion disk (Novikov & Thorne 1973), • a non-thermal power-law component at high energies produced by Compton up-scattering from an optically thin and hot medium, • a soft X-ray excess component produced by Compton upscattering from an optically thick, lower temperature medium.This model is energetically self consistent, in the sense that the primary and only source of energy is the constant mass accretion rate  accretion from a Novikov-Thorne disk onto a black hole (BH).The schematic view of the AGN emission using such a model is shown in Figure 5 of Done et al. (2012).This model has nine parameters, some of which are partially degenerate: (i) black hole mass  • , (ii) AGN luminosity in units of Eddington luminosity / Edd , (iii) black hole spin parameter , (iv) the coronal radius  cor in units of the gravitational radius  g =   BH / 2 , which gives the innermost edge of the disk that is visible, (v) the outer radius of the disk  out in units of  g , (vi) the temperature  of the optically thick gas producing the soft X-ray excess, (vii) the optical depth  of the optically thick gas, 2 The code that can generate this type of spectrum can be found as a local model in the publicly available XSPEC spectral fitting package (Arnaud 1996).
10 3 10 2 10 1 10 0 10 1 10 2 10 (viii) the fraction of the corona energy emitted in the power-law component,  PL (the fraction of the corona energy emitted as the soft X-ray excess is 1 −  PL ) (ix) the high energy power-law index Γ.
Jin et al. ( 2012) used the optxagnf model to fit a sample of 51 low redshift Seyfert 1 AGNs.In this sample, 12 objects are classified as Narrow Lines Seyfert 1 (NLS1), while the rest are Broad Lines Seyfert 1 (BLS1).The BLS1 objects are the closest we can get to the intrinsic SED of AGNs, which is precisely what we need in the simulations.Hence, we choose the SED parameters based on the corresponding distributions for the BLS1 sample.To construct our template SED models we fixed the BH spin to  = 0 and the outer disk radius to  out = 10 4  g , following Jin et al. (2012).For the remaining parameters that our target zoom-in cosmological simulations can not resolve, we used the median values over the BLS1 sample.A few of the AGN in this sample have coronal radii at the upper boundary for  out ,  out = 100 g .For this reason we excluded these AGNs when computing the median  With all these parameters fixed, our template AGN spectra depend only on the black hole mass  • and on the Eddington accretion ratio  = / Edd .
The upper panel of Figure 4 gives the spectra for all cases of radiation sources we consider in this study.Notably, the AGN spectra (blue and orange curves) provide significantly more ionizing photons ( > 1Ryd) than the UVB even at large distances from the quasar (the AGN spectra in the figure correspond to a distance of 1 Mpc).An AGN also emits large number of photons in UV emission lines.Fossati et al. (2021) has shown that the UV line emission spectra of AGN are crucial to explaining the observed line ratios from the CGM of high redshift quasars with reasonable gas metallicity values.For this reason, we added on top of the AGN continuum spectra shown in the upper panel, the stacked emission line spectrum of 53 luminous quasars at  ∼ 2.4 from Lusso et al. (2015).The lower panel of Figure 4 shows the resulting AGN spectra in the UV region marked by the dotted gray rectangle in the upper panel.To obtain only the observed stacked emission lines, we subtracted from the published spectrum the continuum fitted by Lusso et al. (2015):   ∝  −Γ with Γ = 1.70 ± 0.61 at  ≤ 912Å.This spectral slope is, within its uncertainties, compatible with the one we assumed (Γ = 1.81).

UVB spectrum
Our cosmological simulation, like most such simulations in the literature, has been run under the assumption that gas is photoionized by an isotropic and time dependent UVB.Therefore, our default cloudy models use the recent UVB spectra of Khaire & Srianand (2019).To create synthesis models of the extragalactic background light, extending from far infra-red to -rays, these authors used up-to-date values for the cosmic SFR density, dust attenuation, quasar emissivity, and neutral hydrogen distribution in the IGM.Khaire & Srianand (2019) do a particularly detailed treatment of the extreme ultraviolet background, which determines the ionization and temperature of the IGM across time.One practical advantage of using these particular UVB models over the recent alternatives (e.g.Puchwein et al. 2019;Faucher-Giguère 2020), is that they are already incorporated in the latest version of the spectral synthesis code cloudy (version C17.03, Ferland et al. 2017), which we use to create the ancillary data (e.g.cooling and heating tables) needed for cosmological galaxy formation simulations.It is also important to mention that UVB models are generally converged at redshifts lower than the end of reionization (e.g.Khaire & Srianand 2019), and for this reason the mismatch between predictions using our choice of UVB in this work (Khaire & Srianand 2019) and the UVB embedded in the simulations (Haardt & Madau 2012) should be minimal.
In the following sections we describe both the intrinsic and the observable effects that these types of radiation sources have on our simulated quasar host halo, in post-processing.

POST-PROCESSING WITH CLOUDY
In order to create mock observables from this simulation we run large grids of photoionization models with the version C17.03 of cloudy (Ferland et al. 2017), covering the multi-dimensional phase space of the gas in a 4× 200 box centered on the galaxy g3.16e12.These models are plane-parallel models for fixed: total hydrogen density  H , gas temperature , metallicity , and incident radiation field   .The control models are run with the sum of the cosmic microwave background (CMB) and the UVB of Khaire & Srianand (2019) at  = 3.The photoionization models with the AGN spectra of Figure 4 as incident spectra also include the CMB (but not the UVB), and span a range in normalization corresponding to a minimum distance from the quasar of 1 kpc and a maximum distance of 2

√
3 200 (the distance from the center of the simulated galaxy to the farthest corner of the selected box region).We do not include molecules, dust and cosmic rays, because the cosmological simulation does not follow any of them.To take into account the self-shielding of the gas from the incident radiation field, we run the photoionization models with a column density stopping criteria, depending on the local Jeans scale (e.g.Schaye 2001), following the set-up used by Ploeckinger & Schaye (2020) (see section 2.1 of their paper).Other approximations for the local column density, e.g.Sobolev-like based on the local density (Krumholz & Gnedin 2011), would require adding extradimensions to the already many gas properties we need to account for.Instead, using a stopping column density based on the Jeans lenght means that the column depends only on temperature and hydrogen number density.
The outputs of these photoionization models include: heating and cooling rates, electron densities, ionization fractions for a few species for which column densities can be estimated from observations, and line emissivities for a list of optical and UV lines that have been or could be observed in QSO host halos at  ∼ 3.All these properties are saved for the last cloudy zone, and all are functions of ( H , , ,   ).An example of a cloudy script used is explained in detail in Appendix A. Figure 4.
The main reason for including AGN feedback in simulations of massive galaxies is to mitigate the so-called 'overcooling problem' (Larson 1974;White & Rees 1978;Dekel & Silk 1986;Navarro & Benz 1991).Thus, we first look at the impact AGN radiation has on the net cooling rate of the gas.From Figure 2, the hot gas reservoir of the CGM is metal enriched and has densities of the order of 10 −3 cm −3 .In Figure 5 we show how the absolute net cooling rate |Λ − |, where Λ and  are the volumetric cooling and heating rates, for a slab of gas of metallicity 0.1 ⊙ and density 10 −2 cm −3 varies with the temperature.The black curve in the figure corresponds to the case where the gas is only irradiated by the UVB, while the colored curves show |Λ − | as a function of  for gas irradiated by a 10 9 M ⊙ SMBH accreting at ten percent the Eddington rate ( = 0.1) for various distances between the slab of gas and the radiation source (indicated in the legend on the right side).The curves are solid in the region where the gas is cooling (Λ > ) and dashed where it is heating (Λ < ).The most notable thing in this figure is how the equilibrium temperature ( for which Λ = ) moves to higher and higher values as the strength of the incident radiation field increases.While for the UVB case  eq ≈10 4.2 K, gas irradiated by a central quasar and sitting at only 5 kpc (red curve) has  eq ≈10 6 K. Based on this figure and the phase space distribution of CGM gas in Figure 2, we can infer that a significant fraction of the CGM gas will not be cooling, but instead will be heating in the presence of a central AGN radiation source.

GEOMETRIC EFFECTS
As already stated at the beginning of Section 3, the radiation emitted by an AGN does not propagate on halo scales isotropically, but rather within ionization cones, because of obscuring media on small scales.Flavours of the AGN unification model, though very approximate, are quite useful for sub-grid models for galaxy formation simulations, which can not yet resolve the parsec and sub-parsec scales relevant for the SMBH feeding.For example, Fritz et al. (2006) modeled the dust torus of the SMBH as a flared disk (Efstathiou & Rowan-Robinson 1995), and found that roughly half of their type I AGN and ∼70% of the type II were best fitted with a torus opening angle  ≈140 o , while the rest of the sample required  ≈100 o .This means that ionization cone opening angles of AGN  = 180 o −  should be mostly restricted to 40 o ≤  ≤ 80 o .
We concentrate in this work on how the ionization cone opening angle affects not only the intrinsic properties of the CGM gas, but also its observables.We define ionization cones by their opening angle 3 and explore the range between 180 • (isotropic case) and 0 • (no AGN radiation).For each  we choose  random directions passing through the center of the halo to place the cones, and for each of these  directions we place the line of sight (LOS) randomly within the ionization cones.The choice of  is different for the various exercises we do in the following sections.
First, we investigate how much of the CGM would be catastrophically cooling as a function of opening angle.While the value of the cooling time in itself might not be very informative for  cool > 0, comparing it with the free-fall time  ff can provide a good idea of how much of the gas is in (quasi)hydrostatic equilibrium, and how much is catastrophically cooling.The free-fall time is only a function of the total average system's density and thus varies with the position in the dark matter halo: where  is the gravitational constant, and   (< ) is the average mass (dark matter + baryons) density within a sphere of radius :   (< ) = 3 dyn (< )/4 3 .Using the definition of Equation 2 for the free-fall time, Figure 6 quantifies exactly how much of the CGM gas is catastrophically cooling, heating or is in equilibrium, as a function of the quasar spectrum and ionization cone opening angle.The errors on the individual points are the standard deviations of the mass fractions over the =200 random directions chosen to place the ionization cones.In this figure  = 0 corresponds to the case where all gas sees only the UVB, while increasing values of  correspond to the case where only the gas within the double cone of opening angle  is irradiated by the quasar while the rest of it is irradiated by the UVB.The case =180 o means that all gas is directly irradiated by the quasar.Figure 6 illustrates nicely the 'overcooling' problem, 37% of the total CGM gas ( ∼2.43×10 11 M ⊙ ) experiencing catastrophic cooling if the only source of radiation is the UVB.As the ionization cone opening angle increases, an increasingly larger fraction of the CGM is heating (salmon curves), while the cooling fraction (black curves) decreases accordingly reaching almost zero (  cooling ∼0.7%) for the extreme =180 o and the stronger incident spectra.Only the spectra with  • =10 8 M ⊙ and  = 0.1 has a significant  cooling ∼8% for =180 o .On the other hand, the fraction of CGM gas in (quasi)hydrostatic equilibrium (blue curves) is only slightly increasing with , and varies little with the strength of the quasar spectrum.
Figure 7 shows the CGM gas in the density-temperature phasespace for three different opening angles.The color coding is the same as in Figure 6, and the SMBH spectra is for  • =10 9 M ⊙ and  = 0.1.In the left panel corresponding to the UVB only case, we see that most of the gas with temperatures around 10 4 K is cooling and, therefore, it is a reservoir for the ISM.The same is true for some of the gas in the hot CGM.For =60 o shown in the central panel (value within the observational limits, e.g.Fritz et al. 2006), the fraction of CGM gas which is a reservoir for the ISM is roughly halved, while for the extreme =180 o case, virtually all the hot CGM gas is in equilibrium, while the cool CGM is heating under the effect of the central accreting SMBH.Therefore, AGN with large opening angles can completely stall or at least delay the star-formation process as the cool CGM is ultimately its fuel.

CGM OBSERVABLES IN ABSORPTION
For our target galaxies ( ∼ 3 quasar hosts) the best observational data-set in absorption we can compare our simulation with is the one of the project "Quasars probing quasars" (QPQ) introduced by Hennawi et al. (2006).This work pioneered the method of using projected quasar pairs, that are not physically associated with each other but are close on-the-sky, to study quasar environments.Along most sight-lines of the QPQ sample, the measured Lyman- equivalent width is low  ≈1-2Å, meaning that the H i column density can not be well estimated because the systems are on the flat part of the curve-of-growth.As a consequence most observations of the QPQ sample can only place limits on  H i (Prochaska et al. 2013b).Lau et al. (2016) used only a very small sub-sample of QPQ systems with good measurements of  H i to also try quantifying the ionization parameter (representing the ratio between the number of ionizing photons and number of hydrogen atoms) at the absorber location.The foreground quasars in this sub-sample have an average redshift <  >= 2.66, a minimum of 2.00, and a maximum of 4.11.This latter work also provides a formula to approximate the mass of cool gas as a function of projected distance from the quasar.Using their formula we would expect to have ∼1.3×10 11M ⊙ within the  200 =134 kpc of our simulation, which is larger than our cool ( <105 K) CGM mass of ∼0.9×10 11 M ⊙ .As mass fractions, the cool CGM in our simulated galaxy constitutes ∼ 1/5 of the baryons in the halo vs the ∼ 1/3 estimated by Lau et al. (2016) from their 14 quasar pairs.In the following we provide a first test of how the simulation post-processed with AGN radiation compares to these observations.First, the top row panels of Figure 8 show the comparison between observed and simulated H i column density as a function of projected distance from the quasar.To compute the percentiles for  H i for each particular ionization cone opening angle, we proceeded as follows.We generate 10 uniformly distributed directions over the unit sphere, and consider them to give the directions of the ionization cones.Since a quasar is an AGN viewed through its ionization cone, for each  the LOS4 is not centered on the ionization cone axis, but it is chosen randomly among those within the cone 5 .Next, for each LOS we approximate the maps  H i (, ) by: where  j,H i (, ) is the column density contribution of particle  with projected position within the pixel (, ) and position along the LOS, −500kpc << 500kpc.In Equation 3,  is the gas particle mass, X is the hydrogen mass fraction, X H i is the ionization fraction, and Δ  is the 'size' of the particle approximated as Δ  ≡ (  /  ) 1/3 with   the mass density of the particle.To limit the computational time while still having pencil-beam-like information, the maps are done assuming a pixel size of Δ=1.58 kpc (∼0.2 arcsec assuming the halo at  = 3).While the hydrogen mass fraction X is given by the simulation, the ionization fraction X H i ( H , , ,   ) is interpolated from the pre-computed cloudy tables.For particles within the ionization cone of the quasar,   has one of the corresponding spectral shapes shown in Figure 4 (blue or orange curves) and a normalization depending on the radial distance to the center of the halo.For particles outside of the ionization cone,   is given by the UVB spectrum at  = 3 also shown in Figure 4. Finally, we compute the percentiles in radial bins using the  H i (, ) of all 10 LOS together.A more faithful comparison between simulations and observations would require to construct spectra along background sight-lines using a post-processing tool like e.g.trident (Hummels et al. 2017), but this is outside the scope of this paper, and therefore we leave such analysis for a future work.
The radial profiles of  H i in Figure 8 show the simulation to be compatible with the observations if the only source of photoionization is the UVB (all observational data points are between the 2 nd and 98 th percentiles of the  = 0 o case given by the dashed grey curves).This result is in line with previous predictions from simulations run only with stellar feedback (e.g.Faucher-Giguère et al. 2016).As we include ever stronger AGN radiation fields (left to right panels), the simulation is compatible with the observations only for small open-ing angles: the regions between the 2 nd and 98 th percentiles of the  = 60 o case (dashed blue curves) still include all observational data points.The same percentiles for the  = 180 o case (dashed green curves) include all of the observational data points only for the weakest AGN model (left panel).Both these effects -with the strength of the radiation field, and with the ionization cone opening angle -are to be expected, as the ionization parameter  increases with both   and  for fixed number density  H .
It is important to keep in mind that we look at one simulation only, while the observational data points represent the CGM of different  in observations as coming from more massive dark matter hosts (10 12.5 M ⊙ <  h ≤10 13 M ⊙ ) or from the ISM of satellites (e.g. the bump in the higher percentiles simulated curves at  ∼130 kpc is due to the satellite galaxy just entering the virial radius, as can be seen in Figure 1).For the first possibility, it is instructive to look at the profile of total hydrogen column density  H in comparison with the  H ≡  H i /X H i reconstructed from the observations, which we show in Figure 9.While most observational estimates are in agreement with the simulation, a few data points are clearly outside the simulation range.The most natural explanation for these data points is that the corresponding absorbers are coming from dark matter halos of different masses.Given that the total gas density in a dark matter halo decreases with radius, we expect the  H () of the simulation to also decrease with the impact parameter, as seen in the figure.On the other hand, the total hydrogen content of a halo is proportional to the dark matter mass.Therefore, a more/less massive halo would shift the simulated profile up/down.In each of our models, we can also compute the ionization parameter   for each particle  in the simulation: where  is the speed of light, and Φ  () is the ionizing photon flux ( >13.6 eV) at the particle position (see the cloudy documentation), which depends on the AGN SED model and distance to the AGN for particles within the ionization cone, and only on the particle's hydrogen number density  j,H for particles outside (as the UVB just gives a constant Φ() for a given redshift).To compute an average ionization parameter  we weight each particle by its column  j,H i (, ) to construct similar  (, ) maps: The panels in the second row of Figure 8 show the comparison between the radial profiles of  and the values estimated by Lau et al. (2016) for the absorbers in the 14 foreground quasars of their pairs.In their calculation, it is assumed, as in our work, that the ionization structure of the gas is set by photoionization.Specifically, Lau et al. (2016) modeled each absorber as a plane-parallel slab with its corresponding  HI value and a fixed  H = 0.1 cm −3 .The ionization parameter  is determined by varying the metallicity and the input radiation field from the UVB of Haardt & Madau (2012) to higher amplitude assuming a EUV power-law spectrum   ∝  −1.57until the results converge on the available multiple ionization states of individual elements.Same as for the top row, the simulations are in good agreement with the observations if the only photoionization source is the UVB, or if the AGN's ionization cone opening angle is small enough.To make the comparison between observations and simulations clearer, the lower panels of Figure 8 show the 2D log( H i )-log() distribution of sightlines in a logarithmic scale (colored map with red marking the maximum) for three different angles .In these panels it is clear that the observational data favor the AGN models with low accretion rate  = 0.1 (left and center) and low opening angle  = 60 o .On the other hand the UVB model alone (the panels with  = 0) results in a very marked anti-correlation between log() and log( H i ), which, however, is not spread enough in log() to cover well the observational data.
Since Lau et al. (2016) also uses absorption lines of other ions besides H i to constrain , for completeness, in Figure B1 of Appendix B we give our models predictions for various ion columns as a function of impact parameter in comparison to their results.Specifically, we show predictions for column densities of Si ii (16.3 eV), Si iv (45.1 eV), C iv (64.5 eV), N v (97.9 eV) and O vi (138.1 eV) as a function of impact parameter.Given that most of the data are limits, we cannot draw definite conclusions from these plots, i.e. none of the models stand out as preferred.The most important points are (i) the agreement with observations of the predictions for the aforementioned favored model with low accretion rate  = 0.1 and  = 60 o (central column), and (ii) the exclusion of the model with the highest radiation input ( • = 10 9 M ⊙ ,  = 1.0) injected isotropically because of the tension with the Silicon columns (top right panels).
To conclude this section, we condense the behavior of our AGN models into one parameter: the covering fraction   of optically thick absorbers ( H i > 10 17.2 cm −2 ) within the projected CGM region, and compare it with QPQ estimates (Prochaska et al. 2013a,b).In a work using only the QPQ spectra with signal-to-noise ratios S/N>9.5 around the Ly restframe wavelength of the absorbers, Prochaska et al. (2013a) concluded that   ( < 200kpc)=0.64+0.06  −0.07 .This estimate is based on 50 sightlines (32 of which are optically thick), with the foreground quasars having redshifts between 1.63 and 2.97, and a mean redshift <  > QPQ5 = 2.06.Even assuming that quasars sit in halos of fixed mass at all epochs, as usually done, the limit radius of 200 physical kpc will cover different fractions of the projected virial radius at the redshifts of the sample.To correct for such aperture effects given the large  range of the QPQ sample, we take a fixed comoving aperture set as  lim = (1+ <  > QPQ5 )200kpc=612 ckpc and recompute the observational   based on the table 4 of Prochaska et al. (2013b).We get that 51 absorbers with S/N>9.5 have (1 + ) ≤  lim , and 28 of them are catalogued as optically thick.This translates into a new estimate of   (<  lim ) = 0.55 ± 0.07, with the uncertainty given by the binomial distribution.Figure 10 nicely shows how only small ionization cone opening angles are compatible with observations quantified by   (<  lim ) = 0.55 ± 0.07, and that the stronger the AGN's SED, the smaller the allowed  values 6 .The upper limits on  given in Table 1 are computed as the intersections between the upper 2 percentile curves of the four AGN SED models (dashed colored curves) and the lower 2 limit on the observationally estimated   (<  lim ) of 0.41.
It is important, though, not to over interpret these upper limits for various reasons.First of all, we use only one simulation, and because it's a zoom-in we can only consider a limited high-resolution region around the main galaxy.Second, most QPQ foreground quasars have smaller redshifts than the simulation, and it is expected that   increases with redshift (e.g.Rahmati et al. 2015).However, the low number statistics behind the observational estimate of   makes measuring a -trend unreliable.Third, the most generous ranges for dark matter halos hosting quasars span from 10 12 to 10 13 M ⊙ .This mass uncertainty coupled with the wide redshift range of observations means that the probed quasar halos can be in significantly different evolutionary stages, with vastly different amounts of cool gas in the CGM (e.g.Dekel & Birnboim 2006).Finally, the uncertainties in the quasars' redshifts impose lower limits on the wavelength (velocity) window used to model the Ly-line in absorption, which corresponds to physical depths a few times larger than the viral radii of the dark matter host halos.This latter limitation biases high the observationally derived   (<  lim ).All these caveats are discussed by Prochaska et al. (2013a) and Prochaska et al. (2013b), and quantified in some of the numerical works using large-box simulations, like the one of Rahmati et al. (2015).

CGM OBSERVABLES IN EMISSION
A few previous works reproduced the physical scales and level of emission in Ly seen in high- observations using a variety of simulations and modeling approaches, from dealing with radiation transfer effects only in post-processing (e.g.In this work, we compare our simulation with the observations of Fossati et al. (2021), who targeted 27 quasars at a median redshift  med = 3.75 with MUSE, and derived average stacked emission line profiles not only for Ly, but also for the rest frame UV lines He ii1640Å and C iv1549Å.He ii is a recombination line which requires a hard ionization spectrum to excite the gas ( ≥ 54.4eV = 4 Ryd), and therefore can be a good tracer of AGN radiation effects.Ly and C iv are resonant lines, whose modeling is much more complicated (e.g.Dĳkstra 2019).The C iv emission, also requiring high energies ( > 4 Ryd), encodes information on the metallicity of the gas and therefore its extent can provide an estimate of the scale out to which the halo is metal-enriched.
To obtain mock emission line maps we proceed as follows.We make the simplifying assumption that a gas particle of density   , temperature  and metallicity  exposed to a radiation field   has an emergent emissivity at a particular line frequency  0 : where   0 ,intrinsic is the intrinsic line emissivity in units of erg s −1 cm −3 computed with cloudy, and   0 the local optical depth at the line frequency.As our photoionization models have a hydrogen column stopping criteria, and we consider particles to be represented by the physical properties of the last cloudy zone, we compute the local optical depth at the line frequency as: where   0 and   0 are the column and volume number density of the   15, with the maximum likelihood parameters of Table 2.
ionic species producing the line (e.g.H i for Ly),   0 is the absorption cross-section at line center, and  zone is the corresponding depth of the last cloudy zone, which is also a function of ( H , , ,   ).
The Ly absorption cross-section at line center is a function of temperature only (e.g.Dĳkstra 2019): and     =  Ly  H i  zone , while  H i is computed as described in the previous section.
The other emission lines we are interested in are He ii1640Å, C iv1549Å, and H.To compute the absorption cross-sections for these other lines, we use the fact that the line center optical depth is proportional to the oscillator strength  7 .Thus, for H we have with  Ly = 0.416 and  H = 0.6958, the ratio of the statistical weights of the first two levels of hydrogen  2 / 1 = 2/8,   the 7 All the values of atomic physics used in this manuscript are taken from the NIST database https://www.nist.gov/pml/atomic-spectra-database Boltzmann constant, and 10.2 eV the difference in energy between the two levels.The relevant column density is the same as in the Ly case: For the He ii line we have to also consider the difference in atomic mass between helium and hydrogen, such that: where  2 / 1 = 8/2,  He ii = 0.6958 and  He = 4 and  H = 1 are the two atomic masses.In this case the difference in energy levels is 40.8 eV.The C iv line is like the Ly line, but for the triply-ionized carbon atom, such that: with the oscillator strength  C iv = 0.19 and  C = 12.For this line we have to make a further assumption, as the Gasoline2 version we used to run this simulation does not trace the carbon enrichment explicitly.Therefore, we assume that particles have the same carbon to oxygen number density ratio as the Sun's photosphere (Asplund et al. 2009), and use the oxygen abundance of each particle, which is traced by the simulation code, to compute  C iv taking into account the different atomic masses of the two elements.
Once we have the emergent emissivity for each line, we compute for each particle the luminosity at the particular line frequency  0 as: where  and  are the mass and mass volume density of the particle.Finally, the observed flux from a particle becomes: where   is the luminosity distance.For  = 3 and our assumed cosmology   =26070 Mpc.We stress that for the resonant transitions Ly and C iv, we are not performing a full radiative transfer calculation.Such a calculation would require an even higher resolution in the CGM and a realistic dust distribution to converge (e.g.Camps et al. 2021).Nevertheless, our approach takes into account this process in an approximated way.
Following the results presented in Costa et al. (2022) to reproduce extended emission around high-z quasars, we assume that Ly and C iv photons are able to escape to CGM scales through paths of least resistance (the ionization cones in this work), and are able to resonantly scatter in the CGM (inclusion of line emissions in the quasar input spectrum and continuum pumping in the Cloudy calculation).Finally, our calculation of the emergent emissivity (Equation 6) for the resonant lines is conservative as it assumes the cross-section at line center, which maximizes absorption.In reality, local absorption of resonant photons is quickly followed by the re-emission of another photon.These re-emitted photons could escape towards the observer in one or more scatterings.We stress that the local absorption correction of Equation 6 is important when the photoionization source is the UVB, as this spectrum results in larger local column densities of neutral hydrogen than the AGN ones.In fact, in the CGM region, the local absorption correction is minimal if the photoionization source is an AGN, as demonstrated by Figure C1 of Appendix C. Lastly, no spatial redistribution of resonant photons is taken into account.However such phenomenon should not affect strongly the observables used in this work (i.e.stacked images and the SB profiles) as its effects are (i) washed out by the stacking of several objects (resulting in a symmetric glow) and (ii) stronger on scales similar to the radial bins used in the following analysis.
The sample of nebulae observed by Fossati et al. (2021) comprises 27 quasar fields, and these authors stack all of them together in order to increase the signal-to-noise (S/N) of the Ly, He ii and C iv emission lines (see their Figure 6 with the median stacked surface brightness maps).To mimic their observations, and since we have only one simulation, we pick 27 random directions through the halo center to place the AGN ionization bi-cone, and for each of them we randomly select as LOS of our mock observations, a direction within this bi-cone.We then redshift the simulation to the median redshift of the observed quasar sample  med = 3.75 and chose for our mock images the same pixel scale as in observations (0.2 arcsec).
For each pixel in one of these individual mock observations we sum the line flux contribution (Equation 13) from all the particles with 2D projected positions within the pixel, and along a LOS length of 210 kpc (-105 kpc to +105 kpc).The obtained flux in each pixel is then divided by the pixel area to obtain surface-brightness (SB) maps.Each of the 27 maps for a particular emission line is then smoothed with a Gaussian kernel with FWHM = 0.65 arcsec (or ≈ 5 kpc) to take into account the average image quality of the data in Fossati et al. (2021).Finally the 27 SB maps have been median stacked to compare them to the observations.Figure 11 shows how the obtained median stacked SB maps change with the ionization cone opening angle for the AGN SED model with  • =10 9 M ⊙ and  = 0.1.The top panels give the case where only the UVB acts as photoionization source, while the bottom one the extreme case where the AGN illuminates all the CGM gas.The color bar in all panels is the same as the one in Figure 6 of Fossati et al. (2021).In all three lines, the emission is brighter for larger , while for a fixed  the Ly SB is larger than that of C iv, which in its turn is larger than that of He ii.Interestingly, at the sensitivity of Fossati et al.'s observations, a Ly glow would be observed even if the AGN has a minimal impact on the CGM gas ( = 0), but there would be no detectable emission in He ii and C iv.
To compare in a quantitative manner the simulation with the observations, Figure 12 shows the SB profiles as a function of projected radius  for the same AGN SED model with  • =10 9 M ⊙ and  = 0.1, as functions of  (shown only nine profiles with  corresponding to the values of the colorbar ticks).For the Ly (left panel), the simulation reproduces well the shape of the observed stacked profile (black points).The normalization of the simulated profiles scales with the  angle, with  ∼ 45 o falling on top of the observations.For the He ii and C iv (center and right panels), the simulated profiles change not only their normalization, but also their shape with varying .In the case of He ii, for which only the first two radial bins are above the noise level in the stacked profile of Fossati et al., it seems that also the model with  ∼ 45 o agrees best with observations.The observed stacked profile of C iv seems to require a larger  than the one of Ly, for a fixed AGN model.
We fit the Ly profiles with the following type of power-law: where SB 0 is the normalization,   the scale-length, and  the exponent.In a first step, we fitted the Ly profile of each AGN- model using Equation 14.Since the profiles for various  angles are approximately self-similar for a given AGN SED model, we fixed the scale-length and the exponent to their mean values over all s (see values in Table 2).Then, for each curve we computed the normalization as a simple average log(SB 0 )=⟨log(SB data )+(1+/  )⟩.These normalizations vary smoothly with the angle  following: as shown in Figure 13.The parameters of equations 14 and 15, for all four AGN SED models, are given in Table 2.In the literature, the Ly profiles are usually fit with exponentials or power-laws (e.g.Arrigoni Battaia et al. 2019;Farina et al. 2019;Fossati et al. 2021).
We have tried the already published functional forms, but Equation 14is performing better for our simulated profiles, and it is quite a good approximation for the observational profile of Fossati et al. (dashed black curve in Figure 12).In section 8.2 we will discuss in more detail the simulated He ii profiles, which are not as self-similar as the Ly ones.
A tension between observations and the simulated profiles is visible for the C iv emission, with the latter having steeper profiles than the former.While the CGM metal enrichment of the simulation is compatible with absorption studies ( CGM / ⊙ ≥ 0.1, e.g.Prochaska & Hennawi 2009;Lau et al. 2016, and bottom panel of Figure 3), the metallicity gradient is probably too steep in comparison with observations (right panel of Figure 12).This is a possible indication of a missing mechanism to push metals far away enough from the ISM.Thermal AGN feedback can produce such a strong enrichment of the CGM (e.g.Sanchez et al. 2019;Nelson et al. 2019).One should also keep in mind that the details of metal diffusion in numerical experiments are far from being settled (e.g.Wadsley et al. 2008).Also, as we discussed previously, this simulation does not trace explicitly the carbon enrichment.Running it with the new chemical enrichment implementation of Buck et al. (2021) will remove this uncertainty.Lastly, the steeper simulated profiles could be also due to the fact that we do not take into account in our modeling the spatial diffusion due to resonant scattering.However, giving that the models reproduce the Ly profile, which is expected to be even more affected by resonant scattering, indicates that this effect should be marginal for C iv.Because of these uncertainties, we do not provide in this work fits to the C iv profiles.
We can, however, quantitatively compare the simulated profiles with the observed ones, without any fitting, and constrain the ionization cone opening angle for each AGN model.To find the optimal angle α we compute the posterior probability (|data) ∝ (data|)(), where 'data' are, one at a time, the line emission stack profiles of Fossati et al. (2021) for Ly, He ii, and C iv, and we assume a flat prior () = 1/180 o .The likelihood of the observational data given the model is then: where N is the normal distribution8 , SB ,model the model prediction (simulation ⊗ AGN SED) at the impact parameter   of the observed surface brightness SB ,obs , and  ,model the corresponding dispersion of the model in logarithm.In Figure 14 we show the posterior probability distribution of  for the four AGN SED models (different colors), as constrained by the stacked emission profiles above the noise level for the Ly (left), He ii (center), and C iv (right).
Table 1 gives α with its confidence interval constrained for all combinations of SED and observed emission line profiles.Given the resonant scattering effects that can affect the Ly and C iv lines, in a best case scenario one should trust more the constraints provided by recombination lines like He ii.However, He ii has only been detected in some of the Ly nebulae, and typically on much smaller scales.For this reason and also because of the uncertainties on the CGM metal enrichment in our simulation (which affect the predicted C iv profiles) we decided to give separately the constraints for each observed emission line.In an ideal case, maybe feasible in the future, one would combine all the observational constraints into one single likelihood.All these caveats aside, one of the four SED models tested ( • =10 9 M ⊙ and  = 0.1, highlighted in grey in the table) leads to compatible α among all three emission lines, and, within their uncertainties, these values are less than the upper limit provided by the covering fraction of optically thick absorbers.Thus, for this best AGN model we can conclude that α ∼60 o is a value in good agreement with observational studies in emission, and it is also within the limits ].The first three lines give α in units of degrees estimated from fits to the stacked SB profiles of Fossati et al. (2021), while the last line gives the upper limit of  that agrees with the covering fraction of optically thick absorbers of Prochaska et al. (2013b).given by observations in absorption (see Section 6).E.g.Trainor & Steidel (2013) used Ly emitters in the vicinity of hyperluminous QSOs to constrain the activity period to 1 Myr≤  Q ≤20 Myr and the ionization cone opening angle to  ≥60 o .Using a combination of absorption and emission measurements of the QPQ sample, Hennawi & Prochaska (2013) find that  ≤90 o .The constraints on  of both these works overlap in a narrow range with each other as well as with the older study of Fritz et al. (2006).
It is important to mention that using only the AGN continuum, that is the SED without the contribution from the stacked line emission spectrum of Lusso et al. (2015), leads to lower SB levels in Ly and C iv.In its turn, this implies that larger ionization cone opening angles are required to match the observations, but these larger s are in tension with absorption studies.Figure C2 in Appendix C illustrates the differences between SB Ly for the AGN SED with and without the stacked emission lines.Thus, we can conclude that resonant pumping of the AGN spectrum is essential to explain the observed SB levels in the resonant lines, and their ratios with the HeII recombination line, similar to the conclusion reached by Fossati et al. (2021).
To conclude this section we also look at how representative are the four AGN SED models for observed quasars.Figure 15 gives the distribution of SDSS quasars from the catalog published by Rakshit et al. (2020) in the parameter space log( • )-log().To construct the distribution we use only the quasars with the best quality flag (0), and restrict the redshift to  ≥ 2.0, which is the range that covers the observational samples of Lau et al. (2016) and Fossati et al. (2021).These criteria restrict the observational quasar sample to 191073 objects.The four AGN SED models we tested in this work are marked with blue diamonds, and we can notice that our preferred model (log( • /M ⊙ )=9.0 and log()=-1.0) is close to the 2D peak of the SDSS distribution, meaning that it is also the most representative model for observed quasars at  ≥ 2.

Predictions on the covering fraction of Ly𝛼 emitting gas around quasars
In Section 7 we used the observed emission profiles to constrain the central quasar properties and opening angle, favoring an AGN model with  • = 10 9 M ⊙ ,  = 0.1 and α ≈ 60 degrees.Here, we compute the covering fraction of the Ly emitting gas and compare it to observational estimates.Arrigoni Battaia et al. (2019) computed the average covering factor profile  C () of each Ly nebula in their sample of 61  ∼ 3 quasars' fields in logarithmic radial bins, defining  C as the ratio between the area occupied by Ly emission above a SB threshold of ≈ 10 −18 erg s −1 cm −2 arcsec −2 (S/N= 2 in those observations) and the total area of each radial bin.Assuming equal weights for the profile of each individual target, they derived an average covering factor profile.We apply an identical procedure to the same 27 random orientations constructed for the simulated Ly SB profile in Section 7 assuming an AGN SED with  • = 10 9 M ⊙ ,  = 0.1, and repeat the calculation for three critical opening angles 0 o , 60 o , 180 o , corresponding to no quasar ionization, the favored opening angle, and isotropic AGN illumination, respectively.
Figure 16 shows the results of this exercise.The predicted covering factor profile for Ly emission above 10 −18 erg s −1 cm −2 arcsec −2 (blue) is much lower, in agreement, and much higher with respect to the observations (black) for an opening angle of 0 o (left panel), 60 o (central panel), and 180 o (right panel), respectively.The agreement between the observed profile and the simulation for an opening angle of 60 o further confirms the analysis in Section 7, i.e. the quasars illuminate on average a relatively small volume fraction of their halos (2Ω/4 =  2 (/4) ≈ 0.13), and dominate the powering of the Ly emission.
By construction, the covering factor profile depends on the depth of the observational data.It is therefore interesting to make predictions on how such a profile change for more sensitive data sets.We provide forecasts for the average covering factor of Ly emitting gas as a function of SB threshold.Specifically, Figure 16 shows how the covering factor profile varies by lowering the Ly SB threshold to 5 × 10 −19 (orange), 10 −19 (green), 5 × 10 −20 erg s −1 cm −2 arcsec −2 (red).Even ultra-deep observations, reaching SB levels of 5 × 10 −20 erg s −1 cm −2 arcsec −2 (e.g., Bacon et al. 2021), would result in at most an average  C ∼ 0.4 close to the expected virial radius of the DM halo (vertical dashed line).Importantly, these predictions take into account only the radiation from the central quasar.Additional sources would contribute to the extended Ly emission resulting in higher covering factors at larger radii (e.g., see Figure 7 in Arrigoni Battaia et al. 2019).

QSO properties from He ii extended emission
The importance of detecting He ii emission in extended quasars' Ly nebulae and extended Ly nebulae in general has been always emphasized since their discovery (e.g., Heckman et al. 1991b;Villar-Martín et al. 2003;Yang et al. 2006;Prescott et al. 2015;Arrigoni Battaia et al. 2015a;Herenz et al. 2020).He ii is a non-resonant line and it is conveniently located in the observed optical range for high-redshift sources.Because of this, He ii can be used to: assess how strong are the resonant scattering effects on the Ly line, constrain the ionization state of the gas and hardness of the impinging ionization spectrum, and study the gas kinematics.In this work, we have shown that He ii emission is the cleaner tracer of the AGN ionization effects on the surrounding CGM, and that together with other emission and absorbtion lines, it can provide estimates for the central AGN's parameters (SMBH mass, Eddington ratio, and ionization cone opening angle).This first estimate showcases the potential in the use of CGM observables to measure the properties of the central engine and could provide independent and complementary constraints on AGN properties (with respect to the usual methods using quasars emission lines, e.g., Farina et al. 2022), especially at the earliest epochs where these values are most needed to understand the rapid formation of SMBHs (Fan et al. 2022 and references therein).As the He ii emission is very faint, and hard to detect in individual objects (e.g.Guo et al. 2020;Fossati et al. 2021), we provide fits to its SB profile as predicted by our modeling, and quantify how much better we could do with current instruments.
Figure 17 shows the simulated He ii profile for the favored AGN SED model, this time together with the maximum likelihood fits using Equation 14. Contrary to the Ly profiles, the He ii requires varying all three parameters (normalization, scale-length and exponent).Table 3 gives all these parameters for all four AGN SED models, and the root-mean-squared deviation  of the simulated data points from the fitted function, to quantify the goodness of fit.
The center and right panels of Figure 17 show the profiles of SB ratios at fixed accretion rate (=0.1,center) to gauge the effect of changing the black hole mass, and the profiles of SB ratios at fixed black hole mass ( • =10 9 M ⊙ , right) to gauge the effect of changing the accretion rate.We can notice that changing  at fixed  • has only a mild effect on the SB He ii with an increase of 20-25% when increasing  from 0.1 to 1.0 for the larger ionization cone opening angles ( ≥ 90 o right panel of´Figure 17).Increasing  • at fixed  by an order of magnitude has a larger effect on SB He ii (center panel).The SB ratios for reasonable ionization cone opening angles ( ∼ 60 o ) are very sensitive to the projected radius at which they are measured.All this, together with the comparison between our simulated profiles and the stacked data of Fossati et al. (2021), proves that we need more sensitive observations in order to use log(SB 0 ) -19. 29 -16.99 -15.94 -15.77 -16.29 -16.30 -16.41 -16.42 -16.42  29 -16.52 -15.71 -15.70 -16.22 -16.17 -16.28 -16.29 -16.28  29 -16.27 -15.58 -15.68 -16.20 -16.15 -16.25 -16.27 -16.25  In each panel, the horizontal dashed grey line gives the ratio value for Case B recombination SB H /SB Ly =0.12 for =10 4 K and typical CGM metallicities of 0.5 ⊙ .At fixed angle  we collapsed together all the pixels of 27 images, representing 27 random orientations of the ionizing cone, and implicitly the LOS.Bottom: average H SB profile (left), and median SB H /SB Ly profiles (right) for the same AGN SED model and different  angles.The shaded areas in the right panel enclose the regions between the 16 th and 84 th percentiles of ratios, while the horizontal dashed grey line gives the same Case B recombination ratio as in the top panels.
CGM data for constraining firmly the AGN properties.For example, to push the observations to a 2 (3) significance at the SB level of 10 −20 erg s −1 cm −2 arcsec −2 for the He ii transition, one would need to target additional 315 (756) objects to the sample of Fossati et al. (2021) with ∼ 4 hours MUSE observations, following their quoted sensitivity and assuming their redshift distribution9 .Such a large sample of sources should also be able to unveil the targets with the brightest He ii nebulae (e.g., Guo et al. 2020) which could then be used to study the kinematics of the quasar's CGM (e.g., Zhang et al. 2023) and obtain estimates for individual quasar properties (e.g., ionization cone opening angle).

Constraining the Ly𝛼 emission processes with H𝛼 glows
Observing the H6563Å transition from the quasar's CGM can provide similar information as He ii.In particular, H emission should be mainly produced through recombination and should not be subject to resonant scattering effects under most astrophysical conditions.For this reason, it can be used to pin down the powering mechanism of Ly emission, constrain the gas physical properties (e.g., density) and the gas kinematics.Unlike He ii, however, the detection of H emission is not possible from the ground for  ≳ 2.5 due to the Earth atmosphere (e.g., Hayes 2019), and space telescopes such as the James Webb Space Telescope (JWST, Gardner et al. 2006) are needed to overcome this obstacle.To our knowledge, there are yet no published constraints on H emission for the CGM (on scales > 20 kpc as in this study) of  ∼ 3 quasars.Therefore, we provide predictions for the average H SB levels expected in the model favored by our analysis, i.e. a CGM illuminated by an AGN SED with  • = 10 9 M ⊙ and  = 0.1.
We construct H SB maps as done for the other transitions and assuming three different opening angles  = 0 o , 60 o , and 180 o .As for the other emission lines, we pick 27 random directions of the ionization cones and for each of those a random line of sight within the cone.Importantly, these directions are the same as those used for the construction of the Ly SB maps, so that we can now compare the emission levels of the two transitions.First, we focus on the distribution of SB H /SB Ly ratios as a function of SB Ly on CGM scales (radii  > 20 kpc).The top panels in Figure 18 show the distribution of the ratio computed for CGM regions (pixels with area of 0.2 × 0.2 arcsec 2 ) for the three s, 0 o (left), 60 o (center), and 180 o (right).For reference, we indicate the Case B recombination value SB H /SB Ly = 0.12 for  = 10 4 K and a typical CGM metallicity 0.5 Z ⊙ computed with cloudy10 .This approximated regime occurs when the gas is optically thick to Lyman-line photons, which are expected to undergo several scatterings and be converted into Balmer line photons and Ly or two-photon emission (e.g., Ferland 1999).
We find that for all angles  the predicted SB H /SB Ly ratios do not follow uniquely the frequently invoked Case B value.For  = 0 o , i.e. the UVB-only case, most of the CGM follows roughly case B, but there are regions where the ratio is lower and much higher.The regions with a lower ratio are in between case B and case A (according to cloudy at  = 10 4 K case A gives log(SB H /SB Ly ) = −1.8),while the regions with higher ratios are patches where the local absorption correction is the strongest (Figure C1).Introducing the ionization cones of the quasar as described in Section 5 widens the distribution of ratios towards lower values because continuum pumping increases the Ly SB.For the extreme case where the AGN shines on all the CGM most of the regions have values SB H /SB Ly ≲0.01 (top right panel of Figure 18).Removing the quasar's emission lines from the input spectrum (Section 5) would result in higher ratios as the continuum pumping of the resonant lines would be less significant (Figure C3 in Appendix C).
The bottom panels of Figure 18 give the predicted average H SB (left) and the median SB H /SB Ly (right) as functions of radius for the three  angles.The shaded areas in the bottom right panel mark the 16 th to 84 th ratio percentiles.Overall, considering the full distribution of allowed ratios for  = 60 o and the fact that small opening angles could be in place (Table 1), our models suggest that the easiest H detections would be those corresponding to SB H /SB Ly ∼ 0.1, close to a case B recombination scenario.However, statistical samples of observations should result in lower ratios.
With the advent of the JWST Near Infrared Spectrograph (NIR-Spec) observations can test our models and hence obtain firmer constraints on the physical properties of the  ∼ 3 quasars' CGM and on the physical properties of the quasars themselves (i.e.,  • , , and opening angle).While observations of individual systems are necessary as pilot studies, large samples of sources are needed to really pin down the physics of the cool CGM.First results on this front are expected for quasars at the highest redshifts,  ∼ 6, with known Ly nebulae (e.g., Farina et al. 2019) which have been targeted with JWST observations during cycle 1, and will also be targeted during cycle 2.
At lower redshifts, for which both Ly and H are accessible from the ground, a recent work targeted a  ∼ 2.26 quasar and reported a detection for both transitions (Langen et al. 2023).These authors used the IFU KCWI and a longslit spectrum with the MOSFIRE instrument on the Keck telescope to target Ly and H, respectively.A MOSFIRE 1 ′′ slit has been placed on the brightest portion of the Ly nebula, which is also the brightest in their sample.The H detection extends for 2.5 ′′ (or ∼ 20 kpc) on one side of the quasar.After extracting a pseudoslit from the KCWI data corresponding to the MOSFIRE one, they report a ratio of SB H /SB Ly = 0.27, in a region of 1 ′′ × 5.5 ′′ .If this value would be present throughout the whole Ly nebula, it would be in tension with our predictions.However, we find surprising their quoted value of the Ly flux  Ly = 3.5 × 10 −17 erg s −1 cm −2 within the constructed pseudoslit, because it corresponds to only 7% of the total Ly flux of the nebula.This flux clearly does not agree with the SB map presented in their Figure 2: only the 1 ′′ × 1 ′′ region at the high SB peak should already have at least a factor of two higher flux than what they report for the slit.A rough estimate of the flux within the drawn pseudoslit on the same image suggests instead that the Ly flux in that 1 ′′ × 5.5 ′′ region should be ∼ 20% of the total.This would make the ratio SB H /SB Ly ∼ 0.09, a value more consistent with our modeling.We stress that  ∼ 2 quasar's Ly nebulae are known to be on average dimmer than those at  ∼ 3 (e.g., Cai et al. 2019), so the conditions in the CGM of quasars at different redshifts could be different (Arrigoni Battaia et al. 2019).
Once again, the discussion of this example showcases how critical it is to acquire a large statistical sample of observations targeting the H transition from the high-redshift quasar's CGM.Importantly, we need observations of the entirety of the Ly nebulae to constrain radiative transfer effects well.Data achieving SB levels of ∼ 10 −19 erg s −1 cm −2 arcsec −2 should be able to detect on average H emission out to ∼ 30 kpc (Figure 18), if the densities in the simulation are a good representation of reality (Section 8.4).

Density distributions of the cool CGM emitting Ly𝛼
The SB levels of quasars' Ly nebulae require the emitting gas to have high densities (≳ 1 cm −3 ; e.g., Heckman et al. 1991a;Cantalupo et al. 2014;Arrigoni Battaia et al. 2015b;Hennawi et al. 2015), or very broad density distributions (e.g., Langen et al. 2023 and references therein) for a pure recombination scenario.Testing for the presence of such CGM densities can be done in two ways: via sub-grid prescriptions for the unresolved dense cool CGM gas in cosmological simulations (e.g., Hummels et al. 2019), or via very high-resolution simulations that can trace the survival of such dense cool structures in a hot medium (e.g., McCourt et al. 2018;Gronke et al. 2022).On the other hand, if quasar's photons contribute to the emission in resonant lines (Ly, C iv), the required densities would be lowered.Using a suite of high-resolution (maximum spatial resolution ∼80 pc), radiation-hydrodynamic cosmological simulations targeting  ∼ 6 quasar host halos and post-processed with a Ly radiative transfer code, Costa et al. (2022) was able to match the observed Ly SB profile.Importantly, this work showed that quasar's photons can efficiently scatter in the quasar's CGM indicating that extreme densities may not be necessary to match observations.Motivated by this work, in our analysis we used a quasar spectrum including the emission lines to allow continuum pumping to be effective (Section 5), and we show that the predicted observables (column densities and SB profiles) match current observations.Therefore, it is interesting to look at what kind of densities could be recovered from the mock observations in comparison to the ground truth of the simulation.
We constructed the hydrogen number densities radial distributions in Figure 19 by creating 27 random orientation maps (as done for all the results in Section 7) with the line-of-sight mass (black curve) or Ly luminosity (colored curves) weighted  H .For the mass-weighted  H we used only the cool gas particles ( <10 5 K).The estimates for  H,Ly should mimic hydrogen number densities that could be extracted from observations.While most of the CGM emitting Ly is characterized by densities ∼ 10 −2 cm −3 , the distributions in Figure 19 show tails at high densities, up to ∼ 10 cm −2 for the models with an AGN.These high densities are responsible for the brightest Ly emission in the inner CGM.On the contrary, in the model including only the UVB as photoionization source there is almost no contribution for  H > 0.1 cm −3 as those regions are mostly neutral with such a faint impinging continuum.The  H,Ly distributions for the two cases where the AGN is on follow better the intrinsic cool gas mass-weighted density distribution than the UVB case.At large radii the intrinsic densities are lower than the luminosity-weighted ones.
Re-running the simulation at higher resolution will allow resolving higher densities in the CGM, and likely increase the H i column densities by a few tens of percent (e.g.Hummels et al. 2019;van de Voort et al. 2019;Peeples et al. 2019;Suresh et al. 2019;Bennett & Sĳacki 2020).As a consequence, the ionization cone opening angles permitted by the covering fraction of optically thick absorbers would increase.Simultaneously, higher densities would affect the emission levels, requiring smaller angles to match observations and possibly a smaller contribution from continuum pumping.Nevertheless, we should also keep in mind that running the same simulation, but with this new AGN EM feedback model on-the-fly might change the phase space distribution of CGM gas, and therefore, also the distribution of densities.Finally, we should test our predictions against radiation transfer post-processing, which might naturally provide the means to bring in complete agreement the  angles required by the Ly and He ii emission lines.

SUMMARY
We used a high resolution (gas particle mass resolution of 3.6×10 4 M ⊙ , and gravitational softening down to 72 pc) zoom-in cosmological simulation of a massive (10 12.5 M ⊙ ) halo at  = 3 to study in post-processing the direct effects of the radiation field generated by a central accreting black hole on its CGM.The simulation has been run with the galaxy formation model used for the NIHAO project (Wang et al. 2015), without any AGN feedback.For the AGN SED we used the physically consistent continuum model of Done et al. (2012), with all but two parameters fixed to the median values over a sample of type I AGN (Jin et al. 2012), to which we added the stacked quasars line emission spectrum of Lusso et al. (2015).The two free parameters of the AGN SED are the SMBH mass and the Eddington accretion ratio, and for both we explored two different values (10 8 and 10 9 M ⊙ masses, and 0.1 and 1.0 accretion ratios), which together bracket the estimates for a large sample of SDSS quasars at  > 2 (Rakshit et al. 2020).To model the gas response to the radiation field we run large grids of cloudy models with the incident radiation fields given by our choice of AGN SEDs and distances covering all the CGM of the halo, while the control models were run with the  = 3 UVB of Khaire & Srianand (2019).
The radiation from AGN is not emitted isotropically, and to study this anisotropy we positioned bi-cones of given opening angle  at various randomly placed sightlines through the centre of the halo, and stacked any intrinsic or observational CGM property of interest over these sightlines.We consider only the particles within the bi-cones to be illuminated by the AGN, while the ones outside are illuminated by the UVB.Under the radiation of an AGN, a significant fraction of the cool ( < 10 5 K) CGM will not cool catastrophically, but instead heat up.This fraction steadily increases with the angle  reaching almost 100% for the maximum  = 180 o -when all gas is irradiated by the AGN -and for all but the weakest AGN SED studied ( • = 10 8 M ⊙ and  = 0.1), see Figure 6.For angles  = 50 o -60 o , half of the cool CGM is heating, while the other half is cooling catastrophically.From this exercise, it becomes clear that AGN radiation can stop or at least stall the cooling of CGM gas in massive halos (e.g.Ciotti & Ostriker 1997, 2007), and hence provide a viable negative (or preventive) feedback for star formation.
We then turn to observables.In particular, at high redshift, there are CGM observations for quasar host halos both in absorption (e.g. the QPQ survey Hennawi et al. 2006;Hennawi & Prochaska 2013;Prochaska & Hennawi 2009;Prochaska et al. 2013aProchaska et al. ,b, 2014;;Lau et al. 2016Lau et al. , 2018)), and in rest-frame UV emission lines (e.g.Arrigoni Battaia et al. 2019;Fossati et al. 2021).To compare with these observations, we constructed mock column density and ionization parameter profiles, and surface brightness profiles from SB maps in three emission lines (Ly, He ii, C iv).Similar to previous studies using simulations run without AGN feedback (e.g.Faucher-Giguère et al. 2016), we found H i column densities as high as those inferred from the QPQ sample (Lau et al. 2016) if the CGM is only photoionized by the UVB.A similar agreement can be found, also if we do consider a central AGN source, but restrict the ionization cone opening angle  to low values (see Table 1).Lau et al. (2016) also estimated ionization parameters for their absorbers, so we can use the parameter space  vs  H i to check which combination of photoionization sources better explains the observations.In this case, we find that only AGN models with low  are capable of simultaneously explaining the spread in 'observed'  and  H i values (Figure 8).The UVB only model instead leads to a very marked anti-correlation between  and  H i , with a scatter too small to encompass all observations.
Firmer constraints on the AGN model come from comparing the simulation with emission nebulae around quasars.We created mock SB maps for the rest-frame UV lines Ly, He ii, C iv, and from them stacked average radial profiles that can be compared to the stack of 27 quasar fields observed by Fossati et al. (2021).The simulated average profile of Ly has the same functional form as the one in observations (Figure 12 and Equation 14) with scale-length   = 26.8kpcand power  = 5.5, while its normalization depends on the AGN model and angle  (Table 2).This allowed us to use the observed profile for constructing the posterior probability of  for each of the four AGN models.Doing similar comparisons observation-simulation for the other two emission lines we can thus see which AGN model results in compatible  angle for all three emission lines.In our case, the model with  • =10 9 M ⊙ and  = 0.1 has the best compatibility between the three  values (see Table 1), with  = 60 o reproducing well the Ly and He ii observations, and falling within 1 of the angle required by C iv.The tension with the  required by C iv arises from the fact that the average stacked simulated C iv profile is steeper than the observational one, likely indicating that our simulation does not push enough metals beyond the inner CGM region.Rerunning the simulation with this EM AGN feedback on-the-fly will clarify if the CGM can be enriched to larger radii, or if a more explosive type of feedback is needed.
The comparison between the simulation and the CGM observations revealed a few very interesting facts: • It demonstrates how CGM observables can be used to pinpoint properties of the central AGN engine.In this respect, we think that deeper observations targeting non-resonant lines (i.e.He ii, H) will have the greatest power in constraining AGN properties (e.g.mass, accretion rate, opening angle; see Section 8.2).
• The very good match between the brightness (Figure 12) and morphology (Figure 16) of simulated and observed Ly nebulae indicates that the resolution of this simulation is likely enough to numerically model these types of CGM emission (see also Section 8.4 and Costa et al. 2022).The absorption observations, instead, seem to require only slightly higher resolving power for the cool CGM.
• Continuum pumping of the quasar's emission lines is needed to bring in agreement the predictions from emission and absorption observations (Figure C2).
• Detection of C iv nebulae around individual quasars at high significance could trace the metal enrichment to larger distances from the galaxy, and consequently help understand if a realistic CGM enrichment around single quasars can be achieved with only a gentle EM AGN feedback, or if explosive feedback is needed.
• Among the four AGN models, the one resulting in the best agreement with all the observations we considered ( • =10 9 M ⊙ and  = 0.1) is also the most representative of the wider high- population of quasars at the targeted redshifts (Figure 15).
• This favored model implies a balance between the CGM reservoir undergoing heating and cooling (Figure 6), hinting to the importance of self-regulation around AGNs also for the EM feedback, as it has been shown for other AGN feedback mechanisms (e.g., Gaspari et al. 2011).
To validate the exercise we made in this paper, we plan to: implement this EM AGN feedback in the gasoline2 code, run more similarly massive and slightly higher resolution simulations, both without and with this new implementation of feedback, and apply in post-processing a radiation transfer code that can follow also resonant scattering effects.Building up statistically significant samples of high resolution massive galaxies, run with different AGN implementations, is needed in order to fully exploit the potential of CGM observations for constraining not only feedback models, but also the properties of individual observed quasars (e.g.AGN mass, accretion rate, ionization cone opening angle).Future works have to understand whether also the properties of the simulated galaxies match observations, when the simulations are run with a milder or absent AGN ejective feedback together with EM feedback. .Column densities for five ionic species of increasing ionization energy as a function of projected distance from the quasar for three AGN spectra with increasing number of ionizing photons from left to right.From top to bottom we show the column densities of Si ii (16.3 eV), Si iv (45.1 eV), C iv (64.5 eV), N v (97.9 eV) and O vi (138.1 eV).The black points are the observational estimates of Lau et al. (2016) foreground quasars.Each percentile curve was obtained from the stacking of 10 different directions to the simulated halo (each within a randomly placed ionization cone).The line styles and colors represent the same percentiles and s as in Figure 8.

APPENDIX C: EFFECTS OF LOCAL ABSORPTION CORRECTION AND AGN EMISSION LINES ON CGM OBSERVABLES
The local absorption approximation we use has a bigger effect on the gas in the case where the UVB is the only source of ionizing photons.The top panel of Figure C1 shows the median neutral hydrogen local column as a function of distance from the center of the halo, in two extreme cases: all gas particles are photoionized only by the UVB (grey), and all are photoionized only by a 10 9 M ⊙ SMBH accreting at 10% of its Eddington rate (green).The bottom panel shows the corresponding distributions of local absorption corrections as functions of distance from the center.Shaded areas in both panels give the 16 th and 84 th percentiles.For the CGM region ( > 20kpc), the local absorption is negligible for the AGN photoionization case for all gas particles, except the dense ones belonging to the infalling satellite at  ∼  200 .For the UVB photoionization case, the local absorption is very important, especially in the inner CGM region ( ≤ 60kpc).
We tested how the UV line emission maps and derived surface brightness profiles change when we remove the line spectrum contribution (Lusso et al. 2015) from the AGN SED models.As expected, the He ii emission is not impacted at all, while the SBs of Ly and C iv decrease without changing the overall shape of the profiles.To quantify the differences in SBs with and without the line spectrum contribution, we give in Figure C2 the posterior probabilities of  for the three lines, for our preferred AGN model.If the input continuum does not include the line emission spectrum, higher ionization cone opening angles are needed to explain the SB profiles of Ly and C iv. E.g. the most likely  angle needed to explain the Ly SB is ∼ twice as large as our default case.However, such large ionization cone opening angles are excluded by absorption studies as we discuss in Section 6.

Figure 1 .
Figure 1.Dark matter and baryonic component maps for the simulated galaxy g3.16e12 at  ∼ 3.In clockwise direction from top-left we show: dark matter surface mass density, stellar surface brightness in rest-frame V-band, metal surface mass density, and hydrogen column density.The dashed white circles show the position of the 134 kpc virial radius.

Figure 2 .
Figure 2. Density-temperature phase space of simulation g3.16e12 within a 4× 200 box (top), and within the CGM of the main halo, 0.15 200 <  <  200 or 20 kpc <  < 134 kpc (center), and the ISM of the main halo  < 0.15 200 (bottom).The color bar in the left panels quantify the mass in hydrogen, while the one for the right panels the median metallicity in 2D density-temperature bins.

Figure 3 .
Figure3.Top: Density profiles of dark matter (black), cool ( <10 5 K, blue) and hot ( ≥10 5 K, orange) gas, and all gas (purple).Center: number density of hydrogen as a function of radius for the cool and hot gas phases.Bottom: metallicity as a function of radius for the cool and hot gas phases.In the middle and bottom panels, the points and their errors give the mean and standard deviation, the solid lines give the medians, and the shaded regions enclose from the 16 th to the 84 th percentiles (darker shade), and from the 2 nd to the 98 th percentiles (lighter shade).We also show the collapsed density distribution functions of hydrogen number density and metallicity in the right panels.The vertical grey line marks the virial radius  200 .In all panels, we exclude the ISM region  < 20 kpc, where  is the 3D radius.
cor .In summary, based on the BLS1 sample of Jin et al. (2012), our assumed AGN SEDs have seven of nine parameters fixed to the following values:  = 0  out = 10 4  g  cor = 26 g Γ

Figure 5 .
Figure 5. Net cooling rate as a function of temperature for a gas slab with density  H =10 −2 cm −3 and metallicity 0.1  ⊙ , exposed to the  = 3 UVB (black), or to the spectrum of a 10 9 M ⊙ SMBH accreting at 10 per cent the Eddington rate.Colored curves represent various distances between the radiation source and the gas slab, given in kpc within brackets in the legend.The legend on the right gives also the logarithm of the ionization parameter .The dashed/solid parts of the curves represent the conditions where gas is heating/cooling.

Figure 6 .
Figure6.The CGM gas mass budget as a function of the quasar ionization cone's opening angle .The various transparency levels of the curves represent the spectral model for the AGN.The three phases of the CGM represent different thermal states: gas that is heating ( cool < 0, light-red), gas that is cooling (0 <  cool <  ff , black), and gas in (quasi)hydrostatic equilibrium ( cool >  ff , blue).

Figure 7 .
Figure 7. Joint density-temperature phase space distributions for the three thermal states of CGM gas, assuming different quasar ionization cone's opening angles:  = 0 o representing the UVB only case (left),  = 60 o (center), and  = 180 o representing the case where all the CGM gas is impacted by the AGN radiation.The spectral model in all three panels corresponds to an AGN of  • = 10 9 M ⊙ and a  = 0.1 Eddington ratio.The contours of the kernel density estimate enclose equally spaced fixed fractions of particles in each of the three thermal states.

Figure 8 .
Figure 8. H i column density (top) and mean ionization parameter  (center) as functions of projected distance from the quasar  for three AGN model spectra with increasing number of ionizing photons (left to right; see on top of each column for the values of SMBH mass and Eddington ratio).The solid curves are the 50 th percentiles, the shaded areas mark values between the 16 th and 84 th percentiles, while the lower and upper dashed curves give the 2 nd and 98 th percentiles, respectively.Colors encode the ionization cone opening angle , as explained in the legends.Bottom: logarithm of the number of sightlines with given H i column density and mean ionization parameter (colored map) as a function of AGN model (left to right) and ionization cone opening angle (top to bottom).The black points in all panels are the observational estimates of Lau et al. (2016).

Figure 9 .
Figure9.Total hydrogen column density as a function of projected distance from the center of the halo.The median curve and the 1 and 2 shaded areas were obtained by stacking  =10 random LOS.The observational data points in black were computed as  H ≡  H i /X H i , where X H i is the ionization fraction of hydrogen, estimated via photoionization modeling(Lau et al. 2016).

Figure 11 .
Figure 11.Median stacked surface brightness maps for Ly (left), He ii1640Å(center), and C iv1549Å(right) for the AGN model  • =10 9 M ⊙ &  = 0.1 and three different ionization cone opening angles (top row gives the UVB only case, bottom row the AGN only, and the middle row an intermediate ).The black dots mask the central regions (radius of 10 kpc) where the observational stacks of Fossati et al. (2021) are dominated by the point spread function of the quasar.The images span 210×210 kpc 2 and the emission is computed as coming from a halo placed at the median redshift of the observations  med = 3.75.The images are smoothed with a Gaussian kernel of FWHM= 0.65 arcsec, corresponding to the average FWHM in Fossati et al. (2021).

Figure 12 .
Figure 12.Averaged radial surface brightness profiles for Ly (left), He ii1640Å(center), and C iv1549Å(right) for the AGN model  • =10 9 M ⊙ &  = 0.1 for different values of opening angle  corresponding to the ticks of the colorbar.The models are compared with the observations of Fossati et al. (2021).The emission is computed as coming from a halo placed at the median redshift of the observations  med = 3.75.The observational data points marked with an empty black circle are at the noise level or below.The curves in the left panel have the functional form SB fit =SB 0 (1+/  ) − , with   =26.8 kpc and =5.5, while SB 0 is computed averaging the normalization over all data points along a curve: log(SB 0 )=⟨log(SB data )+(1+/  )⟩.The normalization of the black dashed curve was computed from the observed data points, and has exactly the same value as the one of the =45 o model: SB 0 =10 −15.86 erg s −1 cm −2 arcsec −2 .

Figure 13 .
Figure 13.The normalization variations of the Ly profiles (colored points) with the angle  for the four AGN SED models (different colors).The curves follow Equation15, with the maximum likelihood parameters of Table2.

Figure 14 .
Figure 14.Posterior probability distribution of  for the four AGN SED models (different colors), as constrained by the stacked emission profiles above the noise level for the Ly (left), He ii (center), and C iv (right)

Figure 15 .
Figure15.Distribution of high- ( ≥ 2.0) SDSS quasars(Rakshit et al. 2020) in the parameter space of black hole mass and accretion rate.The black contours are kernel density estimates for the 20, 40, 68 and 95 percentiles, while the four black diamonds are the positions of the AGN SED models we tested in this work.The model favored by our analysis is highlighted in cyan.

Figure 16 .
Figure 16.Covering fraction of Ly emission as function of projected radius for various SB thresholds (indicated in the legend in units of erg s −1 cm −2 arcsec −2 ) for our preferred AGN SED model ( • =10 9 M ⊙ and  = 0.1) and three ionization cone opening angles:  = 0 o (left),  = 60 o (center), and  = 180 o (right).Colored curves in all three panels are the predictions from the simulation, while the black one is the average covering fraction profile for the QSO MUSEUM sample with a threshold of 10 −18 erg s −1 cm −2 arcsec −2 (ArrigoniBattaia et al. 2019).The vertical grey dashed line in all panels marks the virial radius  200 .

Figure 18 .
Figure18.Top: 2D distributions of 0.2×0.2arcsec 2 CGM (20<  <134 kpc) pixels in the space of Ly surface brightness and SB ratios of H to Ly for three different ionization cone opening angles (left to right) for our preferred AGN SED model ( • =10 9 M ⊙ &  = 0.1).In each panel, the horizontal dashed grey line gives the ratio value for Case B recombination SB H /SB Ly =0.12 for =10 4 K and typical CGM metallicities of 0.5 ⊙ .At fixed angle  we collapsed together all the pixels of 27 images, representing 27 random orientations of the ionizing cone, and implicitly the LOS.Bottom: average H SB profile (left), and median SB H /SB Ly profiles (right) for the same AGN SED model and different  angles.The shaded areas in the right panel enclose the regions between the 16 th and 84 th percentiles of ratios, while the horizontal dashed grey line gives the same Case B recombination ratio as in the top panels.

Figure 19 .
Figure 19.Median (solid curves), and 16 th (lower dashed curves) and 84 th (upper dashed curves) volume number density of hydrogen weighted with the Ly luminosity of each particle (colored curves), or with the hydrogen mass of cool ( <10 5 K) particles (black curves), as functions of projected distance from the center (left panel).The smaller right panel shows the corresponding  H distributions over the full projected radial range of the left panel.The AGN SED considered is the one for our preferred model, and the three colors represent three different ionization cone opening angles.

Figure C1 .
Figure C1.Local neutral hydrogen column density Δ  as a function of distance to the center of the halo (top), and the corresponding absorption correction for the Ly emission (bottom).Grey and green represent two different sources of photoionization as indicated in the legend and explained in the text.

Figure C2 .Figure C3 .
Figure C2.Comparison between the posterior probability distributions of  as constrained by Fossati's stacked emission line profiles above the noise level for the Ly (solid), He ii (dashed) and C iv (dotted), for our preferred AGN SED model ( • =10 9 M ⊙ &  = 0.1) with (turquoise) and without (black) the contribution of Lusso's stacked emission line AGN spectrum.
Hummels et al. 2019;van de Voort et al. 2019)h column densities >10 17.2 cm −2 for sightlines with projected impact parameters  ≤153 kpc, corresponding to a physical radial limit of 200 kpc at the mean redshift of the observational sample inProchaska et al. (2013a), as function of the ionization cone opening angle.Color curves represent the four AGN SED models, shaded areas mark the regions between the 16 th and 84 th percentiles, and dashed colored curves are the 98 th percentiles.The horizontal grey line gives our new estimate of   from table 4 ofProchaska et al. (2013b), using a fixed comoving aperture  lim =153(1+3)=612 ckpc, with the two shaded areas marking the 1 and 2 uncertainties.systems.Thus, if we consider that our simulation reaches high enough resolution to resolve decently the densities in the CGM (see the effects that resolution can have on  H i in the CGM, e.g.:Hummels et al. 2019;van de Voort et al. 2019), we can interpret the highest values

Table 2 .
Best fit parameters of Equations 14 and 15 applied to the simulated Ly profiles as function of AGN SED model and ionization cone opening angle.The units of SB 1 are erg s −1 cm −2 arcsec −2 .The last three parameters were obtained with the emcee package (Foreman-Mackey et al. 2013), but only   has an asymmetric marginalized posterior distribution.