© The Authors 2024

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; 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-α 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₂) in a Jeans mass-resolving pressurised interstellar medium (ISM; with the hydrogen density n_H ≥ 0.13 cm⁻³), following Davé et al. (2016). The computation of the H₂ 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 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⁵ 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_⊙ 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 cm⁻³) 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⁻¹ Mpc comoving volumes with 512³ gas elements and 512³ dark matter particles with the mass resolution of 1.82 × 10⁷ M_⊙ for gas and 9.6 × 10⁷ M_⊙ for dark matter particles. Cosmology adopted in SIMBA is a standard ΛCDM compatible with results from Planck Collaboration XIII (2016), namely Ω_m = 0.3, Ω_Λ = 0.7, Ω_b = 0.048, H₀ = 68 km s⁻¹ Mpc⁻¹, σ₈ = 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⁻¹ Mpc thick at z = 0.

Fig. 1. Visualisation of full box-size sections (8 h−1 Mpc thick and 50 h−1 Mpc wide in both directions) of the unsmoothed 3D temperature, pressure, total density, HI density, H2 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, nH2, 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.

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, 𝒞_θ = {x ∈ D | f(x)≥θ}, where f denotes a generic field defined over a spatial domain D and θ the threshold (see Figures 2 and 3 for illustration). In three dimensions, the MFs of 𝒞_θ 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, respectively, of the body integrated over its surface. Owing to the Gauss-Bonnet theorem, the Gaussian curvature is proportional to the Euler characteristic χ = G/4π and linearly related to the genus g = 1 − χ/2 of 𝒞_θ, 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 the thresholded field f(x), accounting for the topology of 𝒞_θ. 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).

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 × 103 K (top middle), T ≥ 8 × 103 K (top right), T ≥ 1 × 104 K (bottom left), T ≥ 3 × 104 K (bottom middle) and T ≥ 5 × 104 K (bottom right). For illustrative purposes, 3D maps correspond to fields before applying additional Gaussian smoothing of 1.56 h−1 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.

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 ≥ 7 × 102 K (top middle), T ≥ 104 K (top right), T ≥ 105 K (bottom left), T ≥ 106 K (bottom middle) and T ≥ 107 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.

Here, we adopt the normal based on the so-called intrinsic volumes (V₀, V₁, V₂, V₃) = (V, A/6, H/3π, χ) 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:

$$\begin{array}{cc}\hfill v_0(\theta )& =\frac{1}{\mathcal{V}}{\displaystyle {\int }_{{\mathcal{C}}_{\theta }}}\mathrm{d}V(\mathit{x}),\hfill \end{array}$$(1a)

$$\begin{array}{cc}\hfill v_1(\theta )& =\frac{1}{6}\frac{1}{\mathcal{V}}{\displaystyle {\int }_{\partial {\mathcal{C}}_{\theta }}}\mathrm{d}A(\mathit{x}),\hfill \end{array}$$(1b)

$$\begin{array}{cc}\hfill v_2(\theta )& =\frac{1}{6\pi }\frac{1}{\mathcal{V}}{\int }_{\partial {\mathcal{C}}_{\theta }}\mathrm{d}A(\mathit{x})(\frac{1}{R_1(\mathit{x})}+\frac{1}{R_2(\mathit{x})}),\hfill \end{array}$$(1c)

$$\begin{array}{cc}\hfill v_3(\theta )& =\frac{1}{4\pi }\frac{1}{\mathcal{V}}{\displaystyle {\int }_{\partial {\mathcal{C}}_{\theta }}}\mathrm{d}A(\mathit{x})\frac{1}{R_1(\mathit{x})R_2(\mathit{x})},\hfill \end{array}$$(1d)

where dV(x) and dA(x) are the elementary volume and area of the excursion sets of the six random fields f = T, P, n_gas, n_HI, n_H2, Z (hereafter referred as f-fields), which have boundary surface ∂𝒞_θ with principal local curvature radii R_1, 2(x). Also, 𝒱 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).

3. Morphological analysis

3.1. Guidelines

Figures 4–9 which are discussed in the following subsections show the four MFs per unit volume, v_μ (μ = 0, 1, 2, 3), of the excursion-set of the three-dimensional temperature, pressure, density (total, H I, 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 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 Δ in the bottom sub-panels) as follows:

• impact of X-ray heating: v_μ[Fiducial]−v_μ[NoX],

• impact of jets: v_μ[NoX]−v_μ[NoJets],

• impact of AGN winds: v_μ[NoJets]−v_μ[NoAGN],

• impact of stellar feedback: v_μ[NoAGN]−v_μ[NoFeedback].

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, Δ. 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 ≲ 2), but also by the X-ray heating at both high (z ∼ 4 − 5) and low redshift (z ≲ 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 ∼ 2 − 3.

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.

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.

Fig. 7. HI density. Minkowski curves for the neutral atomic hydrogen number density, or nHI-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.

Fig. 8. H2 density. Minkowski curves for the molecular hydrogen number density or nH2-field, analogous to Figure 4. Its evolution is very similar to that of the total gas density; see Sect. 3.4.

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 ≳ 0.01 Z⊙ for z < 1.5, while at earlier time the morphology is maximally influenced by stellar feedback, especially in regions with Z ≲ 0.1 Z⊙; see Sect. 3.5.

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, $\overline{f}(\mathit{x})$; we will refer to these domains as $\overline{f}$-regions.

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

• The volume filling fraction v₀ is a monotonically decreasing function of the threshold θ, tending to zero if the highest values of the field f are concentrated in point-like regions. Likewise, 1 − v₀ measures the volume fraction of the complementary domain D \ 𝒞_θ = {x ∈ D | f(x) < θ}. We note that for the fixed value of the threshold, v₀ increases (decreases) with time when the corresponding domains expand (contract). When approximated by a step function H(θ_* − θ), like the temperature field at high redshift, v₀(θ) accounts for a homogeneous (single-phase) field with characteristic value f(x) = θ_*. Correspondingly, the density area v₁ is close to a Dirac-delta, indeed v₁(θ)≈dv₀(θ)/dθ. Consistently, a multi-phase field with sharp domains is accounted for by a piecewise constant function, namely, v₀(θ) = ∑_iv_0, iH(θ_i − θ) with ∑_iv_0, i = 1, the ith 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₁(θ) have therefore on average a more regular surface than domains with maximum v₁. Note also that the apparent skewness of v₁(θ) is due to the logarithmic scale of the threshold and typically grows with time.

• The integrated-mean-curvature density v₂(θ) is positive (negative) for domains C_θ that are on average convex (concave), which are dominated by lumps (cavities). The opposite is true for the complementary domains D \ 𝒞_θ.

• 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₃(θ) curve has a characteristic M-shape: it is positive for small values of the threshold θ, which are typically smaller than the mean value $\overline{\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 θ, with 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 >\overline{\theta }$ or $\theta <\overline{\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 rightmost zero of v₃(θ) 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₁(θ) curve peaked at $\theta =\overline{\theta }$, odd v₂(θ) with equal absolute amplitude of extrema and vanishing for $\theta =\overline{\theta }$, and even v₃(θ) with a negative minimum at $\theta =\overline{\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 ∇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₃(θ), especially at low redshift as a result of non-linear gravitational evolution and shot noise (in observed fields, additional deviations are due to redshift space distortions and, if any, primordial non-Gaussianities). Apart from some specific deviations and especially at high redshift, z > 4, their profile resembles that of log-normal random fields.

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

• Isolated regions with field values larger than in the surrounding environment, f > f_env (e.g. hot bubbles in cold or warm IGM), have small v₀, large (small) v₁ if they have smooth (wrinkly) boundary surface, v₂ > 0, and large v₃;

• Isolated bubbles with f < f_env (e.g. cold bubbles in warm IGM) have large values for v₀, surface density v₁ like for isolated regions, while v₂ < 0, and large values for v₃ like for isolated regions;

• Network of filaments with f > f_env (e.g. hot, dense, or metal-rich filaments alimented by local stellar activity) have large (small) v₀ if on average filaments are thick (thin), v₁ as above, v₂ > 0 and smaller for cylindrical shape (the largest radius of curvature diverges), and v₃ ≈ 0 if the network is percolating;

• Network of filaments with f < f_env, namely, an excursion set 𝒞_θ 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₀ if filaments are on average thin (thick), v₁ as before, v₂ < 0 and larger for cylindrical shape, and v₃ < 0 and large in 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⁻¹ 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 v₀(T) is a step-like function with the transition occurring around warm regions with T ≃ 10⁴ 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 ∼10³ 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⁴ 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₁(T) is symmetric and very peaked around 10⁴ 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₁(T) becomes and remains right-tailed until z ∼ 2 and left-tailed afterwards, in favour of regions with temperature smaller than the mean value at late times, namely, for $T<\overline{T}\simeq {10}^6$ K at z = 0 (marked by a circle). These regions are concave (v₂ < 0) and dominated by tunnels (v₃ < 0), especially at redshift z = 1 for T ≳ 10⁴ K and at z = 0.5 for T ≳ 10⁵ K, suggesting a network of cold filaments in a warm percolating IGM. Interestingly, at z = 0 concave regions, namely tunnels and bubbles, with T ≲ 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⁴ K as very special: besides marginal deviations at z < 0.5, at all times they fill about one-third of volume (v₀ ≃ 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 Δ. 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⁴ K, skew the v₁(T) curve by decreasing the surface area of cold (T ≃ 10⁴ K) regions by ∼7% and increasing that of hot regions (T ≳ 10⁶ K) by almost a factor of 10, and almost reverse the convexity of T > 10⁴ 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 ≲ 10⁴ K, without jets they would occupy also cooler regions (the first zero of v₃(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₃(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₃, or tunnels, increases with time, reaching the minimum value at z = 0 for T = 10⁶ 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⁴ K ≲ T ≲ 10^4.2 K) and a much broader range at z = 0 (10^3.3 K ≲ T ≲ 10⁵ K); after z = 2, because of jets, the percolation threshold increases dramatically, up to two orders of magnitude at z = 0 from T ≃ 10^4.5 K to T ≃ 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 ≃ 10⁴ K) and hot (T ≳ 10⁵ K), stellar feedback acts at intermediate redshift (1 < z < 2.5) and on warm regions (10^3.5 K ≲ T ≲ 10⁵ 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⁴ 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₃ and therefore the number of hot clumps (T > 10⁵ K) at z = 0.5, and enhance the peculiar topology of regions with 10^3.7 K ≲ T ≲ 10^4.4 K already set by jets, by increasing (decreasing) the number of tunnels in regions with temperature smaller (larger) than T = 10⁴ 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 $\overline{T}$-regions¹. Focusing on the volume filling fraction, the $\overline{T}$-regions occupy less than half the volume when all the AGN and stellar feedback mechanisms are included, with $v_0(\overline{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⁴ 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 $\overline{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(\overline{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₀ of T ≃ 10⁴ K regions by about 3% and only at z = 0.

In the full model, the density of surface area of $\overline{T}$-regions, $v_1(\overline{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(\overline{T})$ and $v_3(\overline{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 $\overline{T}$-regions with $\overline{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(\overline{T})\simeq 0$); namely, the $\overline{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<\overline{T}$ dominated by tunnels. This corresponds to a network of cold filaments surrounded by hotter 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<\overline{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₀ and v₁ curves at P > 10² k_B cm⁻³ 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₃ for the highest pressure regions (P > 10⁵ k_B cm⁻³ K). However, during this epoch the $\overline{P}$-regions tend to occupy only a small fraction of the total volume ($v_0(\overline{P})\simeq 0.2$ at z = 5 and < 0.05 at z = 4), suggesting that apart from sparse, small, and convex (spherical) high-pressure clumps most of the volume is low pressure with $P<\overline{P}$. For z < 4 and until z ∼ 1.5 the fraction of the total volume occupied by high-pressure domains progressively decreases, the $\overline{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 $\overline{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 cm⁻³ K disappear; regions with $P<\overline{P}={10}^{3.4}{k}_{\mathrm{B}}{\mathrm{cm}}^{-3}$ K occupy more than 60% of the volume and are essentially concave with v₃ ≃ 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 cm⁻³ 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 ≳ 10^4.9 k_B cm⁻³ K to occupy less than 15−20% of the total volume, quenching the time evolution of their shape and topology until z ∼ 2. This is visible by comparing the MFs of the $\overline{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 ≳ 10^2.9 k_B cm⁻³ K, the latter at z < 1.5 (not considering the highest pressure regions with P/k_B ≳ 10^5.9 cm⁻³ 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₀ by 20%, v₁ by more than 40%, v₂ by 50%, and v₃ by up to 100%. At higher redshifts (z ≳ 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 (HI) 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 $\overline{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 $\overline{n}$-regions, or equivalently a network of filaments with $n<\overline{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₃, marked by circles is almost constant in this period); after z = 1 it grows again, at the expenses of higher density clumps (n ≃ 10⁻² cm⁻³) 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₂ of $\overline{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 $\overline{n}$ and, correspondingly, the morphological evolution of $\overline{n}$-domains are instead quenched until z ≃ 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₂). 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 $\overline{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⁻³ ≲ n/cm⁻³ ≲ 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₀ of both overdense and underdense regions at z > 2, increase (decrease) v₀ of overdense (underdense) regions in the range 0.5 < z < 2, and decrease v₀ of both overdense and underdense regions at z < 0.5. The effect on v₁ and v₂ 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 ≤ 10⁻⁵ cm⁻³ (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₀ ≈ 1), very small surface density, and positive integrated mean curvature (correspondingly, the complementary excursion sets with n > 10⁻⁵ cm⁻³ have v₁ ≳ 0 and v₂ ≲ 0), and very small and positive Euler characteristic density, v₃ ≳ 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 cm⁻³ at z = 0 (n ≳ 2 cm⁻³ at z = 5).

3.4.2. H I 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⁻² cm⁻³ < n_HI < 10^−0.5 cm⁻³) 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 cm⁻³ < n_HI < 10^−1.5 cm⁻³) become more spherical (smaller v₁, larger v₂) and more abundant (larger v₃); (iii) an abrupt change of the global morphology occurs at z = 0.5, the volume filling fraction v₀ 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({\overline{n}}_{\mathrm{HI}})\simeq 0.5$) and is on average flat ($v_2({\overline{n}}_{\mathrm{HI}})\simeq 0.5$), namely, it is much more Gaussian although the maxima of v₃ 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₀ by up to 30%, especially for intermediate density regions; it increases (decreases) the surface density v₁ of regions with higher (lower) density than the critical value n_HI ≃ 10^−2.5 cm⁻³ by up to 70%; it does change the density of integrated mean curvature v₂, again with respect to the critical value n_HI ≃ 10^−2.5 cm⁻³; and increases the number of tunnels of excursions sets for a narrow range of thresholds, 10^−2.5 < n_HI/cm⁻³ ≃ 10⁻².

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 cm⁻³ 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_H2-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 ≃ Z_⊙; and (iii) a collection of small bubbles with metallicity increasing with time from 10^−2.5 Z_⊙ at z = 5 to 10^0.5 Z_⊙ 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 $\overline{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_⊙ ≲ 0.018) decreases with time, while for regions with intermediate metallicity (0.01 < Z/Z_⊙ ≲ 1) it attains a maximum and then decreases; at late time (z = 0) both attain negative v₂ and form a network dominated by tunnels. Metal-rich regions with Z/Z_⊙ ≳ 1.8 are instead characterised by v₁ increasing with time, and positive v₂ and v₃, 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_⊙ ≲ 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_⊙ > 0.01, and especially those with Z/Z_⊙ ≃ 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_⊙ < 0.89 (Z/Z_⊙ > 0.89) by 80% at z = 0, changing their mean curvature and topology, enhancing that of metal-rich regions with Z/Z_⊙ > 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_⊙ ≃ 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₀, as Figure 1 would wrongly suggest) at z ≲ 1.5 and z ≳ 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 mechanisms. It is the second most important process in shaping the global morphology of T and n_HI fields at all studied redshift. The n_gas and n_H2 fields show quantitatively similar differences as for n_HI-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 ≲ 1).

Jets. The impact of AGN jets is the largest, and dominant over other feedback mechanisms, on the T-field, in particular at z ≲ 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 ≲ 1.5 − 2. AGN jets impact the morphology of the n_gas, n_HI, and n_H2 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 ≲ 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_gas, n_HI, and n_H2 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⁴ 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 $\overline{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⁴⁵ < L_bol/(erg s⁻¹) < 10⁴⁶, and around z ≃ 1 (z ≃ 2; vertical lines) for AGN with 10⁴⁴ < L_bol/(erg s⁻¹) < 10⁴⁵ (10⁴⁶ < L_bol/(erg s⁻¹) < 10⁴⁸), as measured by Habouzit et al. (2022) (see also Fontanot et al. 2020). Overall, the following can be deduced:

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/kB[cm−3K]) = {3, 3.9, 4.8}; log(n[cm−3]) = { − 3, −2, −1}; log(nHI[cm−3]) = {−4,−3,−2}; log(nH2[cm−3]) = { − 4, −2.5, −1}; log(Z[Z⊙]) = { − 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 1045 < Lbol/(erg s−1) < 1046; vertical lines mark the maximum of the AGN luminosity function for AGN with 1044 < Lbol/(erg s−1) < 1045 (z ≃ 1) and 1046 < Lbol/(erg s−1) < 1048 (z ≃ 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).

Temperature. AGN activity strongly modifies the morphology of the T-field at every redshift (i) by boosting the occupied volume v₀ at increasing times for T > 10⁴ 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₁ at z = 1.5 − 2, especially for cold (T = 10⁴ K) and warm (T = 10⁵ K) isothermal surfaces, which are therefore globally rounder; (iii) by simultaneously increasing the convexity of warm and hot (T = 10⁶ K) domains, which are globally less flat at respectively z > 2 and z > 1 (larger integrated mean curvature v₂) and forming a network dominated by tunnels afterwards (negative v₃); (iv) by changing the sign of the integrated mean curvature for the regions with T > 10⁴ 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 ≤ 10⁴ 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³ k_B cm⁻³ K) is only marginally affected after z ≃ 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, intermediate-pressure domains (P = 10^3.9 k_B cm⁻³ K) instead shrink over time to a much larger extent, already since z ≃ 3, while keeping the same surface until z ≃ 1 as in presence of stellar feedback only. After z ≃ 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⁻³ 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₃), especially when most quasars of intermediate brightness is in place, which coincides with the epoch around the peak of the star-formation activity.

Total, HI, 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 dominate. AGN feedback modifies instead the morphology of high isodensity domains, where n = 0.1 cm⁻³, n_HI = 0.01 cm⁻³, and n_H2 = 0.1 cm⁻³, 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₃ indication of a larger number of tunnels. The sudden increase of v₃ of n and n_H2 (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 ≤ 0.01 Z_⊙) maintain their morphology across time, with minor effect by AGN feedback, which makes them globally smoother (smaller v₁) and slightly flatter (v₂ closer to zero) especially for z < 2 compared to the Z-field in NoAGN model, with a minor change to the topology at z ≲ 1.5. Conversely, the shape and topology of regions with higher metallicity, Z = 0.1 Z_⊙, and Z = Z_⊙, is strongly affected by AGN, especially for redshift respectively z ≲ 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₁, and form a network dominated by filaments, respectively, at z ≲ 2.5 and z ≲ 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 ≲ 1 warm regions with 10⁵ < T/K < 10⁷ resulting after smoothing from undistinguished 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 ≳ 10⁶ K. The T-field is further populated by colder cavities, roughly corresponding to collapsed haloes, which dominate the total volume at 1 ≲ z ≲ 2 with temperature decreasing with time down to 10³ 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 ∝ n^1.2 fitting the polytropic laws for isothermal and ultra-relativistic gases. Half of the volume is occupied by domains with pressure 1 ≲ P/k_B cm⁻³ K #x2272; 10³ 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 ≲ 10 k_B cm⁻³ 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₂ 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_H2-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 ≳ Z_⊙ and filaments with only slightly smaller metallicity when also AGN feedback is present, or with filaments with Z ≃ 0.1 Z_⊙ when only stellar feedback is active.

Altogether, MFs extend over the T − n, P − n, n_HI − n, n_H2 − 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². 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, H I density, H₂ 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_HI, n_H2, 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 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³. 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_H2 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⁵ K (10⁶ 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 cm⁻³ K and the n-field at 0.1 cm⁻³ are instead protected from a temporary topological transition occurring around z = 1.5 and stopping at z ≃ 1 if AGN feedback is active; the topology of the HI, H₂, and metallicity fields for n_HI = 10⁻² cm⁻³, n_H2 > 10^−2.5 cm⁻³, and Z > 0.1 Z_⊙ are also strongly altered by AGN feedback around z = 1.5 − 2. However, the topological and morphological transition of these excursions sets continue until z = 0.

The analysis of the bivariate distributions or phase diagrams T − n, P − n, n_HI − n, n_H2 − 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 thermo-chemical 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₀ and v₃) 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₀ and the large, positive value of v₃). 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 τ ∝ n_HI probed by Lyman-α forest in high-resolution quasar spectra (e.g. McQuinn 2016; Japelj et al. 2019), the X-ray luminosity L_X ∝ T³ probed observing galaxy clusters (e.g. Zou et al. 2016; Bahar et al. 2022; Poon et al. 2023), the Compton parameter Y ∝ ∫dℓ 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-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 complemented 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 complementary 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.

Acknowledgments

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 https://www.cosmicweb23.org. This work is partially supported by the grant https://www.secular-evolution.org 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

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_ν = {x ∈ D|f(x)≥νσ₀} have analytical expectation value given by (Tomita 1986)

$$\begin{array}{cc}\hfill v_0^{\mathrm{GRF}}(\nu )& =\frac{1}{2}\mathrm{erfc}(\nu /\sqrt{2}),\hfill \end{array}$$(A.1a)

$$\begin{array}{cc}\hfill v_1^{\mathrm{GRF}}(\nu )& =\frac{2}{3}\frac{\lambda }{2\pi }{\mathrm{e}}^{-{\nu }^2/2},\hfill \end{array}$$(A.1b)

$$\begin{array}{cc}\hfill v_2^{\mathrm{GRF}}(\nu )& =\frac{2}{3}\frac{{\lambda }^2}{{(2\pi )}^{3/2}}\nu {\mathrm{e}}^{-{\nu }^2/2},\hfill \end{array}$$(A.1c)

$$\begin{array}{cc}\hfill v_3^{\mathrm{GRF}}(\nu )& =\frac{{\lambda }^3}{{(2\pi )}^2}({\nu }^2-1){\mathrm{e}}^{-{\nu }^2/2}.\hfill \end{array}$$(A.1d)

Here erfc(x) is the complementary error function and the amplitude parameter λ² = σ₁²/3σ₀² is the ratio between the variance of the gradient of the field, σ₁² = ⟨(∇f)²⟩, and variance of the field itself, σ₀² = ⟨f²⟩_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

$$\begin{array}{cc}\hfill S_0& =\frac{{⟨f^3⟩}_c}{{\sigma }_0^4},S_1=-\frac{3}{4}\frac{{⟨f^2{\mathrm{\nabla }}^{2}f⟩}_c}{{\sigma }_0^2{\sigma }_1^2},S_2=-\frac{9}{4}\frac{{⟨(\mathrm{\nabla }f·\mathrm{\nabla }f){\mathrm{\nabla }}^{2}f⟩}_c}{{\sigma }_1^4}.\hfill \end{array}$$(A.2)

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

$$\begin{array}{cc}\hfill v_0(\nu )& =v_0^{\mathrm{GRF}}(\nu )+\frac{{\sigma }_0}{{(2\pi )}^{1/2}}\frac{S_0}{2}{H}_2{\mathrm{e}}^{-{\nu }^2/2},\hfill \end{array}$$(A.3a)

$$\begin{array}{cc}\hfill v_1(\nu )& =v_1^{\mathrm{GRF}}(\nu )+\frac{{\sigma }_0\lambda }{9\pi }(\frac{1}{2}{S}_0{H}_3+S_1{H}_1){\mathrm{e}}^{-{\nu }^2/2},\hfill \end{array}$$(A.3b)

$$\begin{array}{cc}\hfill v_2(\nu )& =v_2^{\mathrm{GRF}}(\nu )+\frac{2{\sigma }_0{\lambda }^2}{3{(2\pi )}^{3/2}}(S_0{H}_4+\frac{2}{3}{S}_1{H}_5+\frac{1}{3}{S}_2){\mathrm{e}}^{-{\nu }^2/2},\hfill \end{array}$$(A.3c)

$$\begin{array}{cc}\hfill v_3(\nu )& =v_3^{\mathrm{GRF}}(\nu )+\frac{{\sigma }_0{\lambda }^3}{{(2\pi )}^2}(\frac{1}{6}{S}_0{H}_5+S_1{H}_3+S_2{H}_1){\mathrm{e}}^{-{\nu }^2/2},\hfill \end{array}$$(A.3d)

where v_μ^GRF(θ) are given by Equations (A.1) and {H₀, …, H₅}={1, θ, θ² − 1, θ(θ² − 3),θ⁴ − 6θ² + 3, θ(θ⁴ − 10θ² + 15)} are Hermite polynomials. When the field f is convolved with a Gaussian kernel, as routinely done in numerical computation 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)≡exp(g − σ_g²)−1 related to a GRF g with vanishing mean and variance σ_g² = ln(1 + σ²), the excursion-sets {x ∈ D|g(x)≥νσ_g} and {x ∈ D|f(x)≥F(νσ_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):

$$\begin{array}{cc}\hfill {\nu }_{\mathrm{LN}}(\nu )& =\frac{ln\left[(1+\nu \sigma )\sqrt{1+{\sigma }^2}\right]}{\sqrt{ln(1+{\sigma }^2)}},\hfill \end{array}$$(A.4a)

$$\begin{array}{cc}\hfill {\lambda }_{\mathrm{LN}}(\lambda )& =\lambda {\left[\frac{\sigma }{1+{\sigma }^{2}ln(1+{\sigma }^2)}\right]}^{1/2}.\hfill \end{array}$$(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).

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 = 104 K, blue lines; T = 105 K, green; T = 106 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 σ0 = 0.3, (S​0, S​1, S​2) = (0.2, 1, 0.1), λ = 25.

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

Fig. B.1. Visualisation of the full box-size sections (8 Mpc/h thick and 50 Mpc/h wide in both directions) of the 3D gas temperature, pressure, density (total, HI, H2) 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.

Fig. B.2. Same as Figure B.1, but at z = 5.

Appendix C: Bias of smoothed fields

As stressed in the Guidelines (Sect. 3.1), the effective smoothing 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 outliers with |log T_MF − log T_SIMBA| > 1.5 is quoted in each panel). At z ≤ 1 (first and second column panels) a Gaussian smoothing 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^−4.5cm⁻³ (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.

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 TMF − log TSIMBA| > 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.

Depending on the feedback mechanism, the original P-field in SIMBA simulations is composed by three dominant populations 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, HI, and H₂ 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 n_HI-field, at all redshifts and for the n_H2-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 n_HI-field (not shown). The Gaussian smoothing kernel with σ = 4 cells totally erases sparse cells with density n_HI < 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. Simultaneously, this Gaussian kernel shrinks the high-density tail from the original 10⁻⁷ < n_HI/cm⁻³ < 1 range to 10⁻⁷ < n_HI/cm⁻³ < 10⁻², qualitatively preserving its shape (skewness). Interestingly, contrary to the T-field a higher spatial resolution obtained with a Gaussian kernel with σ = 1 cell is not sufficient to trace the original n_HI-field; it would totally deform the shape of the original probability distribution function, artificially populating the 10^{−10 < n_HI/cm−3 < 10−5} 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 therefore 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).

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_HI, n_H2, 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_* defining excursion-sets with noteworthy MFs: excursion-sets filling half the volume, v₀(f_*) = 0.5 (star); excursion-sets with vanishing integrated mean curvature, v₂(f_*) = 0 (diamond); excursion-sets dominated by isolated regions, namely, with maximum v₃(f_*) (triangle); excursion-sets dominated by tunnels, namely, with minimum v₃(f_*) (dot circle); and excursion-sets with cavities, namely, where v₃(f_*) has the second highest maximum (filled circle). We note that the v₃(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.

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 105). 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 nHI − n bivariate distribution at z = 0. Opaque symbols pinpoint threshold values where v0 = 0.5 (star), v2 = 0 (diamond), v3 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.

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.

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_HI parameter space (on average n_HI ≈ 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^2.8K and T > 10^6.8K), 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 ≤ 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 ≥ 1. Interestingly, without AGN feedback, the MFs defined by threshold values (n_*, T_*) marked by symbols 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_HI: 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_B cm⁻³ K and n > 10⁻⁴ cm⁻³ or n_HI > 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_B cm⁻³ K < 100 and n ≃ 10⁻⁷ cm⁻³ or n_HI ≃ 10⁻¹³ cm⁻³ at z = 0), in which the gas has on average P ∝ n^Γ, with a polytropic coefficient for diatomic molecules, Γ = 7/5; and low-pressure domains in low-density regions (P/k_B < 0.01 cm⁻³ 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_HI​ − ​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_HI ≃ 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_H2​ − ​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_H2 ≈ 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_⊙) 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_⊙ < 10, at the expense of domains with lower values that survive only in the absence of AGN feedback (the relationship between the H₂ 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 ≃ 3.

All Tables

All Figures

Fig. 1. Visualisation of full box-size sections (8 h−1 Mpc thick and 50 h−1 Mpc wide in both directions) of the unsmoothed 3D temperature, pressure, total density, HI density, H2 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, nH2, 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.

In the text

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 × 103 K (top middle), T ≥ 8 × 103 K (top right), T ≥ 1 × 104 K (bottom left), T ≥ 3 × 104 K (bottom middle) and T ≥ 5 × 104 K (bottom right). For illustrative purposes, 3D maps correspond to fields before applying additional Gaussian smoothing of 1.56 h−1 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.

In the text

	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 ≥ 7 × 102 K (top middle), T ≥ 104 K (top right), T ≥ 105 K (bottom left), T ≥ 106 K (bottom middle) and T ≥ 107 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.
In the text

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, Δ. 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 ≲ 2), but also by the X-ray heating at both high (z ∼ 4 − 5) and low redshift (z ≲ 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 ∼ 2 − 3.

In the text

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.

In the text

	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.
In the text

	Fig. 7. HI density. Minkowski curves for the neutral atomic hydrogen number density, or nHI-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.
In the text

	Fig. 8. H2 density. Minkowski curves for the molecular hydrogen number density or nH2-field, analogous to Figure 4. Its evolution is very similar to that of the total gas density; see Sect. 3.4.
In the text

	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 ≳ 0.01 Z⊙ for z < 1.5, while at earlier time the morphology is maximally influenced by stellar feedback, especially in regions with Z ≲ 0.1 Z⊙; see Sect. 3.5.
In the text

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/kB[cm−3K]) = {3, 3.9, 4.8}; log(n[cm−3]) = { − 3, −2, −1}; log(nHI[cm−3]) = {−4,−3,−2}; log(nH2[cm−3]) = { − 4, −2.5, −1}; log(Z[Z⊙]) = { − 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 1045 < Lbol/(erg s−1) < 1046; vertical lines mark the maximum of the AGN luminosity function for AGN with 1044 < Lbol/(erg s−1) < 1045 (z ≃ 1) and 1046 < Lbol/(erg s−1) < 1048 (z ≃ 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).

In the text

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 = 104 K, blue lines; T = 105 K, green; T = 106 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 σ0 = 0.3, (S​0, S​1, S​2) = (0.2, 1, 0.1), λ = 25.

In the text

Fig. B.1. Visualisation of the full box-size sections (8 Mpc/h thick and 50 Mpc/h wide in both directions) of the 3D gas temperature, pressure, density (total, HI, H2) 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.

In the text

	Fig. B.2. Same as Figure B.1, but at z = 5.
In the text

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 TMF − log TSIMBA| > 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.

In the text

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 105). 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 nHI − n bivariate distribution at z = 0. Opaque symbols pinpoint threshold values where v0 = 0.5 (star), v2 = 0 (diamond), v3 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.

In the text

	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.
In the text