# Large-scale geometry and topology of gas fields: Effects of AGN and stellar feedback

Carlo Schimd, Katarina Kraljic, Romeel Davé, Christophe Pichon

► To cite this version:

Carlo Schimd, Katarina Kraljic, Romeel Davé, Christophe Pichon. Large-scale geometry and topology of gas fields: Effects of AGN and stellar feedback. Astronomy & Astrophysics - A&A, 2024, 689, pp.A311. 10.1051/0004-6361/202449276. hal-04704639

HAL Id: hal-04704639

https://hal.science/hal-04704639v1

Submitted on 21 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Large-scale geometry and topology of gas fields: Effects of AGN and stellar feedback

Carlo Schimd<fnref n="1" />,*, Katarina Kraljic<fnref n="2" />, Romeel Davé<fnref n="3" />,<fnref n="4" />,<fnref n="5" />, and Christophe Pichon<fnref n="6" />,<fnref n="7" />,<fnref n="8" />

Received 19 January 2024 / Accepted 4 June 2024

## ABSTRACT

Feedback from stars and active galactic nuclei (AGNs) primarily affects the formation and evolution of galaxies and the circumgalactic medium, leaving some kind of imprint on larger scales. Based on the SIMBA hydrodynamical simulation suite and using the full set of Minkowski functionals (MFs), this study systematically analyses the time evolution of the global geometry and topology of the gas temperature, pressure, density (total, H, and H₂), and the metallicity fields between redshifts z = 5 and z = 0. The MFs show that small-scale astrophysical processes are persistent and manifest on larger, up to tens of Mpc scales, highlighting the specific morphological signatures of the relevant feedback mechanisms on these scales in the last ~12 Gyr. In qualitative terms, we were able to establish a ranking that varies according to the field considered: stellar feedback mostly determines the morphology of the pressure and density fields and AGN jets are the primary origin of the morphology of the temperature and metallicity fields, while X-ray heating and AGN winds play the second most important role in shaping the geometry and topology of all the gaseous fields, except metallicity. Hence, the cosmic evolution of the geometry and topology of fields characterising the thermodynamical and chemical properties of the cosmic web offers complementary, larger scale constraints to galaxy formation models.

**Key words.** galaxies: evolution – galaxies: formation – intergalactic medium – large-scale structure of Universe

## 1. Introduction

Astrophysical thermal and non-thermal processes due to stellar and active galactic nucleus (AGN) activity, jointly with hydrodynamical and tidal gravitational interactions acting on top of gravitational collapse, determine the physical evolution of galaxies and their local environments (Hummels et al. 2013; Tumlinson et al. 2017; Appleby et al. 2023). On such small scales, even though the very intricate physics yields highly non-Gaussian random fields, low-order statistics such as number counts or univariate statistics and two-point correlation functions are the measures routinely adopted to characterise the spatial distribution of matter fields (Scannapieco et al. 2006). Still, the exquisite quality and spatial coverage of current and forthcoming data sets are pushing towards the investigation of the relationship between the local and larger scale dynamics, including possible signatures of feedback physics. This improvement necessarily opens to the use of more sophisticated statistical analyses.

Since the advent of cosmological surveys, many higher order statistics have been used extensively to study the geometry and topology of the matter density field on large scales, from a few megaparsecs up to gigaparsec. Among them, Minkowski functionals (MFs) are of particular interest as they offer a complete and global description of the morphology

of spatial patterns, largely employed in several domains of natural sciences and for image analysis (e.g. Mecke & Stoyan 2002; Mantz et al. 2008). Using integral rather than differential expressions, they are therefore robust against small-scale spatial fluctuations and have a very simple interpretation in two and three dimensions (Adler 1981). Introduced in cosmology by Mecke et al. (1994), MFs supplemented the seminal studies by Gott et al. (1986) and Coles et al. (1993) on the topology of the large-scale structure (LSS) based on the genus curve and Euler-Poincaré characteristic (see also Melott et al. 1989; Gott et al. 1990; Park & Gott 1991; Mecke & Wagner 1991; Park et al. 1992, 2001; Vogeley et al. 1994a; Coles et al. 1996; Schmalzing & Buchert 1997; James et al. 2009; Park & Kim 2010; Zunckel et al. 2011; Choi et al. 2013). Moreover, MFs include the percolation analysis (Shandarin 1983; Shandarin & Zeldovich 1983; Yess et al. 1997; Shandarin & Yess 1998; Colombi et al. 2000) and the void probability function (e.g. White 1979; Vogeley et al. 1994b). They are also intimately related to alpha-shapes and Betti numbers (van de Weygaert et al. 2011; Park et al. 2013; Pranav et al. 2019; Feldbrugge et al. 2019).

Encompassing the two-point correlation functions (e.g. Kerscher et al. 1996, 1998), MFs have been used in cosmology to investigate the primordial non-Gaussianity of the matter field producing CMB radiation (e.g. Schmalzing & Gorski 1998; Schmalzing et al. 2000; Hikage & Matsubara 2012;

<footnote n="1">¹ Aix Marseille Univ., CNRS, CNES, LAM, IPhU, Marseille, France</footnote>

<footnote n="2">² Observatoire Astronomique de Strasbourg, Université de Strasbourg, CNRS, UMR 7550, 67000 Strasbourg, France
e-mail: kraljic@unistra.fr</footnote>

<footnote n="3">³ Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="4">⁴ University of the Western Cape, Bellville, Cape Town 7535, South Africa</footnote>

<footnote n="5">⁵ South African Astronomical Observatories, Observatory, Cape Town 7925, South Africa</footnote>

<footnote n="6">⁶ Institut d'Astrophysique de Paris, CNRS and Sorbonne Université, UMR 7095, 98 bis Boulevard Arago, 75014 Paris, France</footnote>

<footnote n="7">⁷ IPhT, DRF-INP, UMR 3680, CEA, L'Orme des Merisiers, Bât 774, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="8">⁸ Korea Institute for Advanced Study, 85 Hoegiro Dongdaemun-gu, Seoul 02455, Republic of Korea</footnote>

<footnote>* Corresponding author; carlo.schimd@lam.fr</footnote>
Ducout et al. 2013; Modest et al. 2013; Munshi et al. 2013; Buchert et al. 2017) and affecting the galaxy distribution on large scales (e.g. Hikage et al. 2006). They have also been used to study the gravitational non-Gaussianity of the late-time galaxy density field for cosmography applications (e.g. Blake et al. 2014; Wiegand et al. 2014; Fang et al. 2017; Sullivan et al. 2019; Liu et al. 2020; Appleby et al. 2022) and the projected dark-matter field from weak-lensing convergence maps (e.g. Kratochvil et al. 2012; Petri et al. 2013), as well as in extragalactic physics to assess the galaxy assembly history (Hikage et al. 2003; Einasto et al. 2011) and the morphology of the density field during and after the epoch of reionisation (e.g. Gleser et al. 2006; Yoshiura et al. 2017; Chen et al. 2019; Spina et al. 2021).

All the aforementioned applications concerned the morphology of total or baryonic (biased) matter density fields. The question is whether the geometric and topological content of other physical fields that characterise the thermodynamical and chemical properties of the gaseous component offer relevant and complementary information on the process of galaxy formation in general and on the impact of feedback on a larger scale in particular. These fields and their associated observables, such as the X-ray brightness (e.g. Gilli et al. 2007; Ji et al. 2020), Sunyaev-Zel'dovich Compton-$y$ maps (e.g. Bleem et al. 2022; Yang et al. 2022), Lyman-$\alpha$ absorptions (e.g. Peirani et al. 2014; Lee et al. 2018; Japelj et al. 2019; Kraljic et al. 2022), or the CO and $K$-band luminosities (e.g. Bolatto et al. 2013; Garratt et al. 2021) are usually considered on galactic or galaxy clusters scales, since they are associated with local hydrodynamical and feedback processes. The properties of MFs make them privileged statistics to evaluate whether their small-scale effect manifests itself on larger scales and with what delay; namely, they are used to assess how and when the physics of galaxies and circumgalactic medium affect the physics of the intergalactic medium. Indeed, MFs are robust to small-scale fluctuations, so their signal is plausibly measurable even when integrated on scales larger than a few megaparsecs. When compared with suitable observables, this can in turn put additional valuable constraints on often poorly understood baryonic processes, such as stellar and AGN feedback. The latter, namely the energy release from black holes, is expected to have a significant impact on the properties of galaxies, galaxy groups, and galaxy clusters, but also on their outskirts in the circumgalactic medium and beyond, as hot gas bubbles percolate (e.g. Fabian 2012; Sorini et al. 2022; Ayromlou et al. 2023). Indeed, the sources are mainly distributed along the cosmic web, so hot baryons are expected to percolate faster along filaments; the percolation process should then be boosted when the ionisation front enters an already ionised region.

Large-scale cosmological hydrodynamic simulations such as HORIZON-AGN (Dubois et al. 2014), ILLUSTRIS (Vogelsberger et al. 2014), EAGLE (Schaye et al. 2015), ILLUSTRISTNG (Pillepich et al. 2018), SIMBA (Davé et al. 2019), EXTREME-HORIZON (Chabanier et al. 2020), and HORIZON RUN 5 (Lee et al. 2021) are well suited to assess the impact of feedback processes on the large-scale morphology of the matter fields. Among them, the cosmological hydrodynamic simulation suite SIMBA stands out, as it implements several recipes of stellar and AGN feedbacks on top of hydrodynamical interactions; thus, monitoring the gas temperature, pressure, neutral atomic and molecular hydrogen density, and metallicity.

Here, SIMBA is used to address the question whether the MFs of thermodynamical fields can be used to distinguish between various feedback processes, in particular stellar feedback from various components of AGN feedback such as AGN winds, jets,

and X-ray heating. The global MFs of the excursion sets of the 3D gas fields are analysed as a function of the field value (or threshold) and redshift from $z=5$ to $z=0$, separately considering the effect of feedback recipes. Lacking motivated models of the spatial statistics, the analysis is limited to qualitative considerations; some quantitative summary statistics extracted from the Minkowski curves are then considered to describe the morphological time evolution of gas fields, beyond the evident observation that they are not Gaussian. The investigation of observed fields, usually projected on the celestial sphere (therefore requiring MFs in two dimensions) or in the 3D redshift space (and thereby affected by redshift-space distortions) is left for future studies.

The purpose of this study is to set up a physical framework in which to address such questions as how feedback impacts the geometry of the intergalactic medium (IGM); how the thermodynamic state and chemical content of the gas impact its topology; how the large-scale morphology is correlated with the major epochs of the cosmic evolution of star formation and active galactic nucleus feedback; whether there is any delay between the phenomena on large and small scales; and, finally, whether this delay simply reflects a propagation time or whether there are non-trivial percolation effects related to the known anisotropy of the cosmic web.

The remainder of this paper is structured as follows. The SIMBA suite and the mathematics of MFs are introduced in Section 2. The analysis of MFs as a function of redshift, separately for individual feedback processes is presented in Section 3. The overall picture is discussed in Section 4. Finally, Section 5 presents our conclusions.

## 2. Method

### 2.1. The SIMBA suite: baryonic physics on large scales

The SIMBA simulation is fully described in Davé et al. (2019); therefore, here, we provide only its summary and focus on its features relevant to this study. SIMBA was run with a modified version of the gravity plus hydrodynamics solver (Hopkins 2015), employing the GADGET-3 tree-particle-mesh gravity solver (Springel 2005) and a meshless finite mass solver for hydrodynamics.

Radiative cooling and photo-ionisation heating models make use of the GRACKLE-3.1 library (Smith et al. 2017) accounting for metal cooling and non-equilibrium evolution of primordial elements. A spatially uniform UV ionising background follows the Haardt & Madau (2012) model, modified to account for self-shielding based on Rahmati et al. (2013). The neutral hydrogen (HI) content of gas particles is modelled self-consistently.

SIMBA implements star formation based on molecular hydrogen ($H_2$) in a Jeans mass-resolving pressurised interstellar medium (ISM; with the hydrogen density $n_H \ge 0.13 \text{ cm}^{-3}$), following Davé et al. (2016). The computation of the $H_2$ fraction is based on the Krumholz & Gnedin (2011) prescription. The chemical enrichment model tracks eleven different elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe) from type Ia and II supernovae (SNe), and asymptotic giant branch (AGB) stars, following the yield tables of Iwamoto et al. (1999), Nomoto et al. (2006), and Oppenheimer & Davé (2006), respectively. SIMBA also tracks dust growth and destruction for each individual element on the fly (see Li et al. 2019 for a detailed investigation of the dust model, and Donevski et al. 2020 for comparison with observations). The stellar feedback is modelled using decoupled galactic winds, namely, with the hydrodynamics turned off in the
<table><caption>**Table 1.** Summary of main ingredients of SIMBA runs.</caption><thead><tr><th></th><th>Fiducial</th><th>NoX</th><th>NoJet</th><th>NoAGN</th><th>NoFeedback</th></tr></thead><tbody><tr><td>Stellar feedback</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td></td></tr><tr><td>AGN winds</td><td>✓</td><td>✓</td><td>✓</td><td></td><td></td></tr><tr><td>Jets</td><td>✓</td><td>✓</td><td></td><td></td><td></td></tr><tr><td>X-ray heating</td><td>✓</td><td></td><td></td><td></td><td></td></tr></tbody></table>

**Notes.** 'Fiducial' model includes all feedback processes, 'NoFeedback' model includes all the hydrodynamical but feedback processes. See Sect. 2.2.

winds until they leave the ISM, metal-loaded kinetic two-phase galactic winds, and with 30% of wind particles being ejected with the temperature set by the kinetic energy of the wind subtracted from the supernova energy.

SIMBA includes black hole (BH) particles, employing a specific two-mode accretion model. Hot gas ($T > 10^5$ K) is accreted in a spherically symmetric fashion following the Bondi (1952) formula, while cold gas accretion is modelled via a torque-limited sub-grid prescription describing the response of gas inflows near the BH to angular momentum loss due to dynamical instabilities (Hopkins & Quataert 2011; Anglés-Alcázar et al. 2017). As Bondi accretion models gravitational capture from a dispersion-dominated medium, it is more appropriate for hot gas. The torque-limited accretion mode is instead well suited for the growth of BHs in rotationally supported disks. This unique combination of BH accretion modes determines the implementation of feedback from active galactic nuclei (AGNs) in the form of two-mode kinetic feedback, the so-called radiative mode and jet mode feedback:

- Radiative mode: BHs with high accretion rates (above 0.2 times Eddington rate) and mass above $10^{7.5} M_\odot$ eject material in ~1000 km/s winds without changing its temperature. This is consistent with observations of ionised multi-phase gas outflows (Perna et al. 2017).

- Jet mode: As the BH accretion rate drops below 0.2 of the Eddington rate, jet feedback mode turns on and is fully achieved below 0.02. Gas is ejected at much higher velocities compared to the radiative mode, with a velocity increment proportional to the logarithm of the inverse of the accretion rate and capped at 7000 km/s. Another difference related to the radiative mode is the increase of temperature of the ejected particles in the jet mode, consistently with observations (Fabian 2012).

Ten percent of the material accreted into the central region is assumed to fall onto the BH, and the gas elements are immediately ejected in the purely kinetic and bipolar way (namely with zero opening angle w.r.t. the angular momentum of the inner disk) according to these two modes. In addition, SIMBA includes a third AGN-driven feedback channel:

- X-ray heating: X-ray radiation pressure feedback is activated only in galaxies with low cold gas content and when the jet mode is active. The effect of this X-ray radiation pressure feedback is to push outwards the gas surrounding the accretion disk based on the high-energy photon momentum flux generated in the black hole accretion disk. X-ray heating is proportional to the inverse square of the distance of the gas element with respect to the BH. Its implementation broadly follows the model of Choi et al. (2012) and works as follows. Non-ISM gas ($n_H < 0.13 \text{ cm}^{-3}$) is heated by directly increasing its temperature; whereas for ISM gas, one-half of the

X-ray energy is applied kinetically to give the gas particles a radial outwards kick and the second half is added as heat.

The combination of these three forms of feedback produces a population of quenched and star-forming galaxies and their black holes in good agreement with local studies of black hole properties (Thomas et al. 2019), with galactic size-mass relation and radial profiles (Appleby et al. 2020), and with 1.4 GHz radio luminosities (Thomas et al. 2021). Overall, while radiative AGN feedback only mildly affects the galaxy properties, the jet mode is mainly responsible for the quenching of galaxies. The X-ray feedback plays globally a small, but important role in suppressing residual star formation, and seems to be crucial for reproducing the observed centrally suppressed specific star formation rate and gas profiles of star-forming and green valley galaxies. These ingredients were included to improve the model of galactic-scale processes but potentially leave larger-scale signatures, which this paper will quantify using MFs.

## 2.2. SIMBA runs

To examine the effect of stellar feedback and different types of AGN feedback on the global geometry and topology of various gas fields, we use five runs of the SIMBA suite in $50 h^{-1}$ Mpc comoving volumes with $512^3$ gas elements and $512^3$ dark matter particles with the mass resolution of $1.82 \times 10^7 M_\odot$ for gas and $9.6 \times 10^7 M_\odot$ for dark matter particles. Cosmology adopted in SIMBA is a standard ΛCDM compatible with results from Planck Collaboration XIII (2016), namely $\Omega_m = 0.3$, $\Omega_\Lambda = 0.7$, $\Omega_b = 0.048$, $H_0 = 68 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\sigma_8 = 0.82$ and $n_s = 0.97$.

Five different runs of the SIMBA suite correspond to different variants of AGN feedback adopted following a strategy where one input physics element is turned off at a time, to which a model without AGN feedback, and without stellar feedback are added, as summarised below (see also Table 1):

- Fiducial – denotes a run with all forms of feedbacks,

- NoX – denotes a run with X-ray AGN feedback turned off, but including radiative and jet mode AGN feedback and stellar feedback,

- NoJet – denotes a run with both X-ray and jet AGN feedback turned off, but including AGN winds, namely radiative AGN feedback, and stellar feedback,

- NoAGN – denotes a run without any form of AGN feedback, but including stellar feedback,

- NoFeedback – denotes a run without any feedback.

Figure 1 offers an initial insight into the impact of the different feedback mechanisms on the six gas fields, showing two-dimensional sections $8 h^{-1}$ Mpc thick at $z = 0$.

## 2.3. Minkowski functionals

Minkowski functionals (MFs) are integral measures that generalise the notion of volume accounting for the size, shape (geometry) and connectivity (topology) of bodies or spatial patterns fully characterising their morphology, namely being subadditive, invariant under Galilean transformations, and continuous (Hadwiger 1957). Dealing with random fields like the density, temperature, pressure, and metallicity of the gaseous component of the LSS, it is natural to consider their excursion set or isocontours, $C_\theta = \{x \in D | f(x) \ge \theta\}$, where f denotes a generic field defined over a spatial domain D and $\theta$ the threshold (see Figures 2 and 3 for illustration). In three dimensions, the MFs of $C_\theta$ correspond to its volume (V), surface area (A), integral mean curvature (H), and integral Gaussian curvature (G), the two latter measuring the extrinsic and intrinsic curvatures,
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Visualisation of full box-size sections (8 h⁻¹ Mpc thick and 50 h⁻¹ Mpc wide in both directions) of the unsmoothed 3D temperature, pressure, total density, HI density, H₂ density, and metallicity fields (from top to bottom) at z = 0, for different SIMBA models progressively excluding feedback mechanisms (from left to right; see Table 1). It clearly shows a high level of morphological diversity. It also illustrates how individual feedback processes operating on galactic scales leave different morphological imprints on the gas fields at large scales; note for instance the influence of jets (second column) on the morphology of the T, P, nH₂, and Z fields, but not on n and nHI fields. The corresponding 3D smoothed fields at z = 0 and z = 5 as used for the computation of the MFs are shown in Appendix B as well as Figures B.1 and B.2. The purpose of this paper is to compress the information content of all maps into a few numbers extracted from the MFs of their excursion and quantify how feedback impacts them across cosmic time.</figcaption></figure>

respectively, of the body integrated over its surface. Owing to the Gauss-Bonnet theorem, the Gaussian curvature is proportional to the Euler characteristic $\chi = G/4\pi$ and linearly related to the

genus $g = 1 - \chi/2$ of $C_\theta$, both usually adopted as more directly related to the counts of connected regions, tunnels, and cavities of the excursion-set, or maxima, minima, and saddle points of
<figure><img src="image_3.png" /><figcaption>Fig. 2. Excursion-sets of the gas temperature for the fiducial model at z = 5 for different values of the threshold: whole field (top left), T ≥ 6 × 10³ K (top middle), T ≥ 8 × 10³ K (top right), T ≥ 1 × 10⁴ K (bottom left), T ≥ 3 × 10⁴ K (bottom middle) and T ≥ 5 × 10⁴ K (bottom right). For illustrative purposes, 3D maps correspond to fields before applying additional Gaussian smoothing of 1.56 h⁻¹ Mpc. Low values of the threshold reveal cold isolated cavities, at intermediate thresholds an interconnected filamentary structure emerges, while higher thresholds delimit hot regions forming isolated clumps. The Minkowski functionals computed for these thresholded 3D maps quantify the complex morphology of the underlying temperature field. See Figure 3 for analogous excursion-sets at z = 0.</figcaption></figure>

the thresholded field $f(\mathbf{x})$, accounting for the topology of $C_\theta$. Analytical expressions exist for simple configurations such as triaxial ellipsoids, cylinders, or trivial unions thereof describing the idealised isocontours of isolated or merging clusters with filaments or bridges (Schimd & Sereno 2021) and for the expectation value of MFs of Gaussian and Gaussian-related random fields (see Appendix A).

Here, we adopt the normal based on the so-called intrinsic volumes ($V_0, V_1, V_2, V_3$) = ($V, A/6, H/3\pi, \chi$) and compute the MFs as the spatial average of Koenderink invariants using the code presented in Schmalzing & Buchert (1997), which yields the MFs' densities as a function of the threshold, or Minkowski curves, as:

where $dV(\mathbf{x})$ and $dA(\mathbf{x})$ are the elementary volume and area of the excursion sets of the six random fields $f = T, P, n_{gas}, n_{HI}, n_{H_2}, Z$ (hereafter referred as f-fields), which have boundary surface $\partial C_\theta$ with principal local curvature radii $R_{1,2}(\mathbf{x})$. Also, $\mathcal{V}$ is the total sample volume, namely the simulation box. The fields are sampled at 128³ Cartesian regularly spaced lattice points using a cloud-in-cell kernel, to limit the computational load, and folded with a Gaussian kernel with standard deviation corresponding to a resolution of 4 lattice points, equivalent to about 1.56 h⁻¹ Mpc, largely serving the Nyquist limit. Dealing with simple (periodic) boundaries, the results are consistent with those obtained using the Crofton formula (see the discussions in Schmalzing et al. 1999; Hikage et al. 2003).

$$v_0(\theta) = \frac{1}{\mathcal{V}} \int_{C_\theta} dV(\mathbf{x}), \qquad (1a)$$

$$v_1(\theta) = \frac{1}{6} \frac{1}{\mathcal{V}} \int_{\partial C_\theta} dA(\mathbf{x}), \qquad (1b)$$

$$v_2(\theta) = \frac{1}{6\pi} \frac{1}{\mathcal{V}} \int_{\partial C_\theta} dA(\mathbf{x}) \left( \frac{1}{R_1(\mathbf{x})} + \frac{1}{R_2(\mathbf{x})} \right), \qquad (1c)$$

$$v_3(\theta) = \frac{1}{4\pi} \frac{1}{\mathcal{V}} \int_{\partial C_\theta} dA(\mathbf{x}) \frac{1}{R_1(\mathbf{x})R_2(\mathbf{x})}, \qquad (1d)$$

# 3. Morphological analysis

## 3.1. Guidelines

Figures 4–9 which are discussed in the following subsections show the four MFs per unit volume, $v_\mu$ ($\mu$ = 0, 1, 2, 3), of the excursion-set of the three-dimensional temperature, pressure, density (total, HI, and H₂), and metallicity fields for all the models, as a function of the threshold and at different redshift. Each row focuses on the effect of one component at a time, namely X-ray heating, jets, AGN winds, and stellar feedback (from top
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Excursion set of temperature map for the fiducial model at $z=0$ for different values of the threshold: whole field (top left), $T \ge 7 \times 10^2$ K (top middle), $T \ge 10^4$ K (top right), $T \ge 10^5$ K (bottom left), $T \ge 10^6$ K (bottom middle) and $T \ge 10^7$ K (bottom right). As for Figure 2, these maps do not include the additional Gaussian smoothing of $1.56 h^{-1}$ Mpc and adopt the same colour code.</figcaption></figure>

to bottom), by comparing each time the two models (solid and dotted lines on the top panels) differing by this same component. That is, the impact of individual feedback mechanisms is assessed by the differences between MFs densities (quoted as $\Delta$ in the bottom sub-panels) as follows:

- impact of X-ray heating: $v_\mu$[Fiducial] – $v_\mu$[NoX],

- impact of jets: $v_\mu$[NoX] – $v_\mu$[NoJets],

- impact of AGN winds: $v_\mu$[NoJets] – $v_\mu$[NoAGN],

- impact of stellar feedback: $v_\mu$[NoAGN] – $v_\mu$[NoFeedback].

The full model accounting for all the feedback mechanisms and the simplest model with no feedback (only hydrodynamics) correspond to the solid line in the top panels and to the dashed line in the bottom panels, respectively. In each panel, symbols pinpoint the MFs of the excursion set corresponding to the mean value of the field, $\tilde{f}(\mathbf{x})$; we will refer to these domains as $\tilde{f}$-regions.

The Minkowski curves of the excursion sets $C_\theta$ have similar qualitative trends regardless of the gas property one does consider. The following outline guides their interpretation:

- The volume filling fraction $v_0$ is a monotonically decreasing function of the threshold $\theta$, tending to zero if the highest values of the field $f$ are concentrated in point-like regions. Likewise, $1 - v_0$ measures the volume fraction of the complementary domain $D \setminus C_\theta = \{\mathbf{x} \in D | f(\mathbf{x}) < \theta\}$. We note that for the fixed value of the threshold, $v_0$ increases (decreases) with time when the corresponding domains expand (contract). When approximated by a step function $H(\theta_* - \theta)$, like the temperature field at high redshift, $v_0(\theta)$ accounts for a homogeneous (single-phase) field with characteristic value $f(\mathbf{x}) = \theta_*$. Correspondingly, the density area $v_1$ is close to a Dirac-delta, indeed $v_1(\theta) \approx dv_0(\theta)/d\theta$. Consistently, a multi-phase field with sharp domains is accounted for by a piecewise constant function, namely, $v_0(\theta) = \sum_i v_{0,i} H(\theta_i - \theta)$ with $\sum_i v_{0,i} = 1$, the $i$th phase occupying a fractional volume $v_{0,i}$.

- According to the isoperimetric inequality (Schmalzing et al. 1999), the surface area of fixed volume is minimum for a sphere, namely, a wrinkly surface has more area than a smooth one. All the excursions-sets with threshold different from the maximum point of $v_1(\theta)$ have therefore on average a more regular surface than domains with maximum $v_1$. Note also that the apparent skewness of $v_1(\theta)$ is due to the logarithmic scale of the threshold and typically grows with time.

- The integrated-mean-curvature density $v_2(\theta)$ is positive (negative) for domains $C_\theta$ that are on average convex (concave), which are dominated by lumps (cavities). The opposite is true for the complementary domains $D \setminus C_\theta$.

- The Euler characteristic is the alternate sum of Betti numbers, counting the number of topologically invariant domains, namely: connected regions, tunnels, and cavities. Consistently, the $v_3(\theta)$ curve has a characteristic M-shape: it is positive for small values of the threshold $\theta$, which are typically smaller than the mean value $\bar{\theta}$ (marked by symbols), and attains a local maximum that measures the largest possible number of cavities per unit volume in the excursion set; it becomes negative for intermediate values of $\theta$, with
<figure><img src="image_5.png" /></figure>

**Fig. 4.** Temperature. Minkowski curves (columns) of gas temperature for different models (top to bottom) and redshift (colours). Each row shows a comparison of two models differing by one ingredient at a time (indicated in the plot titles), probing its impact on Minkowski functionals. In each row, the upper panel shows the two models (solid and dashed lines) and the bottom panel shows their differences, $\Delta$. Filled coloured circles and plus symbols indicate the mean value of the field at each redshift for two models shown with solid and dashed lines, respectively. First row: Fiducial model is compared to the NoX, showing the effect of X-ray heating. Second row: NoX model is compared to NoJet, showing the effect of jets. Third row: NoJet model is compared to NoAGN, highlighting the effect of AGN winds. Fourth row: NoAGN model is compared to NoFeedback, showing the effect of stellar feedback. The morphology of the T-field is most strongly impacted by AGN jets at low redshift ($z \lesssim 2$), but also by the X-ray heating at both high ($z \sim 4-5$) and low redshift ($z \lesssim 2-3$), though to a somewhat lesser extent. Interestingly, the strongest impact of the stellar feedback (with a strength comparable to X-ray heating) is seen near the epoch of cosmic noon $z \sim 2-3$.

a minimum counting (in absolute value) the largest number density of tunnels in the excursion set, corresponding to a ‘sponge-like’ topology; and it is again positive for larger values of the field, attaining a second (often global) maximum that indicates the largest number density of isolated lumps in the system. When the negative peak occurs for $\theta > \bar{\theta}$ or $\theta < \bar{\theta}$

(for instance, as in the T-field at low redshift for the NoAGN model with stellar feedback only; see Figure 4, bottom-right panel) the topology of the field is respectively interpreted as ‘bubble-like’ (or ‘Swiss-cheese-like’), namely dominated by cavities, or ‘meatball-like’, namely composed of isolated regions (see Coles et al. 1996; Choi et al. 2013). The right-
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Pressure. Minkowski curves for the gas pressure or *P*-field (in units of Boltzmann constant, noted *k*), analogous to Figure 4. Stellar feedback is the first cause of geometrical and topological change of the *P*-field set by gravitational and hydrodynamical interactions, in particular at high redshift (z ≳ 2). At low redshift (z ≲ 1.5), it is AGN winds and AGN jets that have a dominant impact on the morphology of the *P*-field; see Sect. 3.3.</figcaption></figure>

most zero of $v_3(\theta)$ defines an estimation of the threshold for continuous percolation (Mecke & Wagner 1991).

- Gaussian random fields, the usual reference for the LSS, are described by analytical MFs' mean values depending on the variances of the field and its gradient (Tomita 1986), with even $v_1(\theta)$ curve peaked at $\theta = \bar{\theta}$, odd $v_2(\theta)$ with equal absolute amplitude of extrema and vanishing for $\theta = \bar{\theta}$, and even $v_3(\theta)$ with a negative minimum at $\theta = \bar{\theta}$ and two local maxima, the amplitude of $v_{1,2,3}$ being determined by the first spectral moment of the field (equivalent to the rms-

variance of its gradient $\nabla f$). Weakly non-Gaussian and log-normal random fields also admit analytical MFs' expectation values, with coefficients depending respectively on the generalised skewness parameters and on the variance of the field (Matsubara 2003; Hikage et al. 2003, 2006; Gay et al. 2012; Matsubara et al. 2022). More details are given in Appendix A. Not surprisingly, the MFs of the gas fields considered in this study are strongly non-Gaussian, as indicated by the strong deviation from the M-shape around the minimum of $v_3(\theta)$, especially at low redshift as a result of non-
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Total density. Minkowski curves for the total gas number density or *n*-field, analogous to Figure 4. Until *z* ≈ 1 stellar feedback modifies morphology by the largest amount, progressively weakening until *z* ≈ 0.5; at later time AGN feedback mechanisms become the main driver shaping the *n*-field.</figcaption></figure>

linear gravitational evolution and shot noise (in observed
fields, additional deviations are due to redshift space dis-
tortions and, if any, primordial non-Gaussianities). Apart
from some specific deviations and especially at high red-
shift, *z* > 4, their profile resembles that of log-normal ran-
dom fields.

Qualitatively, we can recognise the following spatial patterns
from the values of MFs densities:

- Isolated regions with field values larger than in the surround-
ing environment, *f* > *f*<sub>env</sub> (e.g. hot bubbles in cold or warm
IGM), have small *v*<sub>0</sub>, large (small) *v*<sub>1</sub> if they have smooth
(wrinkly) boundary surface, *v*<sub>2</sub> > 0, and large *v*<sub>3</sub>;

- Isolated bubbles with *f* < *f*<sub>env</sub> (e.g. cold bubbles in warm IGM) have large values for *v*<sub>0</sub>, surface density *v*<sub>1</sub> like for iso-

lated regions, while $v_2 < 0$, and large values for $v_3$ like for
isolated regions;

- Network of filaments with $f > f_{\text{env}}$ (e.g. hot, dense, or metal-rich filaments alimented by local stellar activity) have large (small) $v_0$ if on average filaments are thick (thin), $v_1$ as above, $v_2 > 0$ and smaller for cylindrical shape (the largest radius of curvature diverges), and $v_3 \approx 0$ if the network is percolating;

- Network of filaments with $f < f_{\text{env}}$, namely, an excursion set $C_{\theta}$ dominated by tunnels (e.g. cold filaments in warm IGM, or metal-poor filaments in which metals are expelled by supernovae or stellar winds) have opposite volume filling fraction than in the previous case, namely, large (small) $v_0$ if filaments are on average thin (thick), $v_1$ as before, $v_2 < 0$ and larger for cylindrical shape, and $v_3 < 0$ and large in
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** HI density. Minkowski curves for the neutral atomic hydrogen number density, or $n_{HI}$-field, analogous to Figure 4. They are similar to the Minkowski curves of the *n*-field at *z* > 2.5, and likewise modified at a later time by stellar feedback but in the opposite way; see Sect. 3.4.</figcaption></figure>

absolute value, and vanishing if the network is percolating as before. The directions of the filamentary network and the anisotropy of the excursion set cannot be captured by MFs, but by higher-rank Minkowski valuations (e.g. Beisbart et al. 2002).

Finally, we note that one probes the morphology of fields that are effectively smoothed twice, first by grid sampling and then by the Gaussian filter needed for the calculation of the covariant derivatives in Koenderink invariants. The resulting spatial resolution of the excursion sets is about $1.56 \, h^{-1} \, \text{Mpc}$, largely coarser than the original resolution of SIMBA. The numerical accuracy of MFs' densities is also limited for excursion sets with extreme values of the field at odds with their environment, as expected in very sparse or non-resolved regions. The global morphological

analysis of fields in domains of a size typical of galaxy groups or smaller and in very diffuse domains is computationally demanding and, therefore not addressed in this study. We note that in the following subsections, all the quoted values of fields must therefore be considered carefully, namely, as spatially-averaged values. A more detailed analysis is presented in Appendix C.

## 3.2. Temperature field

Figure 4 shows the MFs of the excursion sets of the smoothed T-field, whose three-dimensional representations are illustrated in Figures 2 and 3.

When all the feedback mechanisms are included (solid lines in first-row panels), at redshift $z > 4$ the volume filling fraction
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** H₂ density. Minkowski curves for the molecular hydrogen number density or *n*<sub>H₂</sub>-field, analogous to Figure 4. Its evolution is very similar to that of the total gas density; see Sect. 3.4.</figcaption></figure>

$v_0(T)$ is a step-like function with the transition occurring around warm regions with $T \approx 10^4$ K, indicating that at that epoch the universe was approximately in a single phase, with hot regions occupying a tiny fraction of the total volume. The main effect of setting this temperature is photo-ionisation heating within the IGM. At low redshifts, due to the dropping meta-galactic photo-ionisation rate, gas can reach temperatures down to $\sim 10^3$ K. However, in the case of full SIMBA physics, the IGM is heated strongly by AGN feedback (Christiansen et al. 2020), resulting in a much smaller fraction of photo-ionised IGM. Another effect is that SIMBA's pressure floor does not allow cooling in the interstellar medium below $10^4$ K, but the fraction of overall baryons

in the ISM is small in any case, particularly with full SIMBA physics.

Correspondingly, the surface density $v_1(T)$ is symmetric and very peaked around $10^4$ K, resembling a Dirac delta, indicating a single-phase T-field at high redshifts. At later times, hotter regions occupy larger and larger volume as time passes due to the combined effect of shock heating on the one hand, and star formation and AGN feedback on the other hand, where the latter typically dominates (e.g. Somerville & Davé 2015). The surface density curve $v_1(T)$ becomes and remains right-tailed until $z \sim 2$ and left-tailed afterwards, in favour of regions with temperature smaller than the mean value at late times, namely, for
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Metallicity. Minkowski curves for the gas metallicity or Z-field, analogous to Figure 4. AGN feedback mechanisms, especially jets, modify the morphology of regions with $Z \gtrsim 0.01 Z_\odot$ for $z < 1.5$, while at earlier time the morphology is maximally influenced by stellar feedback, especially in regions with $Z \lesssim 0.1 Z_\odot$; see Sect. 3.5.</figcaption></figure>

$T < \bar{T} \approx 10^6$ K at $z=0$ (marked by a circle). These regions are concave ($v_2 < 0$) and dominated by tunnels ($v_3 < 0$), especially at redshift $z=1$ for $T \gtrsim 10^4$ K and at $z=0.5$ for $T \gtrsim 10^5$ K, suggesting a network of cold filaments in a warm percolating IGM. Interestingly, at $z=0$ concave regions, namely tunnels and bubbles, with $T \lesssim 10^{5.5}$ K stop to negatively contribute to the Euler characteristic, which is almost vanishing; this is compatible with the late-time appearance of cold bubbles, as illustrated in Figure 3 (top middle panel).

Without AGN and stellar feedbacks, namely with only pure hydrodynamics (dashed line in bottom line panels), the morphological evolution of the temperature field is very different;

namely, it is qualitatively close to a log-normal random field at all redshifts, except $z=0$ (not shown; see Appendix A). MFs indicate regions with temperature slightly larger than $T = 10^4$ K as very special: besides marginal deviations at $z < 0.5$, at all times they fill about one-third of volume ($v_0 \approx 0.35$), have the maximum surface density, vanishing mean curvature, and minimum (negative) Euler characteristic, suggesting a network of cold filaments set by pure gravity and hydrodynamics.

Overall, AGN jets have the strongest impact on the morphology of the temperature field (second-row panels), as can be seen in the amplitude of $\Delta$. Starting at redshift $z=2$, namely ~3.3 Gyr after the Big Bang, they substantially change the
profile of MFs at all the values of temperature explored in SIMBA simulations: they increase by more than 50% the volume filling fraction, especially of regions with $T > 10^4$ K, skew the $v_1(T)$ curve by decreasing the surface area of cold ($T \approx 10^4$ K) regions by ~7% and increasing that of hot regions ($T \ge 10^6$ K) by almost a factor of 10, and almost reverse the convexity of $T > 10^4$ K regions. In the same epoch, AGN jets also strongly alter the relative abundance of clumps, tunnels, and bubbles. In particular, as the universe evolves with time, bubbles are restricted to regions with temperature $T \le 10^4$ K, without jets they would occupy also cooler regions (the first zero of $v_3(T)$ decreases with time, especially for $z > 1.5$, as shown by the dashed line in the second-row or solid line in third-row panels). Because of jets, at the same time, clumps form in regions with larger and larger temperatures (the second zero of $v_3(T)$ increases almost linearly with redshift excluding the redshift range $1 < z < 1.5$, when it is nearly constant). Accordingly, the temperature range concerned by negative values $v_3$, or tunnels, increases with time, reaching the minimum value at $z = 0$ for $T = 10^6$ K regions. Without jets, the topology of the $T$-field would be dominated by tunnels over a very narrow range of temperatures at $z = 5$ ($10^4$ K $\le T \le 10^{4.2}$ K) and a much broader range at $z = 0$ ($10^{3.3}$ K $\le T \le 10^5$ K); after $z = 2$, because of jets, the percolation threshold increases dramatically, up to two orders of magnitude at $z = 0$ from $T \approx 10^{4.5}$ K to $T \approx 10^{6.3}$ K, consistently with a warm-hot (ionised) IGM permeating the largest fraction of the total volume, with less than 10% occupied by hot regions, and very few cold cavities.

X-ray heating and stellar feedback have a secondary role with respect to AGN jets in shaping the morphology, by about the same relative amount with respect to their reference model. However, while X-ray heating operates in two distinct epochs at high and low redshift ($z > 4$ and $z < 1.5$), respectively on cold ($T \approx 10^4$ K) and hot ($T \ge 10^5$ K), stellar feedback acts at intermediate redshift ($1 < z < 2.5$) and on warm regions ($10^{3.5}$ K $\le T \le 10^5$ K). This is because X-ray feedback acts to remove cool gas from the central regions of galaxies with massive black holes, which can alter the morphology of ~$10^4$ K gas via heating or removal at high redshift, while at lower redshift it enables the jet feedback to be more effective at redistributing baryons from within galaxies into the hot IGM. Meanwhile, stellar winds have more moderate velocities than jets, so they remove cool gas from galaxies but are only able to heat it to warm-hot temperatures in the intergalactic medium.

AGN winds have the smallest effect on the morphology of the $T$-field, of the order of a few percent, limited to very low redshift ($z < 0.5$) for regions with temperatures around $10^{4.5}-10^{5.5}$ K. They non-trivially change the topology: AGN winds increase the value of $v_3$ and therefore the number of hot clumps ($T > 10^5$ K) at $z = 0.5$, and enhance the peculiar topology of regions with $10^{3.7}$ K $\le T \le 10^{4.4}$ K already set by jets, by increasing (decreasing) the number of tunnels in regions with temperature smaller (larger) than $T = 10^4$ K. In general, radiative AGN winds are seen to have very minor effects on galaxy properties (Davé et al. 2019), and this appears to be reflected in the global gas morphologies as well.

A complementary way to interpret the Minkowski curves is looking at Figure 4 from top to bottom and considering the excursion set defined by the mean temperature, hereafter dubbed $\bar{T}$-regions<fnref n="1" />. Focusing on the volume filling fraction, the $\bar{T}$-regions occupy less than half the volume when all the AGN and stellar

feedback mechanisms are included, with $v_0(\bar{T})$ decreasing from about 50% at $z = 5$ to 15–20% at $z = 2.5$ and rising back to about 50% at a later time. X-ray heating qualitatively does not change this trend (filled circles and crosses in the first row panels are almost superposed), modified by 10–13% at $z > 4$ for $T = 10^4$ K regions and by less than 10% at $z < 1.5$ in hotter regions. Jets strongly heat the gas, especially at $z < 2$, pushing the $\bar{T}$-regions to occupy more than 20% of the full volume. Stellar feedback and AGN winds do not change qualitatively the trend of $v_0(\bar{T})$ established by pure hydrodynamics, the former mainly affecting the volume filling fraction of $T = 10^{4.2}$ K regions from $z = 3$ down to $T = 10^{3.8}$ K regions at $z = 0.5$. AGN winds only increase $v_0$ of $T \approx 10^4$ K regions by about 3% and only at $z = 0$.

In the full model, the density of surface area of $\bar{T}$-regions, $v_1(\bar{T})$, is maximum at very high redshift ($z > 4$), consistent with a single-phase (uniform) fluid at mean temperature, and at very low redshift ($z < 0.5$), indicating that the corresponding excursion set is jagged, namely, non-spherical, disturbed by late-time or integrated energy injections. Indeed, as reminded in the Guidelines (see Sect. 3.1), spheres are the bounded regions with minimum surface for a fixed volume. Consistently, in these two epochs $v_2(\bar{T})$ and $v_3(\bar{T})$ indicate non-maximal convexity and minimum Euler characteristic, namely, spongy-like topologies dominated by tunnels. The morphological transition occurs between $z = 2.5$ and $z = 1.5$, the $\bar{T}$-regions with $\bar{T} = 10^5$ K at $z = 2$ having a smoother shape with the number of isolated regions being almost equal to the number of tunnels ($v_3(\bar{T}) \approx 0$); namely, the $\bar{T}$-regions percolate at the peak of the star formation rate. At that epoch, the temperature field is strongly non-Gaussian, with regions with $T < \bar{T}$ dominated by tunnels. This corresponds to a network of cold filaments surrounded by hot-ter gas resulting from the energy injection by X-ray heating and AGN jets as expected by the percolation of shock waves. As time goes by, the effect of jets, stellar feedback, and AGN winds (in this order) is to make cool regions with $T < \bar{T}$ less and less round as time goes on (larger surface area, negative integrated mean curvature, negative Euler characteristic), consistent with a spongy-like topology. Finally, by $z = 0$, these cooler regions become very sparse as most of the IGM becomes heated via AGN feedback and gravitational shock heating.

### 3.3. Pressure

Figure 5 shows the MFs of the $P$-field, showing a redshift evolution qualitatively opposite to that of the $T$-field, consistent with a universe that evolves towards a thermodynamical vacuum, namely, with lowering pressure. We note that the multi-phase nature of the gas components invalidates the assumption of a single equation-of-state relating $P$, $T$, and $n$ (see Appendix D). Therefore, the morphology of the $P$-field cannot be deduced from that of the $T$ and $n$ fields.

The Fiducial model (solid lines, first-row panels) manifests a smooth and weak evolution of the morphology until approximately $z = 1.5$, as indicated by a small translation of the Minkowski curves that essentially maintain their shape, followed by a sudden evolution occurring between $z = 0.5$ and $z = 0$ (the abrupt change of $v_0$ and $v_1$ curves at $P > 10^2 k_B \text{ cm}^{-3}$ K, especially visible at $z < 0.5$, is a numerical artefact). At very early times, between $z = 5$ and $z = 4$, the four Minkowski curves do not change, apart from a slight increase of $v_3$ for the highest pressure regions ($P > 10^5 k_B \text{ cm}^{-3}$ K). However, during this epoch the $\bar{P}$-regions tend to occupy only a small fraction of the total volume ($v_0(\bar{P}) \approx 0.2$ at $z = 5$ and $<0.05$ at $z = 4$), suggesting that apart from sparse, small, and convex (spherical)

<footnote n="1">¹ We propose this kind of analysis for the temperature field only, for illustrative purposes.</footnote>
high-pressure clumps most of the volume is low pressure with $P < \bar{P}$. For $z < 4$ and until $z \sim 1.5$ the fraction of the total volume occupied by high-pressure domains progressively decreases, the $\bar{P}$-regions and furthermore the highest pressure domains forming a network with meatball topology. Between $z=1$ and $z=0.5$ the volume occupied by the $\bar{P}$-field and its surface density continue to increase while its convexity decreases, which is compatible with the progressive appearance of tunnels as indicated by the negative Euler characteristic. In the last time step, for $z < 0.5$, all regions with $P > 10^{3.9} k_B \text{ cm}^{-3} \text{ K}$ disappear; regions with $P < \bar{P} = 10^{3.4} k_B \text{ cm}^{-3} \text{ K}$ occupy more than 60% of the volume and are essentially concave with $v_3 \approx 0$, compatible with either a network with an equal number of bubbles and tunnels or a homogeneous distribution; and regions with $10^{3.4} < P/k_B \text{ cm}^{-3} \text{ K} < 10^{3.9}$ form convex, isolated clumps.

Stellar feedback is the main cause of deviation from the morphology of the P-field set by the no feedback case (dotted lines in bottom-row panels). In the absence of any feedback, the geometry and topology of the P-field smoothly evolve over the full redshift range $0 < z < 5$ studied here, the MFs curves shifting from high to low values of P. Stellar feedback limits high-pressure domains with $P \ge 10^{4.9} k_B \text{ cm}^{-3} \text{ K}$ to occupy less than 15–20% of the total volume, quenching the time evolution of their shape and topology until $z \sim 2$. This is visible by comparing the MFs of the P-field, which move towards smaller values only for $z < 2$ when stellar feedback is active.

X-ray heating and jets have a quantitatively similar impact on the MFs, the former in the redshift range $1.5 < z < 3$ and mostly concerning the domains with $P \ge 10^{2.9} k_B \text{ cm}^{-3} \text{ K}$, the latter at $z < 1.5$ (not considering the highest pressure regions with $P/k_B \ge 10^{5.9} k_B \text{ cm}^{-3} \text{ K}$, very likely not resolved because of smoothing; see Appendix D). Among the AGN feedback mechanisms, winds have the largest impact on the global morphology but are limited to the lowest redshift, $z < 1$, changing $v_0$ by 20%, $v_1$ by more than 40%, $v_2$ by 50%, and $v_3$ by up to 100%. At higher redshifts ($z \ge 1.5$) AGN winds modify the geometry of the P-field approximately at the same level as the other AGN feedback mechanisms.

## 3.4. Density fields

### 3.4.1. Total density field

Figure 6 shows the MFs for the total gas number density field $n$. Overall, the redshift evolution of the total and individual density fields, namely neutral (HⅠ) and molecular (H₂) gas, is qualitatively similar to that of the P-field; as suggested by the volume filling fraction, high-density regions progressively disappear in agreement with the global expansion of the universe and the dilution of matter.

When all the feedback mechanisms are included (solid lines in top row panels), the global morphology at high redshift vaguely reminds that of a log-normal random field. Between $z=5$ and $z=0$, the domains with mean density or $\bar{n}$-regions (circles) occupy a fraction of volume that increases from 20 to 40%; starting from $z=1.5$, they have the largest surface density, and approximately at the same epoch their convexity begins to decrease at a lower rate than before, though always being positive. A network of tunnels defined by $\bar{n}$-regions, or equivalently a network of filaments with $n < \bar{n}$, is already in place at $z=5$; it grows until $z=3$ while the abundance of voids is almost constant, then it stops growing between $z=3$ and $z=1$ (the minimum of $v_3$, marked by circles is almost constant in this period); after $z=1$ it grows again, at the expenses of higher density

clumps ($n \approx 10^{-2} \text{ cm}^{-3}$) that occupy a very tiny fraction of the total volume, and at the expenses of voids that instead occupy the largest fraction of the volume.

The largest modification to the global morphology of the $n$-field is caused by stellar feedback (bottom-row panels). This is expected because, while AGNs (particularly jets) dominate the global feedback energy input, stellar winds displace a far larger amount of mass. Without stellar feedback (dotted lines) its morphological evolution would smoothly evolve from $z=5$ to $z=0$; the values of all the MFs but $v_2$ of $\bar{n}$-regions (marked by crosses) are constant over time, and their global topology is dominated by tunnels. With stellar feedback (solid lines), the time evolution of $\bar{n}$ and, correspondingly, the morphological evolution of $\bar{n}$-domains are instead quenched until $z \approx 1-1.5$, namely until the late-end of the star-formation-rate peak (Madau & Dickinson 2014). MFs suggest that at $z > 2$ stellar feedback limits the abundance and size of gas regions to lower values and reduces the number of tunnels, but does not modify their convexity ($v_2$). Below $z=1$, stellar feedback is not morphologically efficient anymore; the geometry and topology of the total density field evolve controlled by AGN feedback mechanisms, especially AGN winds, as also indicated by the MFs of $\bar{n}$-regions (circles in bottom-row panels).

Overall, AGN feedback mechanisms mildly modulate the geometry and topology of the gas density field set by pure hydrodynamics and stellar feedback until $z = 0.5-1$. Afterwards, they start tempering or reinforcing the existing morphology, competing with stellar feedback.

AGN winds and X-ray heating have a qualitatively similar influence on the MFs of regions with density $10^{-3} \le n/\text{cm}^{-3} \le 10^{-0.5}$. Both decrease the volume filling fraction by more than 10% and alter the convexity by ~40% in the same way. Nonetheless, AGN winds mainly operate at $z < 1$ and produce rounder overdensity regions with a smaller surface, while X-ray heating has the same effect also in the redshift range $1 < z < 2$.

Jets modify the volume filling fraction of density regions by only 5% but in a more subtle way: they increase $v_0$ of both overdense and underdense regions at $z > 2$, increase (decrease) $v_0$ of overdense (underdense) regions in the range $0.5 < z < 2$, and decrease $v_0$ of both overdense and underdense regions at $z < 0.5$. The effect on $v_1$ and $v_2$ is qualitatively and quantitatively similar to AGN winds.

As for the topology of the $n$-field, all the AGN feedback mechanisms have qualitatively and quantitatively similar impact on the topology of the density field, limited to $z < 3$. X-ray heating and jets decrease the number of tunnels with a density slightly larger than the mean value and more intensively increase the number of high-density clumps. AGN winds instead decrease the number of high-density clumps. Both jets and winds are moderately more efficient in shaping the topology of the density field at redshift $z < 0.5$.

The morphology of low-density excursion sets with $n \le 10^{-5} \text{ cm}^{-3}$ (not shown in figures) is worth commenting on. Focusing on the Fiducial model at $z=0$, these regions form a pattern made of very small and isolated domains with overall almost vanishing volume-filling-fraction (so that the complementary excursion set has $v_0 \approx 1$), very small surface density, and positive integrated mean curvature (correspondingly, the complementary excursion sets with $n > 10^{-5} \text{ cm}^{-3}$ have $v_1 \ge 0$ and $v_2 \le 0$), and very small and positive Euler characteristic density, $v_3 \ge 0$, namely, they form a diffuse ensemble of disconnected spheres. This morphology is qualitatively unchanged by feedback mechanisms and qualitatively similar for domains with larger density at an earlier time. Although restrained by the spatial smoothing of
the field, similar conclusions concern very dense excursion sets with $n > 10^{-1.5} \text{ cm}^{-3}$ at $z = 0$ ($n \ge 2 \text{ cm}^{-3}$ at $z = 5$).

### 3.4.2. HI density field

The MFs for the number density field of the atomic neutral hydrogen gaseous component are shown in Figure 7. Similarly to the total gas density field, when all the feedback mechanisms are included (solid lines in top row panels) the global morphology of the $n_{HI}$-field is qualitatively similar to that of a log-normal random field at very high redshift, $z = 5$, and evolves in a qualitatively similar way afterwards, though less linearly with redshift. Three epochs can be identified: (i) between redshift $z = 4$ and $z = 2.5$ the surface density, curvature, and number density of very dense clumps ($10^{-2} \text{ cm}^{-3} < n_{HI} < 10^{-0.5} \text{ cm}^{-3}$) increase, indicating an evolution compatible with a local gravitational collapse, like for the total gas $n$-field; (ii) from $z = 2.5$ to $z = 0.5$, unlike the $n$-field, the regions with similar and smaller density ($10^{-2.5} \text{ cm}^{-3} < n_{HI} < 10^{-1.5} \text{ cm}^{-3}$) become more spherical (smaller $v_1$, larger $v_2$) and more abundant (larger $v_3$); (iii) an abrupt change of the global morphology occurs at $z = 0.5$, the volume filling fraction $v_0$ being almost unchanged until this epoch like for the total gas $n$-field. However, unlike the latter, at $z = 0$ the mean HI field occupies half of the volume ($v_0(\bar{n}_{HI}) \approx 0.5$) and is on average flat ($v_2(\bar{n}_{HI}) \approx 0.5$), namely, it is much more Gaussian although the maxima of $v_3$ have again unequal amplitude.

The stellar feedback is the main source of disturbance, like for the $n$-field. However, it works in the opposite way and especially at $1 < z < 3$, not progressively at all redshifts: it does increase $v_0$ by up to 30%, especially for intermediate density regions; it increases (decreases) the surface density $v_1$ of regions with higher (lower) density than the critical value $n_{HI} \approx 10^{-2.5} \text{ cm}^{-3}$ by up to 70%; it does change the density of integrated mean curvature $v_2$, again with respect to the critical value $n_{HI} \approx 10^{-2.5} \text{ cm}^{-3}$; and increases the number of tunnels of excursions sets for a narrow range of thresholds, $10^{-2.5} < n_{HI}/\text{cm}^{-3} \approx 10^{-2}$.

The three AGN feedback mechanisms affect the morphology of the $n_{HI}$-field to a smaller extent and with similar amplitude. X-ray heating modifies the same densities as stellar feedback but with the opposite effect and is limited to $1.5 < z < 3$ and, suddenly, at $z = 0$. Among the AGN feedback mechanisms, it has the largest impact on the topology. AGN winds have an analogous effect but at a later time, in the redshift range $0.5 < z < 1.5$. Jets behave like stellar feedback, apart from a less trivial impact on the topology of the field and a very specific and strong impact on the geometry and topology of domains with $n_{HI} < 10^{-2.5} \text{ cm}^{-3}$ at $z = 0$, like X-ray heating.

### 3.4.3. H₂ density field

Figure 8 shows the MFs for the molecular hydrogen number density or $n_{H_2}$-field. Very likely because of its lower fraction, they are very similar to that of the total gas density (apart from the obvious rescaling, i.e. the shift along the horizontal axis because of the lower density) and the feedback mechanisms operate in a very similar way. Therefore, one can come to the same conclusions discussed above for the $n$-field.

### 3.5. Metallicity

Figure 9 shows the MFs for the gas metallicity or Z-field. The redshift evolution of the MFs, qualitatively similar to that of the

T-field shown in Figure 4 and at odds with that of the $n$-field shown in Figure 6, supports a global metal enrichment over time, with a morphological transition occurring around $z = 1.5$. Overall, the Z-field is composed by (i) small, spheroidal, metal-poor bubbles at early time that progressively disappear; (ii) a network of thin filaments with $Z \approx Z_\odot$; and (iii) a collection of small bubbles with metallicity increasing with time from $10^{-2.5} Z_\odot$ at $z = 5$ to $10^{0.5} Z_\odot$ at $z = 0$.

When all the feedback mechanisms are considered, since the mean metallicity grows with time, the domains with fixed Z continue to expand, occupying larger and larger volumes, with $\bar{Z}$-regions (symbols) filling about 30% of the volume until $z = 2.5$, then slightly shrink until $z = 1.5$, and strongly expanding afterwards up to 50% at $z = 0$. The surface density of metal-poor regions ($Z/Z_\odot \le 0.018$) decreases with time, while for regions with intermediate metallicity ($0.01 < Z/Z_\odot \le 1$) it attains a maximum and then decreases; at late time ($z = 0$) both attain negative $v_2$ and form a network dominated by tunnels. Metal-rich regions with $Z/Z_\odot \ge 1.8$ are instead characterised by $v_1$ increasing with time, and positive $v_2$ and $v_3$, namely, they form a network of sparse and small spherical metal-rich clumps.

For $z > 1.5$, the morphology of the Z-field is not set by feedback (dotted lines in bottom-row panels), and is maximally influenced by stellar feedback, especially regions with $Z/Z_\odot \le 0.1$. Their volume filling fraction of is reduced by 5 to 25%, and their surface density and mean curvature are substantially altered. This is expected since stellar winds are the main distributor of metals into the IGM at high redshifts.

For $z < 1.5$, all the feedback mechanisms coherently establish the morphology of regions with $Z/Z_\odot > 0.01$, and especially those with $Z/Z_\odot \approx 0.1$, inducing modifications that increase with time, qualitatively opposed to stellar feedback whose effect is largest at high redshift. During this epoch, AGN jets play the major role: they increase the volume filling fraction by up to 60% at $z = 0$, decrease (increase) the surface density of regions with $Z/Z_\odot < 0.89$ ($Z/Z_\odot > 0.89$) by 80% at $z = 0$, changing their mean curvature and topology, enhancing that of metal-rich regions with $Z/Z_\odot > 3.2$. Since the mass ejected in jets is small, much of this effect is expected arise due to the interaction of the jets with surrounding gas and pushing it outwards, which is consistent with findings that SIMBA's jets strongly evacuate halos (e.g. Appleby et al. 2020; Sorini et al. 2022). AGN winds have the same qualitative effects as jets, by a factor of approximately two times smaller, and very similar to stellar feedback. X-ray heating only affects the very late-time morphology of metal-rich regions, in a narrow range of metallicities peaked at $Z/Z_\odot \approx 3.2$ and at redshift $0 < z < 0.5$. This is because the direct effect of SIMBA's X-ray feedback is limited to the central regions of halos.

## 4. Discussion

### 4.1. Ranking of feedback mechanisms

Qualitative inspection of MFs amplitudes shown in Figures 4–9 (differences in the lower part of each panel) suggests the ranking of feedback effects reported in Table 2, which takes into account all the four MFs altogether (not only $v_0$, as Figure 1 would wrongly suggest) at $z \lesssim 1.5$ and $z \gtrsim 1.5$; note that it reflects global, integrated effect over all possible threshold values. The following discussion distinguishes the four feedback mechanisms in the same order shown in the aforementioned figures:

**X-ray heating.** None of the explored gas fields is impacted by X-ray heating in a dominant way over other feedback mech-
<table><caption>**Table 2.** Ranking of feedback mechanisms w.r.t. their impact on the global morphology of the gas fields at $z \lesssim 1.5$ and $z \gtrsim 1.5$ (first and second number).</caption><thead><tr><th></th><th><i>T</i></th><th><i>P</i></th><th><i>n</i><sub>gas</sub></th><th><i>n</i><sub>HI</sub></th><th><i>n</i><sub>H<sub>2</sub></sub></th><th><i>Z</i></th></tr></thead><tbody><tr><td>X-ray heating</td><td>2–2</td><td>3–2</td><td>4–2</td><td>2–2</td><td>4–2</td><td>3–4</td></tr><tr><td>Jets</td><td>1–1</td><td>2–3</td><td>3–3</td><td>4–3</td><td>3–3</td><td>1–2</td></tr><tr><td>AGN winds</td><td>4–4</td><td>1–4</td><td>2–4</td><td>3–4</td><td>2–4</td><td>4–3</td></tr><tr><td>Stellar feedback</td><td>3–3</td><td>1–1</td><td>1–1</td><td>1–1</td><td>1–1</td><td>2–1</td></tr></tbody></table>

**Notes.** For instance, the morphology of the *T*-field at $z \lesssim 1.5$ is mostly affected by jets, then by X-ray heating, stellar feedback, and AGN winds in this order.

anisms. It is the second most important process in shaping the global morphology of *T* and *n*<sub>HI</sub> fields at all studied redshift. The *n*<sub>gas</sub> and *n*<sub>H<sub>2</sub></sub> fields show quantitatively similar differences as for *n*<sub>HI</sub>-field, but given the relatively stronger impact of other feedback mechanisms, X-ray heating is ranked only as third and fourth. The morphology of *Z*-field is impacted the least by X-ray heating and this is limited to low redshift ($z \lesssim 1$).

**Jets.** The impact of AGN jets is the largest, and dominant over other feedback mechanisms, on the *T*-field, in particular at $z \lesssim 1.5–2$. This is somewhat expected as jets, when active, deposit large amounts of energy at substantial distances owing to their high velocities and hence large kinetic energies (Christiansen et al. 2020). Jets also impact the *Z*-field very significantly, again dominantly compared to other explored feedback mechanisms, by pushing hot metal-enriched outwards from massive halos (Borrow et al. 2020). As for the *T*-field, the effect of jets on the *Z*-field is strongest at $z \lesssim 1.5–2$. AGN jets impact the morphology of the *n*<sub>gas</sub>, *n*<sub>HI</sub>, and *n*<sub>H<sub>2</sub></sub> fields to a lesser extent compared to other feedback mechanisms, while they are competitive with AGN winds in shaping the *P*-field.

**AGN winds.** According to the implementation of the radiative mode of AGN feedback in SIMBA simulations, AGN winds have a negligible impact on the *T*-field. Their strongest effect is at $z \lesssim 1.5$ on the *P*-field morphology (comparable in strength to the effect of jets, for all the MFs but the volume filling fraction), and subdominant in the other fields, especially at high redshift. The reason is the typical speed of AGN winds, which in 14 Gyr in vacuum, pushes particles up to about 1 Mpc, namely, below the resolution limit set by effective smoothing.

**Stellar feedback.** Star formation-driven winds are found to play a crucial role in modifying the global morphology of the gas even at relatively large scales such as those probed here. While such winds have even lower velocities than the AGN winds, the cumulative amount of mass displaced is very large, substantially exceeding the mass in stars, and weighted towards high redshift since galaxies are small and thus mass loading factors are high. The stellar feedback is the dominant mechanism for all but the *T* and *Z* fields. Its impact on the *P*, *n*<sub>gas</sub>, *n*<sub>HI</sub>, and *n*<sub>H<sub>2</sub></sub> fields is highest at high redshift, progressively decreases, but stays important down to $z = 0$. HI field shows the strongest signature of the stellar feedback between redshift 1.5 and 0.5. As the star formation-driven winds are metal-loaded, it is not surprising to find that they impact the global morphology of the *Z*-field with the strength ranked right after jets although with reversed redshift dependence, namely, the signature of the stellar feedback is strongest at the highest redshift, while that of jets is almost absent at high redshift and becomes strongest at

$z = 0$. The *T*-field is impacted the least by the stellar feedback and is limited to a relatively narrow range of temperatures around $10^4$ K.

## 4.2. Time evolution of MFs

Minkowski curves provide a static picture of AGN and stellar feedback on large scales. The time evolution of the MFs of the gas fields for some specific excursion-set thresholds offers a complementary synthesis of the results. It highlights the possible role of stellar or AGN feedback mechanisms on top of hydrodynamics in driving morphological transitions, namely, abrupt or smooth changes of the MFs following immediately or with some delay the peak of stellar and AGN activity.

This analysis is illustrated in Figure 10, which focuses on three values of the six gas fields approximately spanning the full range covered by mean-field values $\bar{f}$ (indicated by symbols in Figures 4–9 and quoted below) for the Fiducial and NoAGN models. This figure compares the role of all AGN and stellar feedback mechanisms altogether (solid lines) against the role of stellar feedback only (dashed lines) for $z < 4$, namely in the last 11.8 Gyr. For reference, in SIMBA (Davé et al. 2019) the peak of star formation rate density occurs in the redshift range $1.5 < z < 2.5$ (dark shaded band), in agreement with observations (e.g. Madau & Dickinson 2014); the maximum of AGN luminosity function occurs approximately in the same redshift range for AGN with bolometric luminosity $10^{45} < L_{\text{bol}}/\text{erg s}^{-1} < 10^{46}$, and around $z \approx 1$ ($z \approx 2$; vertical lines) for AGN with $10^{44} < L_{\text{bol}}/\text{erg s}^{-1} < 10^{45}$ ($10^{46} < L_{\text{bol}}/\text{erg s}^{-1} < 10^{48}$), as measured by Habouzit et al. (2022) (see also Fontanot et al. 2020). Overall, the following can be deduced:

**Temperature.** AGN activity strongly modifies the morphology of the *T*-field at every redshift (i) by boosting the occupied volume $v_0$ at increasing times for $T > 10^4$ K and with a sudden contraction of colder regions between $z = 3$ and $z = 2$, at the same epoch when the number density of brightest quasars is largest (afterwards, AGN feedbacks and in particular jets become more ubiquitous and enriched by the faint population of quasars, starting to heat the gas on progressively larger and larger scales); (ii) by anticipating the shrinkage of the surface area $v_1$ at $z = 1.5–2$, especially for cold ($T = 10^4$ K) and warm ($T = 10^5$ K) isothermal surfaces, which are therefore globally rounder; (iii) by simultaneously increasing the convexity of warm and hot ($T = 10^6$ K) domains, which are globally less flat at respectively $z > 2$ and $z > 1$ (larger integrated mean curvature $v_2$) and forming a network dominated by tunnels afterwards (negative $v_3$); (iv) by changing the sign of the integrated mean curvature for the regions with $T > 10^4$ K, which is positive at all times when only stellar feedback is active whilst always negative (apart from the epoch $1.8 < z < 2.7$) because of AGN feedback, with the complementary domain, namely regions with $T \le 10^4$ K, forming a network of filaments (filled tunnels) at all probed times.

**Pressure.** The morphology of the *P*-field evolves over time in a non-trivial manner depending on the pressure value and feedback mechanisms. (i) The geometry of the isobaric domains with low pressure ($P = 10^3 k_B \text{ cm}^{-3}$ K) is only marginally affected after $z \approx 2$ by X-ray heating, AGN jets, and winds (solid lines), occupying slightly smaller volume and becoming moderately more irregular (larger surface density) than in presence of stellar feedback only, and their topology is almost unaltered as indicated by the Euler characteristic that remains vanishing at all epochs. (ii) When AGN feedback is included,
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Time evolution of MFs (top to bottom) for all the fields at three specific threshold values spanning the full range, for the Fiducial model (solid line) and for the NoAGN model (dashed line). Threshold values (blue, green, red in increasing order): $log(T[K]) = \{4, 5, 6\}$; $log(P/k_B[cm^{-3}K]) = \{3, 3.9, 4.8\}$; $log(n[cm^{-3}]) = \{-3, -2, -1\}$; $log(n_{HI}[cm^{-3}]) = \{-4, -3, -2\}$; $log(n_{H_2}[cm^{-3}]) = \{-4, -2.5, -1\}$; $log(Z[Z_\odot]) = \{-2, -1, 0\}$. The shaded band indicates the redshift range $1.5 < z < 2.5$ around the peak of star formation rate density, which approximately coincides with the maximum of AGN luminosity function for AGN with bolometric luminosity $10^{45} < L_{bol}/(erg s^{-1}) < 10^{46}$; vertical lines mark the maximum of the AGN luminosity function for AGN with $10^{44} < L_{bol}/(erg s^{-1}) < 10^{45}$ ($z \approx 1$) and $10^{46} < L_{bol}/(erg s^{-1}) < 10^{48}$ ($z \approx 2$). As discussed in Sect. 4.2, these curves display many distinct features depending on the captured feedback processes, supporting the informative value of the cosmic evolution of morphology (each feature likely encodes a specific signature of the underlying process).</figcaption></figure>

intermediate-pressure domains ($P = 10^{3.9} k_B cm^{-3} K$) instead shrink over time to a much larger extent, already since $z \approx 3$, while keeping the same surface until $z \approx 1$ as in presence of stellar feedback only. After $z \approx 0.5$, they start forming a network progressively dominated by isolated regions with positive integrated mean curvature and positive Euler characteristic, namely, with meatball topology at the present time. If AGN feedbacks were not active, intermediate-pressure domains would occupy 30% of the total volume, have a much more irregular shape, and form today a network still dominated by tunnels. (iii) The geometry and topology of high-pressure domains ($P > 10^{4.8} k_B cm^{-3} K$) are instead more strongly affected by AGNs, which firstly operate on small scales by cancelling out the expanding effect of stellar feedback; because of AGN, the high-pressure domains tend indeed to occupy a progressively smaller volume fraction with smaller surface density, namely, they are more clumpy, and form a network with meatball topology (positive $v_3$), especially when most quasars of intermediate brightness is in place, which coincides with the epoch around the peak of the star-formation activity.

**Total, H₁, and H₂ densities.** For the lowest-to-intermediate density thresholds of all the three density fields, the global morphology set by hydrodynamics and stellar feedback is unchanged by AGN feedbacks; the reason is that these densities are typical for very large scales, where gravitational interactions dom-

inate. AGN feedback modifies instead the morphology of high isodensity domains, where $n = 0.1 cm^{-3}$, $n_{HI} = 0.01 cm^{-3}$, and $n_{H_2} = 0.1 cm^{-3}$, with maximum efficiency at the epoch near the peak of star-formation activity and later, especially for the total and HI gas fields. Interestingly, while AGN feedback prompts the formation of a network with meatball topology at $1 < z < 2.5$, it does not modify the topology of the H₂ field in the same epoch, already established by hydrodynamics and stellar feedback while driving a spongy-like topology at $z < 1$ as deduced by the smaller value of $v_3$ indication of a larger number of tunnels. The sudden increase of $v_3$ of $n$ and $n_{H_2}$ (but surprisingly not $n_{HI}$) fields at $z < 0.5$ for intermediate and large values of density could be associated with a recent merging of filaments driven by AGN feedback; however, this claim demands a more careful study w.r.t. smoothing, which here is kept constant through cosmic time.

**Metallicity.** Metal-poor domains ($Z \le 0.01 Z_\odot$) maintain their morphology across time, with minor effect by AGN feedback, which makes them globally smoother (smaller $v_1$) and slightly flatter ($v_2$ closer to zero) especially for $z < 2$ compared to the Z-field in NoAGN model, with a minor change to the topology at $z \lesssim 1.5$. Conversely, the shape and topology of regions with higher metallicity, $Z = 0.1 Z_\odot$, and $Z = Z_\odot$, is strongly affected by AGN, especially for redshift respectively $z \lesssim 2.5$ and $z < 1.8$, mirroring the time evolution of the MFs of the H₂ field with
some delay. These domains progressively occupy larger volume fractions, have more irregular shapes as indicated by the larger $v_1$, and form a network dominated by filaments, respectively, at $z \le 2.5$ and $z \le 0.5$.

### 4.3. Phase diagrams and morphology: Joint analysis

SIMBA simulations allow for a field-based classification of baryonic environments based on both thermodynamical and chemical properties while splitting the role of AGN and stellar feedbacks. Bivariate distributions, hereafter also dubbed phase diagrams, can be used to infer average scaling relations between temperature, density, pressure, and metallicity of gaseous components, including the scatter accounting for degeneracies with hidden variables or some stochasticity. Estimating a quantitative relationship between the bivariate distributions (one-point spatial statistics) and the morphology of the excursion sets encapsulated in MFs (higher-order spatial statistics) is a topic of interest: do the phase diagrams capture the same information as the MFs? Viz., are regions featuring a thermodynamical and chemical state described by some scaling relation characterised by a well-defined typical morphology?

The smoothing adopted here for the computation of MFs only allows for a tentative joint analysis based on the more representative phase diagrams, shown in Appendix D at redshift $z = 0, 1, 2, 3, 5$ and only for the Fiducial and NoAGN models (Figures D.1 and D.2), intended to stress the effect of AGN feedback mechanisms on top of stellar feedback and hydrodynamics. By pin-pointing topologically motivated threshold values of the fields within the phase diagrams, one can attempt to find a relation between the thermo-chemical state of the gas and the percolation level and topology of its geometrical structure. The qualitative summary of results is the following (for more details, see aforementioned appendix):

- When AGN feedbacks are active, at $z \le 1$ warm regions with $10^5 < T/K < 10^7$ resulting after smoothing from undisturbed warm-hot intergalactic medium (WHIM) and warm circumgalactic medium (WCGM; Tumlinson et al. 2017) are almost dominated by tunnels, namely, they leave a network of cold filaments with $T \ge 10^6$ K. The T-field is further populated by colder cavities, roughly corresponding to collapsed haloes, which dominate the total volume at $1 \le z \le 2$ with temperature decreasing with time down to $10^3$ K at $z = 0$. When only stellar feedback is active, the volume is equally partitioned in warm and cold regions, the former dominated by bounded haloes and the latter by cavities, leaving a network of cold filaments characterised by a temperature two orders-of magnitude smaller than in presence of AGN feedbacks.

- The smoothed P-field at $z = 0$ is characterised by a single phase out of the actual three when all the feedback mechanisms are active, with an effective equation of state $P \propto n^{1.2}$ fitting the polytropic laws for isothermal and ultra-relativistic gases. Half of the volume is occupied by domains with pressure $1 \le P/k_B cm^{-3} \le 10^3$ and very noisy topology with vanishing Euler characteristic on average, namely, percolating the full volume. With AGN feedbacks turned off, while the one-point distribution of P values is essentially unchanged qualitatively, the topology of the smoothed P-field is reacher, admitting a network dominated by low-pressure cavities with $P \le 10 k_B cm^{-3}$ K. At higher redshift, especially for $1 < z < 2$, the smoothed P-field seems less dominated by cavities when AGN feedbacks are operational. A more suited smoothing procedure is mandatory to deduce

a solid conclusion about the morphology of the original P-field.

- The total, HI, and $H_2$ smoothed density fields are the most affected by the smoothing process, which destroys the bimodal or trimodal distributions leaving only the high-density regions. This lifting is less severe for the $n_{HI}$-field, at all redshifts and for the $n_{H_2}$-field at $z > 2$, which maintains its almost linear bis relation with the total gas n-field. The sparser regions surviving the smoothing occupy the largest fraction of the volume and are characterised by cavities. However, smoothing prevents the impact of feedback from morphology from being distinguished, which is otherwise quite clear instead in the original bivariate distributions.

- The joint analysis of the Z-field is perhaps the most intriguing. At all redshifts the one-point distribution of the smoothed and original field values are not very different and it does not distinguish between the Fiducial and NoAGN model; this is valid also for the $Z-n$ bivariate distributions. Instead, the MFs of excursion-sets suitably chosen to span the range of metallicities occupy very different regions of the smoothed $Z-n$ phase diagram. This is especially true at $z = 0$, where the joint analysis suggests that the denser medium is an intricate network dominated by metal-rich isolated regions with $Z \ge Z_\odot$ and filaments with only slightly smaller metallicity when also AGN feedback is present, or with filaments with $Z \approx 0.1 Z_\odot$ when only stellar feedback is active.

Altogether, MFs extend over the $T-n, P-n, n_{HI}-n, n_{H_2}-n$, and $Z-n$ phase diagrams tracing the mean values and scaling relations with small dispersion. The MFs analysis on smaller scales would allow us to capture the sparse domains removed by current smoothing (especially regions with small density and pressure, and very high temperature and metallicity) and relate the dispersion of phase diagrams to the fluctuations in MFs.

## 5. Conclusions

This study serves as an initial exploratory examination into the influence of AGN and stellar feedback mechanisms on the global morphology of the cosmic web. The primary objective is to comprehend how these small-scale processes shape the overall geometry and topology of thermo-chemical gaseous fields at circumgalactic and intergalactic scales across cosmic time. Investigating the morphology of gas fields, which constitute multiple facets of the cosmic web, is of significant scientific relevance. It allows us to discern the extent to which subgrid physics influences the geometry of the larger scales, irrespective of achievements at smaller scales, with the ultimate goal of constraining differently feedback physics<fnref n="2" />. To achieve this goal, we analysed the hydrodynamical simulation suite SIMBA, capturing the effect of feedback mechanisms including X-ray heating, AGN jets and winds, stellar, and hydrodynamical feedbacks on six fields, namely: the gas temperature, pressure, total density, HI density, $H_2$ density, and gas metallicity fields. A visual inspection of the diversity of 2D sections, as shown in Figure 1, offers an insight into the topic.

To quantify the morphology of the *T*, *P*, *n*, *n*<sub>HI</sub>, *n*<sub>H₂</sub>, and *Z* fields, we employed the MFs that are able to capture a more complex shape-based picture than that described by low-order statistics. Despite the absence of analytical models for the expectation

<footnote n="2">² The subgrid physics is indeed, in part, tailored to fit the internal and clustering properties of galaxies, aiming at preventing any over-cooling catastrophe.</footnote>
value of MFs of highly non-Gaussian fields, such as the SIMBA output, the numerical estimation of MFs still offer a powerful means to discriminate between the roles of AGN and stellar feedback.

Specifically, we performed an analysis of the MFs of the excursion sets of the ‘smoothed’ fields at different redshifts, spanning the range of $0 < z < 5$ and progressively accounting for all the SIMBA feedback mechanisms<fnref n="3" />. The results were given in Figures 4–9. These Minkowski curves constitute the starting point for all subsequent analyses. The uncertainty of the estimated Minkowski curves (not shown) deserves a comment. An intrinsic source of error is the sample variance, which is expected to not dominate within observed domains as long as the structures traced by the excursion set do not span a significant fraction of the box size. Moreover, the specific estimator of the MFs and the algorithm employed for their computation, including the discretisation and smoothing schemes, introduce an additional error. However, since we are interested in a global description of the morphology (and not in the minute, local details), we have neglected this second contribution and considered the Minkowski curves to be sufficiently precise and accurate to attain our objectives.

From the Minkowski curves we deduced a ranking of feedback mechanisms with respect to their impact on the global morphology of the excursion sets. Taking into account all redshifts, the following qualitative picture emerges (see Table 2): Stellar feedback shapes all but the *T* and *Z* fields, whose morphologies are instead mostly determined by AGN jets that are otherwise marginal for the other fields. Furthermore, AGN winds comprise the second feedback mechanism regulating the morphology of the *P*, *n*, and $n_{H_2}$ fields, while X-ray heating is for *T* and $n_{HI}$ fields, and stellar feedback is for the *Z*-field. The ranking could change if based on MFs for single values (thresholds) of the thermodynamical and chemical fields typical of specific systems, at some fixed redshift or in particular redshift intervals. However, such analysis goes beyond the scope of this study.

The analysis of the time-evolution of the MFs of excursion sets focusing on specific values of the thermodynamical and chemical fields allowed us to evaluate the strength of the feedback mechanisms as a function of cosmic time and to monitor whether they drive or alleviate morphological transitions (indicated by smooth trends or abrupt change of the MFs values; see the key Figure 10). The distinct features seen in that figure highlight the high information content encoded in the geometry and topology of the various physical fields beyond the total density commonly considered in past studies. Each feature likely encodes a specific process (as discussed in this work) that could eventually be used to infer the corresponding physical mechanism. This analysis highlights the potential of MFs as a probe for the cosmic evolution of the chemo-thermal cosmic web. For instance, the (smoothed) *T*-field at $10^5$ K ($10^6$ K) undergoes a topological transition around $z = 1.2$ ($z = 2$) only if AGN feedback is included; the *P*-field at $10^{4.8} k_B \text{ cm}^{-3}$ K and the *n*-field at $0.1 \text{ cm}^{-3}$ are instead protected from a temporary topological transition occurring around $z = 1.5$ and stopping at $z \approx 1$ if AGN feedback is active; the topology of the HI, $H_2$, and metallicity fields for $n_{HI} = 10^{-2} \text{ cm}^{-3}$, $n_{H_2} > 10^{-2.5} \text{ cm}^{-3}$, and $Z > 0.1 Z_\odot$ are also strongly altered by AGN feedback around $z = 1.5-2$.

However, the topological and morphological transition of these excursion sets continue until $z = 0$.

The analysis of the bivariate distributions or phase diagrams $T-n$, $P-n$, $n_{HI}-n$, $n_{H_2}-n$, and $Z-n$ jointly with MFs, shown in Figures D.1 and D.2, illustrates the capability of the latter to capture striking features of the diagrams (for the smoothed fields), for instance using maxima, minima, and vanishing points of the Minkowski curves that account for specific morphologies. Alternatively, and more efficiently, one can choose MFs that correspond to field values (thresholds) spanning the probability distribution functions of the smoothed field. This, in turn, suggests that the geometry of the fields indeed reflects the thermochemical phases of the IGM and, therefore, the micro physics of star formation and feedback as well. Indeed, the location of field values in phase diagrams characterised by average scaling relations (e.g., dense, high-pressure, isothermal regions) corresponding to thresholds associated with specific values of MFs (e.g., values of $v_0$ and $v_3$) suggests relationships between the typical morphology of spatial regions associated with that scaling (in the example, dense, high-pressure, isothermal regions cover a tiny fraction of the volume and are isolated, as indicated by the small value of $v_0$ and the large, positive value of $v_3$). Besides, the covariance matrix of a few values extracted from the MFs needed for data modelling is, in principle, much simpler to estimate than that of the full distribution functions.

The current analysis is confined to SIMBA suite of models, so it would be interesting to compare to other simulations such as ILLUSTRISTNG (Pillepich et al. 2018) or FLAMINGO (Schaye et al. 2023). For instance, it is known that SIMBA’s jet feedback causes wider baryon dispersal (Yang et al. 2022; Gebhardt et al. 2023) and more IGM heating (Christiansen et al. 2020) than other comparable simulations. Since jet feedback appears to drive the strongest features in the MFs, we would expect that comparisons to other large-volume cosmological simulations with different physical models would highlight areas of greatest discrimination between feedback prescriptions.

It would also be of interest to investigate the geometry of the baryonic cosmic web in simulations probing smaller scales such as NEW-HORIZON (Dubois et al. 2021), TEMPEST (Hummels et al. 2019), HESTIA (Libeskind et al. 2020; Damle et al. 2022), ROMULUS (Tremmel et al. 2017, 2019; Saeedzadeh et al. 2023), or Project GIBLE (Ramesh & Nelson 2023), and address the issue of time delays between small-scale trigger and large-scale response in more details. One should parse the range of possible subgrid physics while considering different sets of simulations, ideally spanning the full range of realistic values. Jointly with small-scale sub-galactic analysis, this could provide us with more leverage to disentangle the various processes captured by present-day feedback models. Eventually, variations of the key Figure 10 should be used to fit observed data: the next morphological study will assess the MFs of observable fields instead of elementary thermodynamical or chemical fields, quantifying the corresponding biases. For instance, we can consider the HI optical depth $\tau \propto n_{HI}$ probed by Lyman-$\alpha$ forest in high-resolution quasar spectra (e.g. McQuinn 2016; Japelj et al. 2019), the X-ray luminosity $L_X \propto T^3$ probed observing galaxy clusters (e.g. Zou et al. 2016; Bahar et al. 2022; Poon et al. 2023), the Compton parameter $Y \propto \int d\ell P$ proportional to the line-of-sight integral of the pressure for ideal electron gas, probed by SZ-surveys (Bleem et al. 2015; Kim et al. 2022). We could also consider sub-mm, infrared, or optical luminosities that probe the molecular gas as done within the MUSE Analysis of Gas around Galaxies (MAGG; Lofthouse et al. 2020), which map the environment of high-

<footnote n="3">³ Smoothing is a necessary step for the computation of MFs, which limits the original spatial resolution and biases gas fields. The smoothed distribution and the corresponding excursion sets of the total density are impacted, especially at $z < 3$, as traced by the lost bimodality of the one-dimensional probability distribution function. See Appendices C and D.</footnote>
redshift absorption line systems, or by probing specific metal
contents such as [O/H], [Mg/H] or [Fe/H] (Kewley et al. 2019).

The present shape-based multi-field study could also be com-
plemented by the analysis of the auto and cross-spectra or by the
environment-dependent wavelet power spectrum (Wang & He
2024) of all chemo-thermal fields, to highlight which scale
is most impacted by feedback and when; by studying multi-
persistent homology (Botnan & Lesnick 2022) and the sets
of critical events available in the initial phases (Cadiou et al.
2020), to assess the genuinely feedback-driven and gravity-
driven effects and disentangle the impact of the geometry of
the environment from internal processes. We also used comple-
mentary tools based on persistent homology (Sousbie 2011) and
eventually focusing on each set of critical points by stacking the
information in the frame set by eigenvectors of the curvature of
the fields (Kraljic et al. 2019). This offers a solution that would
jointly provide information on the local shape of the various gas
fields.

Acknowledgements. The authors thank D. Aubert, T. Buchert, G. De Lucia, C. Laigle, P. Monaco, and E. Nezri for their valuable comments and insights. This work has made use of the Horizon cluster on which the simulation was post-processed, hosted by the Institut d'Astrophysique de Paris; we thank S. Rouberol for running it smoothly. CS is grateful to INAF/Osservatorio Astronomico di Trieste for their hospitality. KK and CP thank the KITP for hosting the workshop 'CosmicWeb23: connecting Galaxies to Cosmology at High and Low Redshift'. This work is partially supported by the grant Segal ANR-19-CE31-0017 of the French Agence Nationale de la Recherche and by the National Science Foundation under Grant No. NSF PHY-1748958.

References

Adler, R. J. 1981, The Geometry of Random Fields, Classics in Applied Mathematics (Philadelphia: Society for Industrial and Applied Mathematics (SIAM))

Anglés-Alcázar, D., Davé, R., Faucher-Giguère, C.-A., Özel, F., & Hopkins, P. F. 2017, MNRAS, 464, 2840

Appleby, S., Davé, R., Kraljic, K., Anglés-Alcázar, D., & Narayanan, D. 2020, MNRAS, 494, 6053

Appleby, S., Park, C., Pranav, P., et al. 2022, ApJ, 928, 108

Appleby, S., Davé, R., Sorini, D., Lovell, C. C., & Lo, K. 2023, MNRAS, 525, 1167

Ayromlou, M., Nelson, D., & Pillepich, A. 2023, MNRAS, 524, 5391

Bahar, Y. E., Bulbul, E., Clerc, N., et al. 2022, A&A, 661, A7

Beisbart, C., Dahlke, R., Mecke, K., & Wagner, H. 2002, Lect. Notes Phys., 600, 238

Blake, C., James, J. B., & Poole, G. B. 2014, MNRAS, 437, 2488

Bleem, L. E., Stalder, B., de Haan, T., et al. 2015, ApJS, 216, 27

Bleem, L. E., Crawford, T. M., Ansarinejad, B., et al. 2022, ApJS, 258, 36

Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207

Bondi, H. 1952, MNRAS, 112, 195

Borrow, J., Anglés-Alcázar, D., & Davé, R. 2020, MNRAS, 491, 6102

Botnan, M. B., & Lesnick, M. 2022, ArXiv e-prints [arXiv:2203.14289]

Buchert, T., France, M. J., & Steiner, F. 2017, Class. Quant. Grav., 34, 094002

Cadiou, C., Pichon, C., Codis, S., et al. 2020, MNRAS, 496, 4787

Cen, R., & Ostriker, J. P. 1999, ApJ, 514, 1

Cen, R., & Ostriker, J. P. 2006, ApJ, 650, 560

Chabanier, S., Bournaud, F., Dubois, Y., et al. 2020, A&A, 643, L8

Chen, Z., Xu, Y., Wang, Y., & Chen, X. 2019, ApJ, 885, 23

Choi, E., Ostriker, J. P., Naab, T., & Johansson, P. H. 2012, ApJ, 754, 125

Choi, Y.-Y., Kim, J., Rossi, G., Kim, S. S., & Lee, J.-E. 2013, ApJS, 209, 19

Christiansen, J. F., Davé, R., Sorini, D., & Anglés-Alcázar, D. 2020, MNRAS, 499, 2617

Coles, P., Moscardini, L., Plionis, M., et al. 1993, MNRAS, 260, 572

Coles, P., Davies, A. G., & Pearson, R. C. 1996, MNRAS, 281, 1375

Colombi, S., Pogosyan, D., & Souradeep, T. 2000, Phys. Rev. Lett., 85, 5515

Damle, M., Sparre, M., Richter, P., et al. 2022, MNRAS, 512, 3717

Davé, R., Cen, R., Ostriker, J. P., et al. 2001, ApJ, 552, 473

Davé, R., Thompson, R., & Hopkins, P. F. 2016, MNRAS, 462, 3265

Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827

Donevski, D., Lapi, A., Małek, K., et al. 2020, A&A, 644, A144

Dubois, Y., Pichon, C., Welker, C., et al. 2014, MNRAS, 444, 1453

Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109

Ducout, A., Bouchet, F. R., Colombi, S., Pogosyan, D., & Prunet, S. 2013, MNRAS, 429, 2104

Einasto, M., Liivamägi, L. J., Tempel, E., et al. 2011, ApJ, 736, 51

Fabian, A. C. 2012, ARA&A, 50, 455

Fang, W., Li, B., & Zhao, G.-B. 2017, Phys. Rev. Lett., 118, 181301

Feldbrugge, J., van Engelen, M., van de Weygaert, R., Pranav, P., & Vegt, G. 2019, JCAP, 2019, 052

Fontanot, F., De Lucia, G., Hirschmann, M., et al. 2020, MNRAS, 496, 3943

Fukugita, M., Hogan, C. J., & Peebles, P. J. E. 1998, ApJ, 503, 518

Garratt, T. K., Coppin, K. E. K., Geach, J. E., et al. 2021, ApJ, 912, 62

Gay, C., Pichon, C., & Pogosyan, D. 2012, Phys. Rev. D, 85, 023011

Gebhardt, M., Anglés-Alcázar, D., Borrow, J., et al. 2023, ArXiv e-prints [arXiv:2307.11832]

Gilli, R., Comastri, A., & Hasinger, G. 2007, A&A, 463, 79

Gleser, L., Nusser, A., Ciardi, B., & Desjacques, V. 2006, MNRAS, 370, 1329

Gott, J. R., Melott, A. L., & Dickinson, M. 1986, ApJ, 306, 341

Gott, J. R., Park, C., Juszkiewicz, R., et al. 1990, ApJ, 352, 1

Haardt, F., & Madau, P. 2012, ApJ, 746, 125

Habouzit, M., Somerville, R. S., Li, Y., et al. 2022, MNRAS, 509, 3015

Hadwiger, H. 1957, Die Grundlehren der mathematischen Wissenschaften (Springer)

Hikage, C., & Matsubara, T. 2012, MNRAS, 425, 2187

Hikage, C., Schmalzing, J., Buchert, T., et al. 2003, PASJ, 55, 911

Hikage, C., Komatsu, E., & Matsubara, T. 2006, ApJ, 653, 11

Hopkins, P. F. 2015, MNRAS, 450, 53

Hopkins, P. F., & Quataert, E. 2011, MNRAS, 415, 1027

Hummels, C. B., Bryan, G. L., Smith, B. D., & Turk, M. J. 2013, MNRAS, 430, 1548

Hummels, C. B., Smith, B. D., Hopkins, P. F., et al. 2019, ApJ, 882, 156

Iwamoto, K., Brachwitz, F., Nomoto, K., et al. 1999, ApJS, 125, 439

James, J. B., Colless, M., Lewis, G. F., & Peacock, J. A. 2009, MNRAS, 394, 454

Japelj, J., Laigle, C., Puech, M., et al. 2019, A&A, 632, A94

Ji, S., Chan, T. K., Hummels, C. B., et al. 2020, MNRAS, 496, 4221

Kerscher, M., Schmalzing, J., & Buchert, T. 1996, ASP Conf. Ser., 94, 247

Kerscher, M., Schmalzing, J., Buchert, T., & Wagner, H. 1998, A&A, 333, 1

Kewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019, ARA&A, 57, 511

Kim, J., Golwala, S., Bartlett, J. G., et al. 2022, ApJ, 926, 179

Kraljic, K., Pichon, C., Dubois, Y., et al. 2019, MNRAS, 483, 3227

Kraljic, K., Laigle, C., Pichon, C., et al. 2022, MNRAS, 514, 1359

Kratochvil, J. M., Lim, E. A., Wang, S., et al. 2012, Phys. Rev. D, 85, 103513

Krumholz, M. R., & Gnedin, N. Y. 2011, ApJ, 729, 36

Lee, K.-G., Krolewski, A., White, M., et al. 2018, ApJS, 237, 31

Lee, J., Shin, J., Snaith, O. N., et al. 2021, ApJ, 908, 11

Li, Q., Narayanan, D., & Davé, R. 2019, MNRAS, 490, 1425

Libeskind, N. I., Carlesi, E., Grand, R. J. J., et al. 2020, MNRAS, 498, 2968

Liu, Y., Yu, Y., Yu, H.-R., & Zhang, P. 2020, Phys. Rev. D, 101, 063515

Lofthouse, E. K., Fumagalli, M., Fossati, M., et al. 2020, MNRAS, 491, 2057

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Mantz, H., Jacobs, K., & Mecke, K. 2008, J. Stat. Mech.: Theory Exp., 2008, 12105

Martínez, V. J., Starck, J.-L., Saar, E., et al. 2005, ApJ, 634, 744

Martizzi, D., Vogelsberger, M., Artale, M. C., et al. 2019, MNRAS, 486, 3766

Martizzi, D., Vogelsberger, M., Torrey, P., et al. 2020, MNRAS, 491, 5747

Matsubara, T. 2003, ApJ, 584, 1

Matsubara, T., & Yokoyama, J. 1996, ApJ, 463, 409

Matsubara, T., Hikage, C., & Kuriki, S. 2022, Phys. Rev. D, 105, 023527

McQuinn, M. 2016, ARA&A, 54, 313

Mecke, K., & Stoyan, D. 2002, Morphology of Condensed Matter (Berlin, Heidelberg: Springer)

Mecke, K. R., & Wagner, H. 1991, J. Stat. Phys., 64, 843

Mecke, K. R., Buchert, T., & Wagner, H. 1994, A&A, 288, 697

Melott, A. L., Cohen, A. P., Hamilton, A. J. S., Gott, J. R., & Weinberg, D. H. 1989, ApJ, 345, 618

Modest, H. I., Räth, C., Banday, A. J., et al. 2013, MNRAS, 428, 551

Munshi, D., Smidt, J., Cooray, A., et al. 2013, MNRAS, 434, 2830

Nakagami, T., Matsubara, T., Schmalzing, J., & Jing, Y. 2004, ArXiv e-prints [arXiv:astro-ph/0408428]

Nomoto, K., Tominaga, N., Umeda, H., Kobayashi, C., & Maeda, K. 2006, Nucl. Phys. A, 777, 424

Oppenheimer, B. D., & Davé, R. 2006, MNRAS, 373, 1265

Park, C., & Gott, J. R. 1991, ApJ, 378, 457

Park, C., & Kim, Y.-R. 2010, ApJ, 715, L185

Park, C., Gott, J. R., Melott, A. L., & Karachentsev, I. D. 1992, ApJ, 387, 1

Park, C., Gott, J. R., & Choi, Y. J. 2001, ApJ, 553, 33

Park, C., Pranav, P., Chingangbam, P., et al. 2013, J. Korean Astron. Soc., 46, 125
Peirani, S., Weinberg, D. H., Colombi, S., et al. 2014, ApJ, 784, 11

Perna, M., Lanzuisi, G., Brusa, M., Mignoli, M., & Cresci, G. 2017, A&A, 603, A99

Petri, A., Haiman, Z., Hui, L., May, M., & Kratochvil, J. M. 2013, Phys. Rev. D, 88, 123002

Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077

Planck Collaboration XIII. 2016, A&A, 594, A13

Pogosyan, D., Pichon, C., Gay, C., et al. 2009, MNRAS, 396, 635

Poon, H., Okabe, N., Fukazawa, Y., Akino, D., & Yang, C. 2023, MNRAS, 520, 6001

Powell, L. C., Slyz, A., & Devriendt, J. 2011, MNRAS, 414, 3671

Pranav, P., van de Weygaert, R., Vegter, G., et al. 2019, MNRAS, 485, 4167

Puetter, R., Gosnell, T., & Yahil, A. 2005, ARA&A, 43, 139

Rahmati, A., Pawlik, A. H., Raičević, M., & Schaye, J. 2013, MNRAS, 430, 2427

Ramesh, R., & Nelson, D. 2023, MNRAS, submitted [arXiv:2307.11143]

Saeedzadeh, V., Jung, S. L., Rennehan, D., et al. 2023, MNRAS, 525, 5677

Scannapieco, E., Pichon, C., Aracil, B., et al. 2006, MNRAS, 365, 615

Schaye, J., Crain, R. A., & Bower, R. 2015, MNRAS, 446, 521

Schaye, J., Kugel, R., Schaller, M., et al. 2023, MNRAS, 526, 4978

Schimd, C., & Sereno, M. 2021, MNRAS, 502, 3911

Schmalzing, J., & Buchert, T. 1997, ApJ, 482, L1

Schmalzing, J., & Gorski, K. M. 1998, MNRAS, 297, 355

Schmalzing, J., Buchert, T., Melott, A. L., et al. 1999, ApJ, 526, 568

Schmalzing, J., Takada, M., & Futamase, T. 2000, ApJ, 544, L83

Shandarin, S. F. 1983, Sov. Astron. Lett., 9, 104

Shandarin, S. F., & Yess, C. 1998, ApJ, 505, 12

Shandarin, S. F., & Zeldovich, I. B. 1983, Comments Astrophys., 10, 33

Smith, B. D., Bryan, G. L., Glover, S. C. O., et al. 2017, MNRAS, 466, 2217

Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51

Sorini, D., Davé, R., Cui, W., & Appleby, S. 2022, MNRAS, 516, 883

Sousbie, T. 2011, MNRAS, 414, 350

Spina, B., Porciani, C., & Schimd, C. 2021, MNRAS, 505, 3492

Springel, V. 2005, MNRAS, 364, 1105

Starck, J.-L., & Murtagh, F. 2006, Astronomical Image and Data Analysis (Berlin, Heidelberg: Springer-Verlag)

Sullivan, J. M., Wiegand, A., & Eisenstein, D. J. 2019, MNRAS, 485, 1708

Thomas, N., Davé, R., Anglés-Alcázar, D., & Jarvis, M. 2019, MNRAS, 487, 5764

Thomas, N., Davé, R., Jarvis, M. J., & Anglés-Alcázar, D. 2021, MNRAS, 503, 3492

Tomita, H. 1986, Progr. Theor. Phys., 76, 952

Tremmel, M., Karcher, M., Governato, F., et al. 2017, MNRAS, 470, 1121

Tremmel, M., Quinn, T. R., Ricarte, A., et al. 2019, MNRAS, 483, 3336

Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017, ARA&A, 55, 389

van de Weygaert, R., Pranav, P., Jones, B. J. T., et al. 2011, ArXiv e-prints [arXiv:1110.5528]

Vogeley, M. S., Park, C., Geller, M. J., Huchra, J. P., & Gott, J. R. 1994a, ApJ, 420, 525

Vogeley, M. S., Geller, M. J., Park, C., & Huchra, J. P. 1994b, AJ, 108, 745

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518

Wang, Y., & He, P. 2024, MNRAS, 528, 3797

White, S. D. M. 1979, MNRAS, 186, 145

Wiegand, A., Buchert, T., & Ostermann, M. 2014, MNRAS, 443, 241

Yang, T., Cai, Y.-C., Cui, W., et al. 2022, MNRAS, 516, 4084

Yess, C., Shandarin, S. F., & Fisher, K. B. 1997, ApJ, 474, 553

Yoshiura, S., Shimabukuro, H., Takahashi, K., & Matsubara, T. 2017, MNRAS, 465, 394

Zou, S., Maughan, B. J., Giles, P. A., et al. 2016, MNRAS, 463, 820

Zunckel, C., Gott, J. R., & Lunnan, R. 2011, MNRAS, 412, 1401
## Appendix A: MFs of Gaussian, weakly Gaussian, and log-normal random fields

A relatively small number of random fields admit analytical expressions for the expectation values of the MFs. The standard references are Gaussian random fields (GRF) $f$; in three dimensions, the MFs per unit volume of their excursion-set $C_v = \{x \in D | f(x) \ge v\sigma_0\}$ have analytical expectation value given by (Tomita 1986)

$$v_0^{\text{GRF}}(v) = \frac{1}{2} \text{erfc}(v/\sqrt{2}), \qquad (\text{A.1a})$$

$$v_1^{\text{GRF}}(v) = \frac{2}{3} \frac{\lambda}{2\pi} e^{-v^2/2}, \qquad (\text{A.1b})$$

$$v_2^{\text{GRF}}(v) = \frac{2}{3} \frac{\lambda^2}{(2\pi)^{3/2}} v e^{-v^2/2}, \qquad (\text{A.1c})$$

$$v_3^{\text{GRF}}(v) = \frac{\lambda^3}{(2\pi)^2} (v^2 - 1) e^{-v^2/2}. \qquad (\text{A.1d})$$

Here $\text{erfc}(x)$ is the complementary error function and the amplitude parameter $\lambda^2 = \sigma_1^2/3\sigma_0^2$ is the ratio between the variance of the gradient of the field, $\sigma_1^2 = \langle (\nabla f)^2 \rangle$, and variance of the field itself, $\sigma_0^2 = \langle f^2 \rangle_c$, which scales like the square of the inverse zero-crossing distance (Pogosyan et al. 2009).

The mean MFs densities of weakly non-Gaussian random fields incorporate exponentially damped polynomial corrections into the Gaussian expressions, with coefficient proportional to linear combinations of generalised skewness parameters

$$S_0 = \frac{\langle f^3 \rangle_c}{\sigma_0^4}, \quad S_1 = -\frac{3}{4} \frac{\langle f^2 \nabla^2 f \rangle_c}{\sigma_0^2 \sigma_1^2}, \quad S_2 = -\frac{9}{4} \frac{\langle (\nabla f \cdot \nabla f) \nabla^2 f \rangle_c}{\sigma_1^4}. \qquad (\text{A.2})$$

Up to second-order in perturbation theory, they read (Matsubara 2003; Gay et al. 2012; Matsubara et al. 2022)

$$v_0(v) = v_0^{\text{GRF}}(v) + \frac{\sigma_0}{(2\pi)^{1/2}} \frac{S_0}{2} H_2 e^{-v^2/2}, \qquad (\text{A.3a})$$

$$v_1(v) = v_1^{\text{GRF}}(v) + \frac{\sigma_0 \lambda}{9\pi} \left( \frac{1}{2} S_0 H_3 + S_1 H_1 \right) e^{-v^2/2}, \qquad (\text{A.3b})$$

$$v_2(v) = v_2^{\text{GRF}}(v) + \frac{2\sigma_0\lambda^2}{3(2\pi)^{3/2}} \left( S_0 H_4 + \frac{2}{3}S_1 H_5 + \frac{1}{3}S_2 \right) e^{-v^2/2}, \qquad (\text{A.3c})$$

$$v_3(v) = v_3^{\text{GRF}}(v) + \frac{\sigma_0 \lambda^3}{(2\pi)^2} \left( \frac{1}{6} S_0 H_5 + S_1 H_3 + S_2 H_1 \right) e^{-v^2/2}, \qquad (\text{A.3d})$$

where $v_\mu^{\text{GRF}}(\theta)$ are given by Equations (A.1) and $\{H_0, \dots, H_5\} = \{1, \theta, \theta^2 - 1, \theta(\theta^2 - 3), \theta^4 - 6\theta^2 + 3, \theta(\theta^4 - 10\theta^2 + 15)\}$ are Hermite polynomials. When the field $f$ is convolved with a Gaussian kernel, as routinely done in numerical computation

<figure><img src="image_12.png" /><figcaption>Fig. A.1. Examples of Minkowski curves of Gaussian, weakly non-Gaussian, and log-normal temperature random field (dotted, dashed, and solid lines, respectively) for three values of the mean temperature ($T = 10^4$ K, blue lines; $T = 10^5$ K, green; $T = 10^6$ K, red). The amplitude of the MFs of Gaussian and weakly non-Gaussian random fields are normalised to that of the LN random field; their apparent skewness is due to the logarithmic scale. For illustrative purpose, the parameters are fixed to $\sigma_0 = 0.3$, $(S_0, S_1, S_2) = (0.2, 1, 0.1)$, $\lambda = 25$.</figcaption></figure>

of fields sampled on a regular lattice, the skewness parameters become the function of the smoothing length $R$ and of the linear power spectral density of the field $f$; see Equations (13-16) in Nakagami et al. (2004).

For a log-normal (LN) field $f = F(g) \equiv \exp(g - \sigma_g^2) - 1$ related to a GRF $g$ with vanishing mean and variance $\sigma_g^2 = \ln(1 + \sigma^2)$, the excursion-sets $\{x \in D | g(x) \ge v\sigma_g\}$ and $\{x \in D | f(x) \ge F(v\sigma_g)\}$ are equivalent, namely, they share the same morphology. Accordingly, the average MFs per unit volume can be deduced from Equations (A.1) by substituting the threshold and amplitude parameters with the following expressions (Matsubara & Yokoyama 1996; Hikage et al. 2003):

$$v_{\text{LN}}(v) = \frac{\ln\left[\frac{(1 + v\sigma)\sqrt{1 + \sigma^2}}{\sqrt{\ln(1 + \sigma^2)}}\right]}{\sqrt{\ln(1 + \sigma^2)}}, \qquad (\text{A.4a})$$

$$\lambda_{\text{LN}}(\lambda) = \lambda \left[ \frac{\sigma}{1 + \sigma^2 \ln(1 + \sigma^2)} \right]^{1/2}. \qquad (\text{A.4b})$$

Figure A.1 shows some examples of a temperature random field with parameters fixed to fictitious values for illustrative purposes (see caption). For other non-Gaussian random fields monotonic functions of a GRF, such as chi-square, Student, Fisher, Rayleigh, or Maxwellian (see e.g. Adler 1981).
**Appendix B: 2D sections of smoothed gas fields**

Figures B.1-B.2 illustrate the large-scale imprint of the individual feedback processes operating on galactic scales on 2D

sections (thickness: $8h^{-1}$Mpc) of the gas fields at $z=0$ and $z=5$, when the smoothing by grid sampling and convolution by Gaussian kernel (needed for MFs computations) are taken into account.

<figure><img src="image_13.png" /><figcaption>**Fig. B.1.** Visualisation of the full box-size sections ($8\,\text{Mpc}/h$ thick and $50\,\text{Mpc}/h$ wide in both directions) of the 3D gas temperature, pressure, density (total, HI, H$_2$) and metallicity fields (from top to bottom) at the smoothing scale used for the computation of MFs and at $z=0$, for different SIMBA models progressively excluding feedback mechanisms (from left to right; see Table 1). Similarly to Figure 1, it illustrates how individual feedback processes operating on galactic scales leave different imprints on the gas fields at large scales.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. B.2. Same as Figure B.1, but at z = 5.</figcaption></figure>
**Appendix C: Bias of smoothed fields**

As stressed in the Guidelines (Sect. 3.1), the effective smooth-
ing of 1.56 h⁻¹Mpc resulting from the computation of the
MFs biases the probed regions of the one-point of the fields.
phase-space diagrams. This is illustrated for the T-field only
in Figure C.1, which compares its values at z = 0, 1, 2, 5 for
the Fiducial and NoAGN models before and after smoothing,
shown respectively in colours and contours (the fraction of out-
liers with |log T<sub>MF</sub> − log T<sub>SIMBA</sub>| > 1.5 is quoted in each panel).
At z ≤ 1 (first and second column panels) a Gaussian smooth-
ing kernel with rms-variance σ as large as 4 cells suppresses an
important fraction of cold regions with temperature T < 10⁴ K,
especially those with lowest total density, n < 10⁻⁴.⁵cm⁻³ (first
row panels; contamination for about 40-50 percent of the total),
which accordingly occupy sparse domains of size smaller than
~ 1.56h⁻¹Mpc at this epoch. Warm regions with T ≥ 10⁵ K
instead qualitatively maintain their original temperature also at
z = 0 where they occupy most of the volume, only marginally
polluted by colder and hotter regions. At z ≤ 2 the smoothed and
original T-field are quite similar. Therefore, although the MFs do
not probe exactly the original T-field, they still yield important
information about its global morphology for specific range of
underlying density and especially at z > 1. Smaller kernel width
of size σ = 1 cell would preserve the trace of cold domains,
although diminishing their number by a factor of about 2.

Depending on the feedback mechanism, the original *P*-field
in SIMBA simulations is composed by three dominant popula-
tions at low redshift and two at high redshift. Smoothing mixes
the intermediate and high pressure populations, suggesting that
they occupy sparse domains with characteristic size smaller than
1.56 *h*⁻¹Mpc. The difference is less important at *z* > 2, where the
smoothed *P*-field reasonably well represent the actual one.

The total, H1, and H2 smoothed density fields are the most
affected by smoothing, which destroy the bimodal or trimodal
distributions leaving only the high-density regions. This lifting is
less severe for the nHI-field, at all redshifts and for the nH2-field
at z > 2 which maintains its almost linear bis relation with the
total gas n-field. In particular, a similar suppression to the T-field
occurs for the nHI-field (not shown). The Gaussian smoothing
kernel with σ = 4 cells totally erases sparse cells with density
nHI < 10⁻⁷cm⁻³, which represent about 95 percent of the total
sample volume; domains with such a low value of density behave
as noise with a characteristic size smaller than ~ 1.56h⁻¹Mpc
at z = 0, embedded in higher-density domains. Simultane-
ously, this Gaussian kernel shrinks the high-density tail from the
original 10⁻⁷ < nHI/cm⁻³ < 1 range to 10⁻⁷ < nHI/cm⁻³ < 10⁻²,
qualitatively preserving its shape (skewness). Interestingly, con-
trary to the T-field a higher spatial resolution obtained with a
Gaussian kernel with σ = 1 cell is not sufficient to trace the
original nHI-field; it would totally deform the shape of the orig-
inal probability distribution function, artificially populating the
10⁻¹⁰ < nHI/cm⁻³ < 10⁻⁵ density range.

The one-point distribution of the Z-field is almost unchanged
by smoothing, making it a privileged and challenging field for
morphological analyses.

The computationally conservative choice of a Gaussian
smoothing kernel with rms-variance σ = 4 cells will there-
fore allow us to probe the global morphology of the WCGM
and, to a smaller extent, of the network of (cold and dense)
gaseous haloes. More sophisticated mass assignment schemes
and smoothing techniques such as adaptive kernels or wavelet
denoising (Puetter et al. 2005; Starck & Murtagh 2006) could be
more appropriate to maintain the morphology of the original
field, at the cost of a likely more ambiguous interpretation of
the probed scales (e.g. Martínez et al. 2005).
<figure><img src="image_15.png" /><figcaption>**Fig. C.1.** Scattering plot mapping the point-wise (cell) values of original SIMBA and smoothed temperature fields for the Fiducial and NoAGN models (colours and contours), as a function of redshift (left to right) and for three selections of the underlying original HI density field (top to bottom for low, mid, and high density; see legend). The fraction of outliers with $|\log T_{\text{MF}} - \log T_{\text{SIMBA}}| > 1.5$ is marked by dotted lines and indicated in each panel (in parenthesis, the fraction with respect to the full sample). Contours account for 100, 1 000, and 10 000 counts, not shown if sparse.</figcaption></figure>
# Appendix D: Phase diagrams of original and smoothed fields

We describe here the analysis of the principal bivariate distributions or phase diagrams of the *T*, *P*, *n*, *n*<sub>HI</sub>, *n*<sub>H₂</sub>, and *Z* fields for the Fiducial and NoAGN models at redshift *z* = 0, 1, 2, 3, 5, before and after the smoothing required by the computation of MFs. This supports the tentative joint analysis of (smoothed) phase diagrams and morphology. Figures D.1-D.2 report the logarithmic counts computed on the 128³ grid points both from the output of raw SIMBA fields without additional smoothing (contour lines) and after smoothing (colours). The over-plotted symbols pinpoint field values *f*<sub>*</sub> defining excursion-sets with noteworthy MFs: excursion-sets filling half the volume, *v*<sub>0</sub>(*f*<sub>*</sub>) = 0.5 (star); excursion-sets with vanishing integrated mean curvature, *v*<sub>2</sub>(*f*<sub>*</sub>) = 0 (diamond); excursion-sets dominated by isolated regions, namely, with maximum *v*<sub>3</sub>(*f*<sub>*</sub>) (triangle); excursion-sets dominated by tunnels, namely, with minimum *v*<sub>3</sub>(*f*<sub>*</sub>) (dot circle); and excursion-sets with cavities, namely, where *v*<sub>3</sub>(*f*<sub>*</sub>) has the second highest maximum (filled circle). We note that the *v*<sub>3</sub>(*P*) curve at *z* = 0 is very noisy for low-pressure values and vanishing on average, showing only a global maximum and no minimum and no secondary maximum; the corresponding symbols are absent.

A detailed description of the phase diagrams follows:

* *T* − *n*: The temperature-density parameter space is commonly adopted to define the baryonic cosmic-web environments (Fukugita et al. 1998; Cen & Ostriker 1999; Davé et al. 2001; Cen & Ostriker 2006; Powell et al. 2011). Referring to the *z* = 0 phase diagram of the Fiducial model (Figure D.1, top-left panel), one can identify as hot medium the regions with *T* > 10⁷ K, typically under-dense or diffuse; warm-hot intergalactic medium (WHIM) and warm circumgalactic medium (WCGM; Tumlinson et al. 2017) are regions with 10⁵ < *T*/*K* < 10⁷ and, respectively, *n* ≲ 10⁻⁴ cm⁻³ and 10⁻⁴ ≲ *n*/cm⁻³ < 0.13; diffuse intergalactic medium and collapsed haloes are regions with *T* < 10⁵ K and same bounds on density as WHIM and WCGM; and star forming regions with *T* < 10⁷ K and *n* > 0.13 cm⁻³. These bounds are consistent with Martizzi et al. (2019, 2020), who performed a similar classification based on ILLUSTRISTNG simulations in the *T*−*n*<sub>HI</sub> parameter space (on average *n*<sub>HI</sub> ≈ *n* for *n* > 10⁻⁵ cm⁻³; see later). In agreement with the analysis of bias of smoothed field reported in the Appendix D (see Figure C.1), smoothing erases the lower and higher-density regions (*n* < 10⁻⁶ cm⁻³ and *n* > 10⁻² cm⁻³ at *z* = 0) and the cold and hot regions (*T* < 10².⁸ K and *T* > 10⁶.⁸ K), so that MFs finally probe the morphology of excursion-sets that barely corresponds to the WCGM approximately at all redshift. The conclusion is qualitatively similar when X-ray heating is not included (not shown).

If only stellar feedback is taken into account (see Figure D.2, first line, contour lines), the hot medium and the WHIM are less prominent, and the WCGM has an average lower temperature, especially at $z \le 1$; the phase-space diagrams of the Fiducial and NoAGN models are similar only at high redshift (see panel at $z = 5$) when the AGN population is not sufficiently energetically powerful to affect the one-point statistics on large scales. The phase-space after smoothing for MFs (shaded areas in the figure) biases the raw bivariate distribution to a smaller extent in favour of a broader range of temperatures and densities than in the Fiducial model, here including colder and less dense regions especially at lower redshift, while it is qualitatively similar at $z \ge 1$. Interestingly, without AGN feedback, the MFs defined by threshold values (*n*<sub>*</sub>, *T*<sub>*</sub>) marked by sym-

bols cover a wider portion of the phase-space, pinpointing also the collapsed haloes. This correspondence suggests that when investigating the feedback mechanisms operating on gas fields, a small number of MFs computed for appropriate excursion sets of *T* and *n*-fields (namely, a few simple numbers, with associated small-size covariance matrix) provide complementary, valuable, and perhaps simpler information than one-point (bivariate) distributions (namely, continuous functions, with large non-trivial covariance matrices).

* *P* − *n*<sub>HI</sub>: At all redshifts, and regardless of the feedback mechanism, the pressure-density diagram shown in Figure D.1 (second row, contours) indicates three phases: high-pressure domains in high-density regions (*P*/*k*<sub>B</sub> > 100cm⁻³K and *n* > 10⁻⁴ cm⁻³ or *n*<sub>HI</sub> > 10⁻⁵ cm⁻³, at *z* = 0), in which the gas is almost isothermal, namely, *p* ∝ *n*; intermediate-pressure domains in low-density regions (0.1 < *P*/*k*<sub>B</sub>cm⁻³K < 100 and *n* ≈ 10⁻⁷ cm⁻³ or *n*<sub>HI</sub> ≈ 10⁻¹³ cm⁻³ at *z* = 0), in which the gas has on average *P* ∝ *n*<sup>Γ</sup>, with a polytropic coefficient for diatomic molecules, Γ = 7/5; and low-pressure domains in low-density regions (*P*/*k*<sub>B</sub> < 0.01cm⁻³K and *n* ≈ 10⁻⁷ cm⁻³ at *z* = 0), in which the gas behaves on average like monoatomic molecules with Γ = 5/3. The smoothing for MFs computation (shaded in the figure) removes a large part of the second phase and completely erases the third one, allowing for a morphological description of the high-pressure domains only. However, this bias is less crucial for the smoothed *P*-field at higher redshift, almost absent at *z* = 5. Very similar conclusions can be deduced for the NoAGN model.

* *n*<sub>HI</sub> − *n*: This diagram shows the relative bias of the neutral hydrogen with respect to the total gas density, which is almost monotonic; see the third rows in Figures D.1-D.2. Three phases coexist when all the AGN feedbacks are active; the one with lower density (*n*<sub>HI</sub> ≈ 10⁻¹⁵ cm⁻³ at *z* = 0) is a genuine product of AGN activity, especially the X-ray heating and AGN jets, absent only when stellar feedback is considered. Like in the pressure-density diagram, regardless of the feedback mechanism considered, the highest HI density domains survive to smoothing for MFs and without any significant bias; MFs can therefore accurately sample this phase, which on average is linearly biased to the underlying total gas density field (dashed line in the diagrams).

* *n*<sub>H₂</sub> − *n*: Analogously to the previous phase diagram, this one yields the relative bias of molecular hydrogen with respect to the total gas density; see the fourth rows in Figures D.1-D.2. As expected, proportionality holds only at high redshift and high density. AGN feedback mechanisms, especially X-ray heating and jets, increase the low-density phase, producing a third phase at *z* < 2 with mean density *n*<sub>H₂</sub> ≈ 10⁻⁹ cm⁻³ in the densest regions with 10⁻⁴ cm⁻³ < *n* < 10⁻² cm⁻³. Smoothing only drops the phase at a higher density, biasing the original population at *z* = 0 but not at higher redshift. Consistently, MFs accurately probe this H₂ phase without important bias in the redshift range 1 < *z* < 5.

* *Z* − *n*: The metallicity-density relation is shown in bottom panels of Figures D.1-D.2. Before smoothing, hydrodynamical interactions possibly augmented by stellar feedback produce two phases occupying regions with average density *n* ≈ 10⁻⁷ cm⁻³ and 10⁻² cm⁻³. A third phase with high metallicity (*Z* > *Z*<sub>☉</sub>) is produced in the low-density regions (10⁻⁸ cm⁻³ < *n* < 10⁻⁶ cm⁻³) since *z* = 2 by AGN feedback, specifically X-ray heating and jets. At *z* = 0 this last population becomes one of the two dominant ones, with average metallicity 0.1 < *Z*/*Z*<sub>☉</sub> < 10, at the expense of domains with lower values that survive only in the absence of AGN feedback (the relationship between the
<figure><img src="image_16.png" /></figure>

**Fig. D.1.** Phase diagrams for the Fiducial model at redshift $z = 0, 1, 2, 5$ for the original and smoothed fields (contours and density plot, respectively; lines and colours account for log-counts over five decades between $10$ and $10^5$). As a reference, some polytropic laws are superposed to the $P-n$ diagram at $z = 0$ (see legend) and the linear bias relation is indicated on the $n_{HI} - n$ bivariate distribution at $z = 0$. Opaque symbols pinpoint threshold values where $v_0 = 0.5$ (star), $v_2 = 0$ (diamond), $v_3$ is maximum (triangle), minimum (dot circle), and has a second global maximum (filled circle), illustrating the relationship between thermodynamical and chemical states, possibly described by some scaling relation, and the average morphology of the corresponding excursion-set.

$H_2$ and the Z phases pumped by AGN feedback at $z < 1$ will be considered elsewhere). The distribution of the smoothed Z-field is almost unchanged, with phase mixing derived from smoothing

of the *n*-field. Also in this case the MFs of well-chosen Z-fields pinpoint the high-density phase while fairly well reproducing the range of Z values, at least until $z \approx 3$.
<figure><img src="image_17.png" /><figcaption>**Fig. D.2.** Same as Figure D.1 but for the NoAGN model. As expected, galactic nuclei impact jointly both the phase diagrams and their geometric and topological counterparts.</figcaption></figure>