© The Authors 2025

1. Introduction

In the current cosmological paradigm, baryons make up around 15% of the total matter in the Universe, while the rest is in dark matter (Planck Collaboration XXIV 2016). Dark matter is collisionless and interacts only through gravity. Thus, it is hard to detect but relatively easy to model analytically on large scales or numerically in the non-linear regime. On the contrary, baryons are relatively easy to observe, but they are subject to hydrodynamical forces and a variety of astrophysical processes, which makes them difficult to model theoretically, especially on small scales, where their effects are more significant.

The unknown mapping between gas and underlying dark matter poses strong challenges to fully extracting the cosmological information from ongoing and future large-scale structure (LSS) surveys. To constrain our cosmological models at the level of precision required, we need to take both these ingredients into account, or severe biases could appear in our cosmological analyses (e.g. Semboloni et al. 2011). For example, weak gravitational lensing (WL), which directly probes the cosmic matter field (The Dark Energy Survey Collaboration 2005; de Jong et al. 2012; Aihara et al. 2017), is particularly sensitive to the effects of baryons.

Nevertheless, since cosmic gas traces the underlying dark matter, it can be used as a powerful cosmological probe. For instance, galaxy clusters detected through the thermal Sunyaev-Zeldovich (tSZ, Sunyaev & Zeldovich 1972, 1980) effect or X-ray emission by their hot gas can be used to constrain the structure growth, the total neutrino mass, and the background expansion of the Universe (e.g. Bocquet et al. 2024; Ghirardini et al. 2024). Additionally, the kinetic Sunyaev-Zeldovich (kSZ, Sunyaev & Zeldovich 1972, 1980) power spectrum provides information on the large-scale gas mass and velocity field of the Universe, which is sensitive to primordial non-Gaussianities (Münchmeyer et al. 2019), neutrino mass (Mueller et al. 2015; Roncarelli et al. 2017), dark energy (Hernandez-Monteagudo et al. 2006; Bhattacharya & Kosowsky 2008), the law of gravity (Bianchini & Silvestri 2016; Roncarelli et al. 2018), and the duration of reionization (Zahn et al. 2005; McQuinn et al. 2005).

For all these scientific cases, it is necessary to carefully quantify the relationship between the spatial properties of gas and the underlying cosmological models, considering the impact of star formation, cooling, and feedback. This detailed knowledge is paramount with respect to optimally exploiting the data from current and next-generation cosmological surveys, increasing the robustness of cosmological analyses, or using it directly to constrain cosmological parameters.

Predicting the cosmic matter field in the presence of baryons is challenging. The gravitational evolution of small matter perturbations into a highly non-linear field can be accurately described by N-body simulations (for a review, see Angulo & Hahn 2022). These simulations have relatively well-understood convergence and range of validity (Power et al. 2003; Power & Knebe 2006; Knebe et al. 2009; Ludlow et al. 2019), and comparisons between different codes (Schneider et al. 2016) agree at the per cent level. However, baryons are subject to hydrodynamical forces and astrophysical processes that are much harder to model.

Hydrodynamical simulations evolve dark matter and gas particles under the simultaneous effect of gravitational and hydrodynamical forces. Unfortunately, simulations featuring cosmological volumes cannot resolve galaxy formation and astrophysical processes ab initio because of the vast dynamic range of scales involved. Thus, effective sub-grid prescriptions need to be adopted, for example, to model star formation, black hole seeding and growth, stellar winds, supernovae explosions, and active galactic nucleus (AGN) feedback, among others (for a review, see Vogelsberger et al. 2020). These astrophysical processes can impact the gas distribution up to relatively large scales. Consequently, the predictions of hydrodynamic simulations depend on the choices of the sub-grid prescriptions and their effective parameters (e.g. van Daalen et al. 2011).

The sub-grid parameters that regulate AGN feedback have the most significant impact on the large-scale distribution of the gas around high-mass halos (Ayromlou et al. 2023). These parameters vary the details of the energy injection from the accretion of supermassive black holes to the surrounding gas. Massive galaxy clusters retain the universal baryon fraction within their gravitational wells, even though they are not closed systems (see e.g. Mitchell & Schaye 2022). On the contrary, in galaxy groups, astrophysical processes are sufficiently energetic to eject gas from the halo, accumulating gas mass in the halo outskirts and making its density profile much more extended than that of dark matter (e.g. Velliscig et al. 2014; van Daalen et al. 2014; Davies et al. 2019; Tollet et al. 2019; Angelinelli et al. 2022; Sorini et al. 2022, 2024; Ayromlou et al. 2023). This mass redistribution has a non-negligible impact on the summary statistics of the cosmic matter field, dominated by the impact of AGN around high-mass halos on the scales of interest (e.g. Salcido et al. 2023).

The baryonic effects on the two-point statistics of the matter field have been extensively studied with hydrodynamical simulations (e.g. van Daalen et al. 2011, 2020; Mummery et al. 2017; Springel et al. 2018; Chisari et al. 2018; Hernández-Aguayo et al. 2023; Schaller et al. 2024). They can be modelled with flexible approaches such as baryonification (Schneider & Teyssier 2015; Schneider et al. 2019; Aricò et al. 2020, 2021), a halo model (Semboloni et al. 2011, 2013; Debackere et al. 2020; Mead et al. 2020), or other models that are directly based on hydrodynamic simulations (e.g. Salcido et al. 2023). The amplitude of the baryonic effect on the matter power spectrum seems to be tightly connected to the baryon fraction retained in high-mass halos (van Daalen et al. 2020; Salcido et al. 2023), even if secondary effects such as the shape of the gas density profiles need further investigations.

Nevertheless, the analysis of the cosmic gas statistics has received comparatively little attention. On large scales, Park et al. (2018) showed that in Illustris simulations (Vogelsberger et al. 2013; Nelson et al. 2015), gas is ‘anti-biased’ compared to the dark matter distribution. On small scales, Springel et al. (2018) showed that in IllustrisTNG simulations (Pillepich et al. 2018), gas clustering is suppressed; however, it is more clustered than the dark matter at higher redshifts. Mead et al. (2020) presented various auto- and cross-correlations of the gas field for BAHAMAS, which they used to calibrate the HMx halo model (Mead et al. 2021). More recently, Ni et al. (2023) compared the gas power spectra of ASTRID (Bird et al. 2022), IllustrisTNG (Weinberger et al. 2017), and SIMBA (Davé et al. 2019) models in small simulations of 25 h⁻¹ Mpc box size, finding that SIMBA (ASTRID) predicts the strongest (weakest) baryonic feedback among them. We note that, as pointed out by Schaller et al. (2024), these power spectra may not be converged due to the limited volume of the simulations, which likely tends to overestimate the baryonic suppression.

The impact of baryons on the cosmic velocity field has been studied even more seldom than gas density. Kuruvilla et al. (2020) found that gas pairwise velocities are suppressed on small scales relative to a gravity-only Universe in different hydrodynamic simulations. Strikingly, they observe that on large scales, r≥10 h⁻¹ Mpc, gas velocities are still suppressed by a constant factor of a few percent. Recently, Kwan et al. (2024) did not find any indication of a velocity bias in the radial velocity profiles around halos in the BAHAMAS simulations (McCarthy et al. 2017). Moreover, according to Kwan et al. (2024), the radial velocity profiles around halos of different masses are almost unaffected by feedback, both for the total matter and the gas. Contrary to Hellwing et al. (2016), who analysed the EAGLE simulation (Schaye et al. 2015), they found that baryonic physics leaves a non-negligible imprint on the redshift-space power spectrum of the total matter field, probably due to the stronger feedback in BAHAMAS compared to EAGLE.

This work aims to study the cosmic gas density and velocity fields relative to the case where only gravity is considered. We built upon existing works and extended them by leveraging the FLAMINGO cosmological-hydrodynamical suite of simulations (Schaye et al. 2023; Kugel et al. 2023). Due to their unprecedentedly large boxes and many sub-grid parameter variations, these simulations are well suited to studies of the impact of baryonic physics on cosmological scales.

Instead of focusing on any specific observable, we studied the impact of galaxy formation on several statistics of the gas field and identified the physical processes that drive those differences. We have sought to build a physical picture of the relevant baryonic processes that should be considered when modelling the cosmic gas density and velocity fields. We focus on the impact of high-mass halos, as they are expected to be the most significant for cosmological studies. We demonstrate that: (1) on large scales, k<0.1 h/Mpc, gas moves as dark matter, but it is less clustered than dark matter due to the most highly biased gas being converted into stars; (2) on small and intermediate scales, k>0.1 h/Mpc or r<10 h⁻¹ Mpc, baryonic effects suppress both gas clustering and velocities; (3) feedback changes the gas density and velocity profiles up to the outskirts of group-size halos; and (4) the outflowing gas reaches several virial radii and its properties depend on feedback strength as well as feedback implementation scheme.

The paper is organised as follows. We present the FLAMINGO simulations in Sect. 2. In Sect. 3, we study the cosmic gas density and velocity fields and characterise the main effects that drive the differences with respect to the gravity-only case. In Sect. 4, we focus on the astrophysical processes that intervene in halos of different masses, disentangling features that arise from different sub-grid parameters. Finally, in Sect. 5, we summarise our main results and present our conclusions.

2. Simulations

This work is based on the state-of-the-art FLAMINGO cosmological-hydrodynamical simulations (Schaye et al. 2023). The suite features different box sizes and resolutions. Here, we have employed the intermediate resolution, namely, N = 1800³ dark matter and baryonic particles and N = 1000³ neutrino particles in boxes of L = 1000 Mpc = 681 h⁻¹ Mpc. The average gas particle mass is m_g = 1.07·10⁹ h⁻¹ M_⊙, while the dark matter particle mass is 5.65·10⁹ h⁻¹ M_⊙ in the full-physics run and 6.72·10⁹ h⁻¹ M_⊙ in the gravity-only run. Unless stated otherwise, all the results in this paper are at the intermediate resolution and fiducial cosmology at z = 0. Details of the initial conditions, neutrino implementation, hydrodynamic solver and sub-grid prescriptions are given in Schaye et al. (2023). In the following, we briefly discuss the sub-grid models for AGN feedback due to their relevance for this work.

The FLAMINGO suite includes two different sub-grid models for AGN feedback: thermal and isotropic energy injection (Booth & Schaye 2009), and kinetic or jet-like energy injection (Huško et al. 2022; Schaye et al. 2023). Most of the simulations adopt the thermal implementation. While the black hole (BH) is accreting gas, a fraction (ϵ = 0.015) of the accreted rest mass energy is added to the BH energy reservoir. Once the BH has collected enough energy to increase the temperature of a gas particle by ΔT_AGN, the thermal energy is injected into the closest SPH particle to the black hole. Statistically, this becomes an isotropic energy injection. The main parameter controlling the strength of AGN feedback is ΔT_AGN, which is left free in the calibration.

The kinetic or jet AGN model implemented in FLAMINGO is a simplification of the sub-grid model presented in Huško et al. (2022). The same energy as in the thermal implementation is available to the BH, the only difference being how the energy is injected back into the surrounding gas. Once the energy reservoir of the BH exceeds $2\times \frac{1}{2}\times m_{\mathrm{g}}\times v_{\mathrm{jet}}^2$, two gas particles are randomly selected in a cone around the BH spin direction (each on one side), and kicked with velocity, v_jet. Analogously to the ΔT_AGN in the thermal model, v_jet, the velocity of the jet controls the AGN feedback strength, and its value is set in the calibration.

A powerful aspect of the suite is that it contains several astrophysics variations (i.e. choices for the sub-grid parameters) that are run with the same code, resolution, box size, and cosmology; thus, the impact of astrophysics can be isolated. The fiducial parameter set were calibrated to agree with the measured gas mass fractions in low-redshift galaxy clusters and the galaxy stellar mass functions at z = 0 (Kugel et al. 2023). Other variations were calibrated to reproduce deviations from observed gas fractions in clusters (e.g. fgas+2σ and fgas−8σ) or from the observed stellar mass functions (e.g. M_*−σ). We notice that the calibrations are not done ‘by hand’, but using Gaussian process emulators (noting that the calibration procedure is explained in detail in Kugel et al. 2023). Four sub-grid parameters, two controlling stellar feedback and two controlling AGN feedback, are fitted in the calibration; however, the parameters that control the gas fractions in clusters are mainly related to the AGN feedback efficiency. The parameters that control the stellar mass function are instead primarily related to supernova feedback (see Table 1 of Schaye et al. 2023). The suite also contains a ‘non-radiative’ run that was first presented in McCarthy et al. (2024). This run does not form stars and, consequently, does not feature feedback, which helps disentangle the effects of galaxy formation from other hydrodynamical processes¹. Moreover, the simulations with different AGN feedback prescription (but calibrated to the same observables) are useful for testing the impact of the adopted sub-grid model and its numerical implementation on the final observables.

Structure finding was done using the VELOCIRAPTOR code (Elahi et al. 2019), which builds on top of the friends-of-friends groups (FoF, Press & Davis 1982; Einasto et al. 1984; Davis et al. 1985), finding dynamically distinct groups by performing an iterative 6D FOF search in the phase space of the original FoF groups. The centre of the halos was established as the position of the most bound particle. Then, Spherical Overdensity and Aperture Processor (SOAP) was run on top of the previously identified halos and sub-halos to obtain 3D or projected properties, among other options (Schaye et al. 2023). In this work, we employ the spherical overdensity outputs (SO) with enclosed density Δ = 200ρ_m, where ρ_m is the mean density of the Universe. All the results shown are for isolated main halos without considering satellite sub-halos.

3. Gas statistics: mass distribution and velocity

In this section, we characterise the gas mass and velocity fields in the FLAMINGO suite. We study how non-gravitational physics impacts several summary statistics on different scales. In particular, we investigate the density power spectrum (Sect. 3.1), velocity divergence power spectrum (Sect. 3.2), and mean pairwise velocities (Sect. 3.3).

3.1. Density power spectrum

The impact of baryonic processes on the distribution of a given component is encoded in the ratio of its power spectrum with the gravity-only power spectrum, S_δδ(k) = P_δδ/P_{δδ, GrO}(k). The power spectrum is defined as $P_{\delta \delta ,i}(k)=\langle |\widehat{{\delta }_i}(k){|}^2\rangle$, where $\widehat{{\delta }_i}$ is the Fourier transform of ${\delta }_i=({\rho }_i-{\overline{\rho }}_i)/{\overline{\rho }}_i$, while ρ_i and ${\overline{\rho }}_i$ are the density and mean density of the corresponding species, respectively. If the ratio P_{δδ, i}/P_{δδ, GrO} is 1, the clustering is identical to the gravity-only run. All the power spectra were shot noise-corrected, following Eq. (B.11) from Mead et al. (2020). In the left panel of Fig. 1, we plot the results for different species: gas, stars and black holes, dark matter, and total matter (gas + dark matter + stars + black holes, dotted lines). The different colours denote the different FLAMINGO simulations.

Fig. 1. Ratio of density power spectra between hydrodynamic and gravity-only runs for several species. The colours indicate different FLAMINGO models. For comparison, we add the non-radiative run. Left panel: Results at z = 0 for gas (solid lines), dark matter (dashed lines), stars and black holes (dashed double-dotted lines), and total matter (dotted lines). Middle panel: Zoom of the large-scale suppression of four representative FLAMINGO models in the gas (solid lines) and total baryons (dashed-dotted lines) components at z = 0. Right panel: Redshift evolution of the large-scale gas bias, defined as the ratio between gas and gravity-only power spectra at k = 0.01 h/Mpc at each redshift.

The two main baryonic processes that shape the distribution of matter with respect to a purely gravitational Universe are (i) the formation of massive galaxies in the centre of halos and (ii) AGN feedback, which pushes gas outside halos. Consequently, the dark matter relaxes according to the modified matter distribution and gravitational field. The relevance of these baryonic processes varies depending on the matter species and scales we are looking at.

Stars and black holes (gas) are the most (least) clustered species on small scales. The strength of the AGN feedback mainly drives the clustering suppression of gas over the gravity-only run. For models calibrated to low gas fractions, such as fgas−8σ, the gas clustering is suppressed by 95% at k = 3 h Mpc⁻¹, while models with higher gas fractions, such as fgas+2σ, attain the same suppression only on smaller scales, k>10 h Mpc⁻¹. By k = 7·10⁻¹ h Mpc⁻¹, gas clustering is suppressed by more than 20% in all models considered.

The dark matter back-reacts to the combined effect of all the baryonic components. Its clustering differs from a gravity-only Universe by P_dm≈1.05−0.95 P_GrO on scales k<10 h Mpc⁻¹, being affected by both the creation of stars (enhancing the clustering) and the expulsion of gas (lowering the clustering). As gas is much more abundant than stars and black holes, the total matter clustering is suppressed between P_matter≈0.8−1.0 P_GrO.

Large-scale bias

On large scales, the dark matter and total matter clustering in the full-hydro run trace the gravity-only run. However, on the same scales (k≤2·10⁻² h Mpc⁻¹), the gas is biased low by a constant 4−5%, irrespective of the specific AGN feedback. Similarly, on the same scales, stars and black holes are biased high, namely, by 1.8−2.2. This is unrelated to the different masses of the gravity-only, gas, stellar, and black hole fields. As explained at the beginning of this section, the power spectra have been normalised, so their ratios would go to 1 when they have the same clustering, irrespective of the total mass of each species.

In the middle panel of Fig. 1, we look more closely at the large-scale gas bias. We plotted the suppression for four relevant FLAMINGO models: the fiducial model (solid blue, fgas), one with extremely low gas fractions in clusters (solid purple, fgas−8σ), and one with fewer stars (solid brown, M_*−σ). We also added the non-radiative run, which does not form galaxies and, therefore, it has no feedback (solid magenta). For the fiducial model, we added the suppression from a larger box with L = 2800 Mpc, but with the same resolution (thin green line); on scales of k = 10⁻² h Mpc⁻¹ the suppression is converged.

We observe that FLAMINGO models sharing the same stellar mass function converge to the same suppression of ≈5%, with sub-percent differences irrespective of the gas fractions (fgas vs fgas−8σ). The gas bias is slightly weaker in models calibrated to a reduced stellar mass function, reaching 4% (M_*−σ). These differences disappear in the non-radiative run or once we take into account all baryons (dashed-dotted lines): baryons are seen to perfectly trace the gravity-only clustering at k<10⁻² h Mpc⁻¹. On smaller scales, the differences between baryons’ and dark matter's primordial distributions are not completely washed out yet by gravity (Angulo et al. 2013) and a residual suppression of 1% remains, consistent with the linear prediction from CLASS at k<10⁻¹ h Mpc⁻¹ (black dashed line).

We interpret the large-scale gas bias as a necessary consequence of the large-scale stellar bias. As gas turns into stars in high-density regions, presumably halos (which are more clustered and biased compared to the dark matter distribution), the stars preferentially form and reside in more clustered environments relative to the dark matter field (with bias larger than one). Conversely, gas mass is reduced in those high-density environments, or, in other words, it preferentially resides in lower density and less clustered environments (with a bias of less than one). In agreement with this interpretation, the large-scale gas (stellar) bias is also smaller (larger) in the model calibrated to a reduced stellar mass function and there is no gas bias in the non-radiative run.

In the right panel of Fig. 1, we plot the redshift evolution of the large-scale gas bias, defined by the power-spectrum suppression at k = 0.01 h/Mpc at each redshift. The non-radiative run closely follows the linear theory prediction for the total baryons. Until z≈6, the gas distribution is in agreement with it, but once the stars start forming, the gas distribution becomes biased. At all redshifts, the bias in the M_*−σ run is lower compared to the fiducial run. The discrepancy between fgas and fgas−8σ is larger at higher redshift, which could be due to their stellar mass functions being different at higher redshift (noting that the calibration was done to reproduce the same stellar mass functions at z = 0, which does not guarantee the same at higher redshifts).

Interestingly, gas is most biased in all models at z≈1, becoming less biased at lower redshifts. We speculate that this could be again a consequence of the star formation history, the stellar bias becoming smaller at lower redshifts (e.g. Fry 1996; Tegmark & Peebles 1998; Springel et al. 2018), simultaneously reducing the bias in the gas field. We explore this more in Appendix A.

Earlier studies by Shaw et al. (2012) and Park et al. (2018) identified a significant gas bias in other hydrodynamic simulations, including Illustris (Vogelsberger et al. 2013; Nelson et al. 2015). We confirm their findings and extend them to the larger scales explored in the FLAMINGO simulations. Although the redshift evolution and the amplitude of the bias vary quantitatively, the qualitative trends and their interpretation remain consistent.

Thus, care is required when interpreting observations that employ gas as a tracer. For instance, in halo model approaches, the gas bias should be considered on top of the halo bias when analysing the thermal Sunyaev-Zeldovich (tSZ) power spectrum measurements (Bolliet et al. 2018; Salvati et al. 2018; Douspis et al. 2022), kSZ power-spectrum measurements (Shaw et al. 2012; Park et al. 2018), or halo bias-weighted mean electron pressure of the Universe (Chiang et al. 2020, 2021; Chen et al. 2023, 2024). The analyses based on hydrodynamic simulations inherently account for the gas bias (e.g. Battaglia et al. 2012a, b; McCarthy et al. 2014).

3.2. Velocity-divergence power spectrum

In linear theory, the peculiar velocity field follows a pure potential flow (i.e. without any rotational component), fully described by

$$\theta =\nabla ·v=-\delta \mathit{afH},$$(1)

where θ is the velocity-divergence, δ is the density contrast, a is the expansion factor, $f=\frac{\mathrm{d}lnD(a)}{\mathrm{d}lna}$ is the linear growth rate and H is the Hubble parameter. Thus, given the density field, the velocities are known: particles move convergently towards high-density regions and divergently away from low-density regions. However, the Universe is no longer linear at low redshifts, and the velocity field becomes uncorrelated with the density field as we approach small scales, creating a divergent component that is uncorrelated with the density field, and a curl component (Zhang et al. 2013). Gravitational collapse and several baryonic processes contribute to this decorrelation. In this section we focus on the latter².

We computed the velocity-divergence using the Delaunay-tessellation field interpolator (DTFE) technique (Schaap & van de Weygaert 2000; Cautun & van de Weygaert 2011). The DTFE package performs a Delaunay triangulation of particles and then interpolates the fields to a regular grid. Even though it has some limitations, velocity-divergence power spectra are well-behaved (Pueblas & Scoccimarro 2009; Hahn et al. 2015). However, as the velocity-divergence is a coarse-grained quantity, as we increase the resolution of the mesh, the value of the velocity-divergence can change dramatically. We computed the velocity-divergence and its power spectrum using a mesh of N = 256³ cells. By comparing the power spectra with higher and lower resolution meshes, we find that it is converged only at k<6·10⁻¹ h/Mpc, which defines the smallest scale we investigate.

In Fig. 2, we plot the ratio of the velocity-divergence power spectra of hydrodynamic and gravity-only simulations, S_θθ(k). The power in potential flows of gas (solid lines) is suppressed compared to that of dark matter, both for the dark matter in the full-physics run and in the gravity-only run (dashed lines). The dependence on the baryonic model is similar to for S_δδ. Hydrodynamical forces alone suppress potential flows in the non-radiative run, while galaxy formation and feedback processes further enhance this suppression. It is interesting to note that while the gravitational non-linearities act oppositely on densities and velocities (increasing the mass clustering and suppressing potential flows) (Pueblas & Scoccimarro 2009; Zhang et al. 2013; Zheng et al. 2013), the baryonic effects act in the same direction, suppressing both the clustering and potential flows.

Fig. 2. Ratio of velocity-divergence power spectra between hydrodynamic and gravity-only runs for several species at z = 0. The different colours indicate different FLAMINGO models. We display the results for gas (solid lines) and dark matter (dashed lines).

Furthermore, unlike density clustering, velocity divergence is not biased on large scales in any FLAMINGO model³. We notice that this aligns with the expectations from linear theory (Eq. 1): the gas feels the gravitational force from the entire matter field (not just the gas), which remains unbiased on large scales. Finally, the scale-dependent suppression of power starts on larger scales for the velocity divergence than for the density: for instance, for fgas−8σ at k = 4 h/Mpc, S_δδ≈0.8, while S_θθ≈0.6. Therefore, even if a fraction of the suppression of power in velocity-divergence can be understood as a result of different mass distributions in the full-physics and gravity-only runs (via Eq. (1)), it is clear that some baryonic processes further destroy the potential flows of gas and are efficient on very large scales.

Correlation between density and velocity divergence

Figure 3 illustrates the correlation between velocity-divergence and density for the gravity-only simulation (black lines) and gas in several FLAMINGO models (coloured lines). The density and velocity divergence are computed on a mesh of N = 256³ cells, with a cell size of 2.66 h⁻¹ Mpc using the Delaunay-tessellation field interpolator (Cautun & van de Weygaert 2011). The dark matter in the full-physics runs yields similar results to the gravity-only run, so we do not plot it for clarity. The dotted and dashed lines represent the 68.3 and 99.7 percentiles of 2d histograms, while the solid line shows the mean of the relation.

Fig. 3. Contours showing the correlation between velocity divergence and density contrast in the gravity-only (black lines) and gas in the hydrodynamic runs (coloured lines). Different FLAMINGO models are represented with different colours. The solid lines show the mean of the distribution, while the dotted and dashed lines define the 68.3 and 99.7 percentiles of the distribution, respectively.

At low densities, log₁₀δ<−1, the relation is tight for gas and dark matter, and the mean value is very similar – linear theory still holds. As we go to higher-density regions, the scatter around the mean increases, and by log₁₀δ>0, the mean starts to differ between dark matter and gas. High-density regions become shell-crossed, potential flows are destroyed, and transversal motions; namely, the velocity field curl component – arise (Zhang et al. 2013; Zheng et al. 2013; Hahn et al. 2015). This purely gravitational effect contributes to the scatter in the dark matter and gas. However, baryonic processes affect the gas on top of it. On the one hand, there is a clear distinction in the high-density tails of the distribution, visible in the P = 99.7 percentiles: the convergent motion into high-density regions stops faster for gas than in dark matter (at θ_gas≈−250 vs θ_DM≈−300), almost independently of the model, including the non radiative run without net cooling. This is consistent with the explanation that gas gets shocked when infalling into high-density environments, and its kinetic energy is converted into heat, contributing to the decrease of convergent motions. On the other hand, the dependence on the AGN strength appears in the divergent tail of the distribution (θ>0) at densities δ<10. In this regime, the velocity-divergence of the non radiative run agrees with the gravity-only case, while the rest of the models have increasingly divergent tails as feedback is stronger. We interpret this as gas that has been ejected by the AGN and exhibits positive radial velocities, namely, the gas is directed away from the high-density regions.

3.3. Mean pairwise velocity

Another common way to study the velocity field is using pairwise velocities. Pairwise velocities tell us how two particles are moving with respect to each other as a function of the separation scale. Even though the distribution is broad (Kuruvilla et al. 2020, see Figure 1), on average, particles are getting closer to each other due to the attractive nature of gravity. In other words, the mean pairwise velocity ω(r) is always negative,

$$\omega (r)=\frac{\langle (1+{\delta }_1)(1+{\delta }_2)(v_2-v_1)·\hat{r}\rangle }{\langle (1+{\delta }_1)(1+{\delta }_2)\rangle }<0,$$(2)

where r=r₂−r₁, δ_i=δ(r_i), v_i=v(r_i) and the average is over all particle pairs with a separation denoted as r. We notice that this is a mass-weighted quantity because we sample the velocity field only where particles are located. At leading order, Eq. (2) is then reduced to (a full derivation is given in Appendix B):

$$\omega (r)=[\langle {\delta }_1{v}_2\rangle -\langle {\delta }_2{v}_1\rangle ]·\hat{r}.$$(3)

This expression shows we are in fact dealing with the pairwise momentum rather than pairwise velocity. As defined here, both the density and velocity fields contribute to the pairwise velocities.

In Fig. 4, we plot the ratio of the mean pairwise velocity in hydrodynamic and gravity-only simulations, R_ω=ω(r)/ω_GrO(r). The mean pairwise velocities were computed using the HALOTOOLS package (Hearin et al. 2017) with some modifications. We computed the pairwise velocities among 100³ randomly selected particles. We checked that the pairwise velocities using this sample size are converged for r>0.7 h⁻¹ Mpc. We found that the gas pairwise velocities are suppressed on all scales, separated into two regimes. On small scales, r<10 h⁻¹ Mpc, pairwise velocities are slower in a model-dependent way, reaching to 25−30% suppression at r = 1 h⁻¹ Mpc. The models with stronger AGN feedback predict stronger suppression. At large scales, r>10 h⁻¹ Mpc, all models converge towards a constant bias ≈2−3%, though some model dependency remains. The pairwise velocities of the non radiative run are less suppressed on small scales and are not biased on large scales.

Fig. 4. Ratio of mean pairwise velocities between hydrodynamic and gravity-only runs for the gaseous (solid lines) and dark matter (dashed lines) components at z = 0. The colours represent different FLAMINGO models.

Similar results were found in other hydrodynamic simulations by Kuruvilla et al. (2020). Since the non radiative run does not have feedback, Kuruvilla et al. (2020) argued that the reason for the large-scale bias was that stellar feedback prevents gas from infalling. However, the interpretation of pairwise velocities in general and large-scale ‘velocity bias’ in particular is not straightforward. As Eq. (3) shows, one cannot tell whether gas is infalling slower than dark matter or if the measured bias is an effect of the mass-weighted nature of the velocity field.

We disentangle these two contributions in the upper panel of Fig. 5. Our approach consists of keeping the gas distribution but assigning the velocity field of dark matter. Effectively, per gas particle in the simulation, we find the closest dark matter particle and assign its velocity to the original gas particle⁴. This can be understood as sampling the dark matter velocity field in gas positions. Dashed-dotted lines show the resulting pairwise velocity. Similarly, the result for dark matter distribution with gas velocities is shown using dotted lines.

Fig. 5. Ratio of mean pairwise velocities between fiducial hydrodynamic and gravity-only runs at z = 0. The gas and dark matter ratios are displayed in both panels in solid and dashed green lines, respectively. Upper panel: Ratio of mean pairwise velocities with dark matter velocity field sampled at gas positions (dashed-dotted lines) and gas velocity field sampled at dark matter positions (dotted lines). Lower panel: Predictions from linear theory assuming different boosts for velocity and density power spectra, as indicated by the legend. The details of the procedure are explained in the main text.

On intermediate and small scales, r<10 h⁻¹ Mpc, both gas distribution (dashed-dotted lines) and gas velocities (dotted lines) differ from those of dark matter, suppressing the mean pairwise velocity of the gas. On these scales, gas particles indeed fall slower than dark matter. This is probably due to many effects, such as losing kinetic energy through shocks and AGN feedback. On large scales, r>10 h⁻¹ Mpc, the bias of gas pairwise velocities is dominated by the gas distribution: using dark matter or gas velocities (solid vs dashed-dotted lines) makes little difference. Conversely, there is a slight bias if one assumes the dark matter distribution, no matter the velocity field (dashed vs dotted lines). This means that, at large separations, gas particles infall towards each other at the same velocity as the dark matter. The bias in gas pairwise velocities mostly comes from weighting the (unbiased) velocity field with the (biased) density field. This interpretation is consistent with Kwan et al. (2024) and the analysis presented in Sect. 4.1, where we do not find any significant velocity bias in either the FLAMINGO or BAHAMAS simulations.

Mean pairwise velocities in linear theory

Using the insights from linear theory, we can make more quantitative statements about the relation between pairwise velocities and density and velocity fields. It can be shown that in linear theory (see Appendix B for the full derivation), we have

$${\omega }_l(r)\propto \int {\mathit{kP}}_{\delta \delta ,\mathrm{l}}(k)j_1(\mathit{kr})\mathrm{d}k,$$(4)

where P_{δδ, l}(k) is the linear matter power spectrum and j₁ is the Bessel function of the first kind. Taking into account:

$$-ik·\tilde{v}=\tilde{\delta }\mathit{afH}.$$(5)

Then, the following relation exists between velocity, velocity-divergence and matter power spectra:

$$P_{\mathit{vv},\mathrm{l}}(k)=P_{\theta \theta ,\mathrm{l}}(k)k^2=P_{\delta \delta ,\mathrm{l}}(k)k^2(\mathit{afH}{)}^2.$$(6)

Thus, we can re-write Eq. (4) as

$$\omega (r)\propto \int k\sqrt{P_{\delta \delta ,\mathrm{l}}{P}_{\theta \theta ,\mathrm{l}}}{j}_1(\mathit{kr})\mathrm{d}k\propto \int \sqrt{P_{\delta \delta ,\mathrm{l}}{P}_{\mathit{vv},\mathrm{l}}}{j}_1(\mathit{kr})\mathrm{d}k.$$(7)

As in Eq. (3), the dependence on density and velocity is explicit.

These expressions are strictly valid only in linear theory and for a collisionless fluid. However, we want to apply them to gas, which is subject to other forces and phenomena, such as pressure or feedback, on top of purely gravitational forces. To that end, we modify the linear-theory expression, introducing the “baryonic suppression” S(k),

$$S(k)=\sqrt{S_{\delta \delta }(k)S_{\mathit{vv}}(k)}=\sqrt{\frac{P_{\delta \delta }}{P_{\delta \delta ,\mathrm{GrO}}}}\sqrt{\frac{P_{\mathit{vv}}}{P_{\mathit{vv},\mathrm{GrO}}}},$$(8)

where now all power spectra are non-linear. We notice that with Eq. (6), S_vv=S_θθ. Because we multiply the linear theory power spectra with the ratio between non-linear power spectra, the modification is of order 1 (although we note this is an ad hoc term). Finally, the pairwise velocities are expressed as

$$\omega (r)\propto \int k{P}_{\delta \delta ,l}(k)S(k)j_1(\mathit{kr})\mathrm{d}k.$$(9)

We can directly test the impact of density and velocity fields in the pairwise velocity, changing the suppression factor in Eq. (8) via the S_δδ and S_θθ measured in simulations (Fig. 1). The predictions for different combinations are shown in the lower panels of Fig. 5 in grey. Assuming the gas mass distribution and velocities yields a fairly good description of pairwise velocities on all scales (solid line). Replacing gas velocities with dark matter velocities (dashed-dotted line) makes the suppression smaller on small scales. Thus, density and velocity fields impact the pairwise velocities on small scales by 10−20%. Equation (7) is based on linear theory and expected to fail on small scales. Hence, the small-scale behaviour should be interpreted qualitatively.

The large-scale behaviour is set by the choice of S_δδ: assuming the gas distribution yields a ≈2.5−3% bias (solid and dashed-dotted lines). In contrast, the dark matter distribution does not give a significant bias (dashed and dotted lines), independent of whether a gas or dark matter velocity field is used. Hence, the large-scale bias in the pairwise velocities can be explained by the large-scale bias in the gas power spectrum (Fig. 1).

The qualitative interpretation of such a bias (see Sect. 3.1) is consistent with this bias showing up in pairwise velocities too. Halos are clustered, and clustered regions fall towards each other faster than unclustered regions (Tinker 2007; Shirasaki et al. 2021). According to our interpretation of the large-scale bias in the gas power spectra, because gas forms stars in the centres of halos, it preferentially resides in lower-density regions, creating a bias b<1. In pairwise velocities, this implies that, relative to dark matter, the gas field weights halos or clustered regions less, that is, it weights less the regions with the largest (or most negative) pairwise-velocities, yielding on average lower (less negative) pairwise-velocities on large scales.

Although pairwise kSZ observations do not directly measure the mean pairwise velocity field of gas as computed here and are, instead, focussed on halo positions (Hand et al. 2012; Calafut et al. 2021; Chen et al. 2022; Gong et al. 2024; Li et al. 2024; Schiappucci et al. 2024), the gas bias may still introduce a notable effect. Following the same argument presented in Sect. 3.1, we found it helpful to consider that gas underweighs high-density regions and introduces a gas bias on top of the halo bias to interpret the measurements accurately. Nevertheless, in the case of pairwise kSZ measurements, the overall impact is expected to be smaller since it scales with the square root of the baryonic suppression (Eq. 8).

4. Gas density and velocity in halo outskirts

In the previous section, we analyse the summary statistics of gas density and velocity fields and found that several astrophysical processes significantly impact them. Notably, feedback shapes both fields on small scales. In this section we study the source of those effects: halos. We restrict the analysis to the FLAMINGO runs calibrated to varying gas fractions within clusters and exclude the variations of the stellar mass function.

As an example, Figs. 6a and 6b offer a visualisation of three halos of mass: M₂₀₀≈8.7·10¹² h⁻¹ M_⊙ (upper row), M₂₀₀≈9.92·10¹³ h⁻¹ M_⊙ (middle row), and M₂₀₀≈3.97·10¹⁴ h⁻¹ M_⊙ (lower row) in the fiducial run. The left and right columns display density and radial velocity, respectively, where blue colours denote positive (outflowing) radial velocities, and red colours negative (inflowing) radial velocities, normalised to $\pm v_{200}=\pm \sqrt{{\mathit{GM}}_{200}/r_{200}}$. Figure 6a corresponds to the dark matter, and Fig. 6b to the gas. We used the DTFE algorithm (Cautun & van de Weygaert 2011) to estimate density and velocities on a grid and then projected a slice of ±2r₂₀₀ along the line of sight direction. The cell size is 70, 80, 140 h⁻¹ kpc in each halo mass. A region of ±10r₂₀₀ around the halo is shown in each case, and the white circle denotes the radius of 5r₂₀₀.

Fig. 6. Density (left) and radial velocity (right) fields around three halos of mass M200 = 8.7·1012, 9.92·1013, 3.97·1014 h−1 M⊙ and radius r200 = 0.50, 1.12, 1.77 h−1 Mpc, from top to bottom. Figures 6a and 6b display the respective dark matter and gas fields of the fiducial (fgas) run. The panels show a region of ±10r200 around the halo centre, projecting a slice of thickness 4r200 in the z direction. The white circles represent a 5r200 radius. Positive (negative) radial velocities are shown in blue (red).

A quick visual inspection readily shows the impact of feedback on the density and velocity fields. Compared to dark matter, the most significant differences appear at the lowest mass halos, where AGN feedback can overcome the halo's gravitational potential and efficiently redistribute gas well beyond the boundary of the halo. As a result, the gas is more extended than the dark matter. A complementary picture appears in velocity space: while the bulk of the dark matter is infalling (red) and gets virialised within the halo (white central region), AGN feedback events push gas with positive velocities out of the halo (blue bubbles). In the fiducial FLAMINGO model depicted in the figures, the outflows reach up to ≈6r₂₀₀ with velocities of the order of v₂₀₀ in the lowest mass halo. In the following section we quantify such effects.

4.1. Density and radial velocity profiles

We computed the spherically averaged density and radial velocity profiles of gas around halos for different FLAMINGO models. The halo positions are defined as the minimum of the potential, while the halo velocity is the centre of mass velocity of the halo. We cross-matched⁵ the halos in the gravity-only and full-physics runs so that we compare the same halos among different models. Depending on the model, the halos may have lost some mass due to feedback, so we classified them in mass bins using their gravity-only mass. Moreover, we only considered isolated halos. A halo is isolated if no halo more massive than itself is found within 20r₂₀₀. This is an important criterion mainly for the lowest mass bin, where ≈90% of halos fall towards a larger one. The merging signal dominates in those cases, while the density and velocity profiles change dramatically. As we are only interested in the smooth gas component rather than in the effects of feedback in the merging history of halos, we imposed a somewhat extreme isolation criterion (for a similar approach, see Lau et al. 2015). Thus, we select halos preferentially living in low-density environments. Finally, we select 100 cross-matched halos in each mass bin and stack the profiles. In Appendix C, we show the results if no isolation criterion is imposed.

Figure 7 shows (i) the density; (ii) the ratio of the enclosed mass of gas and dark matter of the gravity-only run; and (iii) the mean radial velocity profiles in three halo mass bins, increasing halo mass from left to right. The mean radial velocity is computed using all the particles that fall in a given spherical shell. The profiles are computed per halo, and then we take the average of the profiles of all halos in each mass bin. Black dotted lines represent the gravity-only profile, and solid lines represent the gas. Gas density and mass profiles are normalised to their cosmic mass fraction f_gas for easier comparison. We checked that there are no noticeable differences between the dark matter on the full-physics run and that of the gravity-only run; hence, we only show the gravity-only result. Different colours denote different FLAMINGO models, as indicated by the legend.

Fig. 7. Gravity-only (dotted black lines) and gas (solid lines) spherically averaged radial profiles around isolated halos of three mass bins at z = 0. Top row: Density profiles, where the gas is normalised by the cosmic gas fraction. Middle row: Enclosed mass profiles over the gravity-only enclosed mass, where the gas is normalised by the cosmic gas fraction. Lower row: Radial velocity profiles over v200.

For all halo masses, stellar and AGN feedback push the gas away from the central region of the halo, lowering the gas density. In low mass halos, the ejected gas reaches r>r_{200, GrO}. More than half of the gas is ejected and remains outside r_{200, GrO}, and it is only around r≈10r_{200, GrO} that the enclosed gas-to-matter ratio approaches the cosmic value. At higher halo masses, feedback cannot overcome the gravitational potential of the halo and the differences among models are restricted to the central regions, r<r_{200, GrO}. However, the mass distributions of gas and dark matter are not the same at the boundary of the halo. This is visible in the ratio of their enclosed masses at the highest mass bin; specifically, around r≈(0.6−1.0) r_{200, GrO}, dark matter and gas seem to trace each other. However, at r>r_{200, GrO}, there is apparently more gas than dark matter once this ratio is normalised by the gas fraction.

Mean radial velocities show a consistent picture. In the highest mass bin, both gas and dark matter are being accreted by the halo with negative velocities. Outside the halo, dark matter and gas infall at the same velocity (Lau et al. 2015; Albæk et al. 2017). When they reach the boundary of the halo, the average velocity rapidly reaches zero. Dark matter and gas particles draw orbits around the halo centre, averaging inward and outward movements. On top of this, gas is collisional, and when it gets to the boundary of the halo, it feels ram pressure and is shocked. The combination of both effects slows the gas while the dark matter pierces through the halo without any resistance. As a consequence, gas velocities are less negative and converge faster to a mean of zero (e.g. Lau et al. 2015). Additionally, gas piles up at the boundary of the halo, creating the step-like feature in the enclosed mass profiles. Thus, in addition to using tSZ observations to measure the shock location in clusters (Baxter et al. 2021; Adhikari et al. 2021; Anbajagane et al. 2022; Towler et al. 2024), we speculate that kSZ observations could serve as a complementary tool for this purpose.

At lower halo masses, the average gas radial velocity profiles are shallower and model-dependent, becoming less negative the stronger the feedback (as parametrised through the gas fraction within r_500c). Relative to the underlying dark matter, on average, gas is accreted 50−80% more slowly. The range of scales affected is wider too, ≈[2−10]r_{200, GrO}. It is not straightforward to understand this. There could be several things happening. A possibility is that, as feedback ejects gas outside the halo (more or less effectively, depending on the feedback strength), the infalling gas encounters different environments when it arrives at the boundary of the halo. Thus, in stronger feedback scenarios, the infalling gas finds ejected gas at larger radii and gets shocked, losing its kinetic energy, driving the average radial velocity towards zero faster; namely, by reducing the negative radial velocity (infalling) tail of the distribution. This is consistent with the discussion of Kwan et al. (2024), where they argued that the multistream region of the gas is more extended than dark matter's at low halo masses. Another possibility is that the positive velocities of outflows contribute in a non-negligible way out to large radii and drive the average velocity towards zero – i.e. by increasing the positive radial velocity (outflowing) tail of the distribution. Most probably, it is a combination of both. In the next section we explore the role of the latter.

Nevertheless, there are clear differences between the velocity profiles of Kwan et al. (2024) and ours. First, the velocity profiles in Kwan et al. (2024) have a minimum around the virial radius, and a maximum at slightly larger radii than the minimum. At larger radii, the velocity profiles drop again, reaching even more negative velocities than at the minimum, mainly in low-mass halos. We have not been able to reproduce this feature, even though we speculate that it may be related to the isolation criterion (see Appendix C).

Second, Kwan et al. (2024) found that the several BAHAMAS variations with different feedback strengths do not impact noticeably the velocity profiles. However, in FLAMINGO the changes between different feedback strengths can get to 40% differences around 2r₂₀₀ (lower left panel of Fig. 7). This difference is probably a consequence of Kwan et al. (2024) not having isolated halos: the particles of the larger halo dominate in the computation of the mean radial velocity, and the differences in the smooth component become invisible (see Appendix C). Finally, Kwan et al. (2024) do not match halos across different simulations when comparing their velocity profiles. This likely hides the impact of feedback on individual halos too.

4.2. Outflows

In Fig. 8, we plot density and radial velocity fields around the same M₂₀₀ = 1.00·10¹³ h⁻¹ M_⊙ halo: dark matter in the gravity-only run (upper row) and gas in several FLAMINGO models (from second to last row), with increasing feedback strength. The lower panels show the jet implementation of AGN feedback, while the other runs use the thermal implementation. A region of ±10r₂₀₀ around the halo is shown, where r₂₀₀ = 500 h⁻¹ kpc. Both densities and velocities share colour normalization in all cases. In particular, radial velocities are coloured in the range ±v_{200, GrO}.

Fig. 8. Density (left) and radial velocity (right) fields around the same halo of MGrO = 1.00·1013 h−1 M⊙ and r200 = 500 h−1 kpc. The upper panel shows the halo in the gravity-only simulation, while all the others show the gas component of the halo in different feedback scenarios, with increasing strength from top to bottom. The panels show a region of ±10r200 around the halo centre, projecting a slice of thickness of 4r200. White circles represent a 5r200 radius. Positive (negative) radial velocities are shown in blue (red).

The outflows vary dramatically with feedback strength. On the one hand, the strength sets the scales the outflows reach: in the weakest feedback scenario of the FLAMINGO models (second row), they get out to ≈5r₂₀₀, while for the strongest feedback, they reach further than ≈10r₂₀₀. We notice that these outflows, identified by the radial velocity, extend to much larger scales than the high-entropy outflows previously found in simulations (Tremmel et al. 2019; Nobels et al. 2022), which are confined to ≈0.1r₂₀₀. In SIMBA, the outflows reach ≈4r₂₀₀ (Yang et al. 2024). On the other hand, the outflow speed also scales with feedback strength. In fgas−8σ, gas is outflowing with velocities comparable to v_{200, GrO}, while in fgas+2σ, gas is slower than ≈0.5v_{200, GrO}.

4.2.1. Anisotropic feedback

Feedback is not an isotropic process. The radial velocities clearly demonstrate that the outflows exhibit a biconical bubble structure. Despite the fact that the energy injection in the thermal AGN implementation is isotropic (Booth & Schaye 2009), the gas takes the path of least resistance and expands towards low-density regions, creating buoyant bubbles (e.g. Tremmel et al. 2019; Nobels et al. 2022). In the jet implementation, the injection is happens via the kicking of two particles in the direction of two cones with 7.5° opening angles around the BH spin axis (Huško et al. 2022; Schaye et al. 2023); however, that does not prevent them from losing the dependence of the initial direction at large scales and forming buoyant bubbles (Huško et al. 2024). In fact, visually, the bubbles in jet and thermal implementations do not show obvious differences.

The outflows transport gas from the halo centre towards the outskirts, making the gas less dense in the centre for increasing feedback strength. Presumably, due to the anisotropic structure of the outflows, the redistribution of gas is also anisotropic, transporting more gas to low-density regions than to high-density ones and making the overall gas distribution more spherical. In Fig. 8 the visual comparison of gas density around ≈5r₂₀₀ with increasing feedback strengths suggests so. We speculate that stacking kSZ measurements based on orientation could provide a new method to constrain feedback effects (see Hadzhiyska et al. 2025).

Similar structures have been found in other hydrodynamic simulations with different sub-grid models. In IllustrisTNG and EAGLE (Péroux et al. 2020; Truong et al. 2021b,a, 2023) the outflows were typically confined to the halo, while in SIMBA they reach ≈4r₂₀₀ with the jet-variant (Yang et al. 2024). These works define the major and minor axis of the galaxy, and analyse the anisotropy of feedback according to those directions. They analyse several gas properties, such as metallicity, temperature, X-ray luminosity, velocity and densities along the minor and major axis directions. All properties except the density are enhanced in the minor axis direction compared to the major axis due to feedback activity. Hence, they propose oriented X-ray observations to detect such bubbles in the real Universe, and forecast that they should lie on the detection threshold of eROSITA (Truong et al. 2021a, b; Yang et al. 2024). Alternatively, Martín-Navarro et al. (2021) reported a complementary signature in anisotropic quenching of galaxies in Sloan Digital Sky Survey (SDSS) data (Ahn et al. 2014). However, later works showed that the anisotropic quenching is a prediction of hierarchical structure formation rather than being an indicator of anisotropic AGN feedback (Karp et al. 2023).

Whereas previous works find a decrease in gas density in the outflow direction, we argue that the visual inspection of Fig. 8 suggests the opposite. There could be several possibilities for this discrepancy: on the one hand, we define the anisotropy with respect to the large-scale density distribution rather than the galaxy's orientation. However, those two are statistically correlated in observations (Rodriguez et al. 2022) and simulations (Rodriguez et al. 2023). On the other hand, the scales that we are probing are also different: typically residing within the halo versus at several virial radii. We will further explore and quantify this in future works.

4.2.2. Mass and radial scaling of outflows

To quantify the relative importance of outflows, we need criteria that identify them. One could define any particle with positive radial velocities as outflowing, as done in Muratov et al. (2015) and Ayromlou et al. (2024), for instance. However, this simple criterion makes the interpretation of the results more complex. Because both dark matter and gas are virialised within the halo, and they have some velocity dispersion associated with them, this criterion identifies as outflows ≈50% of the material within the halo, both in gas and dark matter components (white region – average velocity 0 – in the centre of the upper right panel in Fig. 8). Similarly, the material far away from the halo, probably falling into another halo (boundaries of the upper right panel in Fig. 8) is also considered as outflow⁶. Thus, at least in the small and large-scale limits, the results would not have a clear meaning (see Appendix D.2).

Our approach is to define outflowing or inflowing regions and then tag the particles in those regions accordingly. We do so by using only the kinematic information of gas and dark matter particles. Using the DTFE algorithm, we compute the radial velocity in a grid around each halo, with a cell size of l = 30r₂₀₀/128. From the lowest (highest) mass halos in our sample, this gives a cell size of 120 h⁻¹ kpc (420 h⁻¹ kpc). In Appendix D.1 we show that this cell size yields converged results. Unlike other methods, the DTFE algorithm provides a value for the velocity at each point in space. Additionally, computing velocities on a grid by construction removes a part of the intrinsic velocity dispersion. Finally, a grid cell is considered as outflowing if the following conditions are satisfied:

$$\begin{array}{cc}\hfill V_{\mathrm{gas}}(x,y,z)>0,& \hfill \\ \hfill V_{\mathrm{dm}}(x,y,z)<0,& \hfill \end{array}$$(10)

where V_gas(x, y, z) and V_dm(x, y, z) are the radial velocity grids of gas and dark matter in the full-physics run. Once we have identified the outflowing regions, we tag the particles that fall in those regions as outflowing. Finally, we describe outflows with the mass that they carry (M_out) and their median velocity (v_out) in a given radial bin. In Appendix D.2 we compare our findings with more standard outflow definitions.

The results are shown in Fig. 9. We take the mean and standard deviation of 100 isolated halos in each mass bin. We display the amount of mass in outflows over the total gas mass (upper row) and its velocity (lower row) in a given radial shell. In the left column, we compute this quantity as a function of halo mass at a fixed radius, while in the right column, we keep the halo mass fixed and change the radius. The different colours denote different FLAMINGO models.

Fig. 9. Radial and halo mass scaling of gas outflows at z = 0. Upper row: Gas mass fraction in outflows, defined as the ratio of the outflowing gas mass over total gas mass in a given radial shell. Lower row: Mean velocity of outflows in a given radial shell. Left column: Halo mass scaling of outflows at fixed radius r = 5r200, GrO. Right column: Radial scaling of outflows at fixed halo mass M200 = 1013 h−1 M⊙. The averages are taken over 100 isolated halos in each mass bin between M200∈(1013, 1015) h−1 M⊙, and the error bars denote the standard deviations. The colours stand for different FLAMINGO models.

We notice that even though the trends are robust to the details of the selection criteria, the exact numbers do change. For instance, by setting the minimum velocity to be considered as outflowing to 0 (Eq. 10), we effectively get an upper limit on the mass in outflows and a lower limit on their velocities. Furthermore, feedback from infalling halos is also expected to contribute to the total outflowing gas. The isolation criterion used in the halo selection likely reduces this contribution. As a result, we may underestimate the relative significance of outflows for halos of a given mass. Thus, we mainly consider the overall trends and relative differences between FLAMINGO models. Additionally, as mentioned above, in virialised regions, it is hard to reliably distinguish outflowing gas from regions with positive radial velocities within their orbit using only kinematic information. Therefore, we make a conservative cut and show the results only at several virial radii.

Overall, the selection shows the expected behaviour with halo mass and distance from the halo. At fixed feedback strength, outflows carry relatively more mass at low halo masses (top left) and closer to the halo (top right). As we go further from the halo, the outflows slow down (bottom right). It is not clear how the velocities scale with halo mass (bottom left): in some runs, they become slower in larger halos (fgas−8σ, Jet_fgas−4σ), in others they get faster (fgas−2σ, Jet_fgas), and in others, they are nearly constant. This is highly dependent on the distance from the halo centre too.

The stronger the feedback, the more gas is outflowing and the higher its velocity. This generally agrees with the qualitative trends inferred from Fig. 8. For low halo masses, where feedback is most effective, the percentage of outflowing gas over total gas mass can vary by a factor of a few between the weakest (fgas+2σ) and strongest (fgas−8σ) feedback runs at a fixed radius. Moreover, outflowing gas reaches larger scales in stronger feedback scenarios.

Interestingly, the radial and halo mass profile of the outflows in the Jet runs seems to differ from those in thermal AGN implementations. In particular, compared to the thermal AGN runs, the outflows in Jet_fgas−4σ are relatively less massive and slower in the inner regions (low halo masses) and then become relatively more massive and faster in the outer regions (high halo masses). This suggests that the radial profile and halo-mass scaling of outflowing gas (and, therefore, the impact of feedback on the total gas distribution around halos) depend noticeably on the AGN feedback implementation scheme, even for models calibrated to the same cluster gas fractions and galaxy mass function.

In Fig. 10, we show the average halo baryon fraction within r₂₀₀ as a function of halo mass for different FLAMINGO models. We computed the halo baryon fraction of isolated halos with M₂₀₀∈(10¹³−10¹⁵) h⁻¹ M_⊙ in 10 logarithmically spaced mass bins, randomly selecting 1000 halos in each bin. We selected only cross-matched the halos and classified them in mass bins using their gravity-only mass, as described in Section 4.1. We computed the halo baryon fraction as the baryonic to total halo mass within r₂₀₀. We notice that this plot is similar to Fig. 7 of Kugel et al. (2023); that is, it features the gas fractions as a function of the halo mass of the data used to calibrate the simulations. Due to the limited box size of the calibration runs, only halo masses up to M₅₀₀ = 2.5·10¹⁴ M_⊙ were used in the calibration. Hence, the halo mass scaling is shown to be at least partially driven by the calibration data.

Fig. 10. Average halo baryon fraction at r200 as a function of halo mass at z = 0. The results are shown for isolated halos with M200>1013 h−1 M⊙. Colours denote different FLAMINGO models.

For the same calibration data, it becomes evident that jets are comparatively less efficient at lowering the baryon fraction within r₂₀₀ at lower halo masses, in agreement with Fig. 9. However, as seen in the mass profiles of Fig. 7, the gas ejected by Jet_fgas−4σ reaches scales as large as the ejected gas in fgas−8σ. This is a simulation prediction: at a given halo baryon fraction, the extent of feedback differs for the thermal and jet AGN feedback implementations.

Complementary information is provided by previous works analysing the reach of the ejected gas in several hydrodynamical simulations (Tollet et al. 2019; Davies et al. 2019; Angelinelli et al. 2022; Sorini et al. 2022; Ayromlou et al. 2023). In particular, Ayromlou et al. (2023) found that different hydrodynamical simulations redistribute gas to different scales at fixed halo baryon fraction within r₂₀₀ – the so-called closure radius. A comprehensive analysis of the closure radius in FLAMINGO simulations will be presented in Denison et al. (in prep.).

5. Summary and conclusions

In this paper, we employ the state-of-the-art FLAMINGO simulations (Schaye et al. 2023) to study the impact of hydrodynamical forces and galaxy formation physics on the density and velocity fields of cosmic gas. This suite has been calibrated with observed gas and stellar fractions within galaxy clusters and groups (Kugel et al. 2023) and is therefore particularly suited to the purposes of this study.

In the first part of the paper, we have explored several clustering statistics of cosmic gas (e.g. the density power spectra, velocity-divergence power spectra, and mean pairwise velocities) and characterised the differences relative to expectations from gravity-only simulations. Furthermore, we have focussed on the properties of the gas around halos of mass between M₂₀₀∈(10¹³, 10¹⁵) h⁻¹ M_⊙, disentangling different physical processes at play, and quantifying the effect of sub-grid parameters and AGN feedback implementations. Our main findings are as follows:

• On large scales, the gas velocity is identical to that of dark matter (Fig. 2). However, the gas density field is biased low relative to that of dark matter even on scales as large as k = 0.01 h/Mpc (Fig. 1). At z = 0, this bias is b²≈5% in the power spectrum, with a maximum amplitude of b²≈8% at z≈1. Indirectly, this also creates a large-scale bias in the mean pairwise velocities of gas, which could be misinterpreted as a sign of gas moving more slowly than dark matter at scales as large as r≈10 h⁻¹ Mpc (Fig. 4). We show that the bias is almost independent of the gas fraction in clusters, however, does noticeably depend on the stellar mass function (Fig. 1). We interpret this as a consequence of star formation in dense and clustered regions, preferentially leaving gas in less dense regions.

• On small scales, gas is greatly affected by non-gravitational forces. The power spectrum of gas is strongly suppressed relative to the gravity-only case (Fig. 1), as are the potential flows (Figs. 2 and 3) and the mean pairwise velocities (Figs. 4 and 5). While this is arguably the combined effect of several non-gravitational processes, it strongly correlates with the feedback strength of the FLAMINGO simulations.

• In agreement with previous works, we find that in the FLAMINGO simulations, AGN feedback reshapes the density field by pushing gas away from halos, being most effective in group-scale halos with M₂₀₀≈10¹³ h⁻¹ M_⊙ (Fig. 7).

• Gas radial velocity profiles differ from the gravity-only profiles in two ways. First, in massive galaxy clusters, gas gets shocked at the halo boundary and loses kinetic energy, driving its mean radial velocity towards zero faster than dark matter's. Second, for group-scale halos (M₂₀₀≈10¹³ h⁻¹ M_⊙), the radial velocity profiles noticeably depend on feedback at radial distances as large as ≈10 r_{200, GrO}≈5 h⁻¹ Mpc (lower panels of Fig. 7).

• Feedback-dependent changes in the density and radial velocity profiles can be explained mainly in terms of large outflows sourced by the accretion onto supermassive black holes (Figs. 6 and 8). Irrespective of thermal or jet implementations of the AGN feedback, the surrounding gas eventually expands, creating buoyant bubbles that can reach tens of virial radii with radial velocities comparable to the escape velocity of the halo (for a visual impression, cf. Fig. 8). These bubbles show a clear biconical structure, implying that in FLAMINGO AGN feedback becomes anisotropic on resolved scales.

• The fractional mass in outflows is larger at lower halo masses and near the halo boundary, whereas it becomes less important at larger distances (Fig. 9). Baryonic effects on density and velocity scale with the AGN feedback strength.

• The radial profile and halo mass scaling of outflows is sensitive to the AGN feedback implementation (Fig. 9). Jets become less efficient than thermal AGN at lowering the halo gas fractions within r₂₀₀ for lower halo masses (Figs. 7, 9 and 10). However, they redistribute gas to larger distances (Figs. 7, 9 and 10).

In summary, the connection between cosmic gas and dark matter is modulated by galaxy formation processes across a wide range of scales. In particular, gas around group-size halos is especially susceptible to feedback processes, significantly impacting the density and velocity fields out to large scales. Thus, in principle, measurements of the gas density and velocity in the outskirts of these halos can serve as reliable proxies for AGN feedback.

To this end, the kSZ effect is a powerful probe of the gas distribution around group-sized halos. The kSZ signal is proportional to the line-of-sight integral of the gas momentum, so it directly measures the gas density and velocity fields without the need to model the gas pressure or temperature (unlike thermal Sunyaev-Zeldovich or X-ray measurements). Therefore, we find that kSZ analyses are equipped to provide independent constraints on AGN feedback, which can be compared and combined with other complementary datasets, informing hydrodynamic simulations conducted in future studies. This will be highly valuable for weak lensing surveys, as it helps isolate the effects of galaxy formation on the matter field and allows for the extraction of robust cosmological information.

Acknowledgments

We thank Daniele Sorini and Marcel van Daalen for their helpful comments. LO acknowledges the support of “la Caixa” Foundation (ID 100010434) for the fellowship with code LCF/BQ/DR21/11880028. REA acknowledges support from project PID2021-128338NB-I00 from the Spanish Ministry of Science and support from the European Research Executive Agency HORIZON-MSCA-2021-SE-01 Research and Innovation programme under the Marie Skłodowska-Curie grant agreement number 101086388 (LACEGAL). IGM acknowledges the support of the Science and Technology Facilities Council (grant number ST/Y002733/1). This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 769130).

References

Appendix A: Redshift evolution of gas bias

In this appendix, we further explore the redshift evolution of the gas bias. In Fig. A.1 we represent the gas bias, splitted on two components. In solid lines, we show the suppression of gas power spectra relative to gravity-only power spectra evaluated at k = 0.01h/Mpc. Unlike in Sect. 3.1, we normalise gas power spectra relative to the baryon mass. In dotted lines, the total gas mass is relative to the total baryon mass squared. That is, the gas bias, as shown in the right panel of Fig. 1, is given by the ratio of the two (solid/dotted) – if there was no bias, the dotted and solid lines should overlap. The colours represent different FLAMINGO runs, as indicated by the legend.

Fig. A.1. Split of the gas bias on its components: suppression of the gas power spectra (S(k = 0.01h/Mpc)) and the evolution of the total gas to baryon mass squared ($m_{\mathrm{gas}}^2/m_{\mathrm{b}}^2$). If the gas was not biased with respect to the gravity-only field, the solid lines should lie on top of the dotted lines. Solid lines: Redshift evolution of the suppression of gas power spectrum, as defined by the ratio of gas and gravity-only power spectra at k = 0.01hMpc. The gas power spectrum is normalised to the total baryon mass. In dashed lines, we plot the linear theory prediction for the baryon suppression. Dotted lines: Redshift evolution of the gas to baryon mass fraction squared, defined as the total mass in gas over total baryon mass. The colours denote different FLAMINGO runs.

First, we focus on the non-radiative run. Its gas mass fraction is always one because it does not form stars. However, at high redshift gas shows a small bias, consistent with the linear theory prediction coming from the different primordial clustering of baryons and dark matter (black dashed line) (for a detailed discussion, see Angulo et al. 2013).

The rest of the FLAMINGO runs do form stars. Thus, the fraction of their gas-to-baryon mass ratio decreases (dotted lines), starting around z≈4−5. Approximately, that is the time when the suppression in the gas power spectrum starts deviating from the linear theory expectation for the baryon fluid (coloured solid lines vs black dashed line). The ratio of both gives the gas bias (right panel of Fig. 1), increasing as more stars form in biased regions of the gas field at lower redshifts.

At low redshift z≈1, gas bias turns and becomes less biased (right panel of Fig. 1). In Fig. A.1 we see the reason for that: at lower redshifts than z≈1, stars still form, and the gas-to-baryon mass ratio keeps decreasing. However, the ratio of the gas and gravity-only power spectra is flattened.

We interpret this again as the consequence of the evolution of stellar formation. At lower redshift, the new stars form in less biased regions compared to high redshifts, just due to the redshift evolution of structure formation – at lower redshift, more and more ”common” (less biased) peaks cross the density barrier for collapse or star formation. This was already theorised in early works of galaxy bias (Fry 1996; Tegmark & Peebles 1998) and found in simulations in Springel et al. (2018). The indirect consequence is that the gas field also becomes less biased.

Appendix B: Mean pairwise velocity in linear theory

In this appendix we will derive the linear theory expression at leading order of the mean pairwise velocity, defined by the following equation

$$\langle \omega (r)\rangle =\frac{\langle (1+{\delta }_1)(1+{\delta }_2)(v_2-v_1)\rangle ·\hat{r}}{\langle (1+{\delta }_1)(1+{\delta }_2)\rangle },$$(B.1)

where r=r₂−r₁ and the average is over all particle pairs at separation |r|=r. Doing the multiplication,

$$\langle \omega (r)\rangle =\frac{\langle v_2-v_1+{\delta }_2{v}_2-{\delta }_1{v}_1+{\delta }_1{v}_2-{\delta }_2{v}_1+{\delta }_1{\delta }_2(v_2-v_1)\rangle ·\hat{r}}{\langle 1+{\delta }_1+{\delta }_2+{\delta }_1{\delta }_2\rangle }.$$(B.2)

We study each term individually:

• 〈v₁〉 = 0 and 〈v₂〉 = 0 because otherwise there would be a net motion.

• 〈δ〉=〈(ρ−〈ρ〉)/〈ρ〉〉 = 0.

• δ₁δ₂(v₂−v₁) is neglected at O(3).

• δ₁v₁ and δ₂v₂ have no real part in linear theory. This is clear when we change r₂ by r₁ in Eq. B.7.

• $\frac{1}{\langle 1+{\delta }_1{\delta }_2\rangle }\approx \langle 1-{\delta }_1{\delta }_2\rangle \approx 1$ at O(3).

Thus, at leading order O(2), the mean pairwise velocity reduces to

$$\langle \omega (r)\rangle =[\langle {\delta }_1{v}_2\rangle -\langle {\delta }_2{v}_1\rangle ]·\hat{r}.$$(B.3)

Let's expand the first term in the equation. Writing the density and velocity in Fourier space,

$$\langle {\delta }_1{v}_2\rangle =⟨\int d^{3}k{\tilde{\delta }}_{k}exp(ik·r_1){\left[\int d^{3}k\prime {\tilde{v}}_{k\prime }exp(ik\prime ·r_2){\delta }_D(k-k\prime )\right]}^{*}⟩$$(B.4)

Applying the Dirac delta:

$$\langle {\delta }_1{v}_2\rangle =⟨\int d^{3}k{\tilde{v}}_k^{*}{\tilde{\delta }}_{k}exp(-ik·r)⟩$$(B.5)

Moreover, in linear theory

$$-ik·\tilde{v}=\tilde{\delta }\mathit{afH},$$(B.6)

where a is the expansion factor of the Universe, f is the linear growth rate, and H is the Hubble parameter. We replace $\tilde{\delta }$,

$$\langle {\delta }_1{v}_2\rangle =\frac{1}{\mathit{afH}}⟨\int d^{3}k(-ik·{\tilde{v}}_k){\tilde{v}}_k^{*}exp(-ik·r)⟩$$(B.7)

As only the part in the radial direction survives, we can write ${\tilde{v}}_k={\tilde{v}}_k\hat{r}$. Thus, $k·{\tilde{v}}_k^{*}=k{\tilde{v}}_k^{*}cos\theta$, θ being the angle between $\hat{r}$ and $\hat{k}$. Besides, we can decompose differential d³k

$$\int d^{3}k={\displaystyle {\int }_0^{2\pi }}d\varphi {\displaystyle {\int }_0^{\pi }}d\theta sin\theta {\displaystyle {\int }_0^{\infty }}\mathit{dk}{k}^2=-2\pi {\displaystyle {\int }_{-1}^1}dcos\theta {\displaystyle {\int }_0^{\infty }}\mathit{dk}{k}^2,$$(B.8)

Putting everything together, we get

$$\langle {\delta }_1{v}_2\rangle =\frac{2\pi \hat{r}}{\mathit{afH}}{\displaystyle {\int }_0^{\infty }}\mathit{dk}{k}^2\langle {\tilde{v}}_k{\tilde{v}}_k^{*}\rangle {\displaystyle {\int }_{-1}^1}dcos\theta [\mathit{ik}cos\theta ]exp(-\mathit{ikr}cos\theta )$$(B.9)

which can be recast into the matter power spectrum in linear theory (see Eq. 6),

$$\langle {\delta }_1{v}_2\rangle =2\pi \mathit{afH}\hat{r}{\displaystyle {\int }_0^{\infty }}\mathit{dkP}(k){\displaystyle {\int }_{-1}^1}dcos\theta [\mathit{ik}cos\theta ]exp(-\mathit{ikr}cos\theta )$$(B.10)

Expanding the exponential

$$exp(-ik·r)=exp(-\mathit{ikr}cos\theta )=cos(\mathit{kr}cos\theta )-isin(\mathit{kr}cos\theta ),$$(B.11)

and saving only the real part, we get

$$\langle {\delta }_1{v}_2\rangle =2\pi \mathit{afH}\hat{r}{\displaystyle {\int }_0^{\infty }}\mathit{dkkP}(k){\displaystyle {\int }_{-1}^1}dcos\theta cos\theta sin(\mathit{kr}cos\theta )$$(B.12)

Finally, we do a change of variables u = cos θ:

$$\langle {\delta }_1{v}_2\rangle =2\pi \mathit{afH}\hat{r}\int \mathit{dkkP}(k){\displaystyle {\int }_{-1}^1}\mathit{duu}sin(\mathit{kru})$$(B.13)

Doing the integral on u,

$$\langle {\delta }_1{v}_2\rangle =-4\pi \mathit{afH}\hat{r}{\displaystyle {\int }_0^{\infty }}\mathit{dkkP}(k)\left[\frac{sin\mathit{kr}}{(\mathit{kr}{)}^2}-\frac{cos(\mathit{kr})}{\mathit{kr}}\right]$$(B.14)

where, introducing the Bessel function definition

$$\langle {\delta }_1{v}_2\rangle =-4\pi \mathit{afH}\hat{r}{\displaystyle {\int }_0^{\infty }}\mathit{dkkP}(k)j_1(\mathit{kr})$$(B.15)

We note that, 〈δ₁v₂〉=−〈δ₂v₁〉 due to the imaginary number in Eq. B.6. Thus, we have:

$$w(r)\hat{r}=\langle {\delta }_1{v}_2\rangle -\langle {\delta }_2{v}_1\rangle =2\langle {\delta }_1{v}_2\rangle .$$(B.16)

Furthermore,

$$w(r)=-8\pi \mathit{afH}{\displaystyle {\int }_0^{\infty }}\mathit{dkkP}(k)j_1(\mathit{kr}).$$(B.17)

Appendix C: Isolation criterion

In this section, we evaluate how selecting isolated halos versus the global halo population affects our results. Similar to Fig. 7, Fig. C.1 shows the spherically averaged profiles around halos in two mass bins: density profiles (top panels), enclosed mass profiles relative to the gravity-only expectation (middle panels), and radial velocity profiles (bottom panel). The dashed lines represent results without any isolation criteria, while the solid lines correspond to halos that have no more massive neighbour within r_{200, GrO}.

Fig. C.1. Gravity-only (black lines) and gas (coloured lines) spherically averaged radial profiles around halos of two mass bins at z = 0. The solid lines correspond to exclusively isolated halos, while in the dashed lines we do not impose any isolation criterion. Top row: Density profiles, where the gas is normalised by the cosmic gas fraction. Middle row: Enclosed mass profiles over the gravity-only enclosed mass, where the gas is normalised by the cosmic gas fraction. Lower row: Radial velocity profiles over v200.

In the highest halo mass bin, M = 6.31·10¹⁴h⁻¹M_⊙, where ≈60% of halos are isolated, we have checked that the density and enclosed mass, and radial velocity profiles match those of the global population. As halo mass decreases, differences become more noticeable. For r>r_{200, GrO}, isolated halos are less dense than the global population, indicating that isolated halos reside in lower-density environments by definition. The radial velocity profile at r>r_{200, GrO} is more negative for the global population.

The isolation criterion has the strongest impact on the radial velocity profiles. In the lowest halo mass bin, where only ≈10% of halos are isolated, radial velocities around the global halo population are much more negative than those around isolated halos at r>r_{200, GrO}. The global minimum of the radial velocity shifts to a larger distance, r≈4r_{200, GrO}, in both the gravity-only simulation (for dark matter) and the full-physics runs (for gas). Additionally, at the position of the velocity minimum around isolated halos, a distinct ’knee’ appears, most notably in the gravity-only run. Lastly, the impact of feedback on gas radial velocity profiles in the range r∈[1−10]r_{200, GrO} is partially smoothed out in the global halo population.

We interpret this as the merging signal of low-mass halos moving toward larger halos. Close to the virial radius, the particles that are indeed falling to the halo dominate the profile. At larger radii, though, the low mass halos themselves are falling into the potential of a more massive companion, which dominates over the small halo particles and creates misleading negative velocities ’towards’ the small halo in the centre (notice that the amplitude of the infall velocities scales roughly as $v_{200}\propto \sqrt{M_{200}/r_{200}}$).

Appendix D: Outflow selection

D.1. Convergence

In this appendix we test the convergence of our outflow selection criterion (Eq. 10). To do so, we double the resolution of the DTFE mesh in which we compute the velocity fields. We show the results in Fig. D.1 for the fiducial physics run. The halo mass and radial scaling of the outflows at the fiducial resolution (cell size $l=\frac{30}{128}{r}_{200\mathrm{Gro}}$) are converged within error bars.

Fig. D.1. Convergence of outflow selection criterion with cell-size. As in Fig. 9, but only for the fiducial fgas simulation. In solid and dashed lines, the outflows identified with $l=\frac{30}{128}{r}_{200,\mathrm{GrO}}$ and $l=\frac{30}{256}{r}_{200,\mathrm{GrO}}$ cell-sizes, respectively. The results are well converged within errorbars.

D.2. Definition

Previous studies in the literature (e.g. Muratov et al. 2015; Ayromlou et al. 2024) have defined outflowing particles as those with a positive radial velocity,

$$v_r=(v-v_h)·r>0,$$(D.1)

where v and v_h are the 3D velocities of the particle and the halo, respectively, and r is the radial vector between them. Figure D.2 shows how the mass in outflows and their average velocity scale with radius and halo mass using this definition (similar to Fig. 9).

Fig. D.2. Dependence on outflow selection criterion. As in Fig. 9, but with the outflow defined as any particle with a positive radial velocity.

Overall, the trends between different models are similar to those in Fig. 9: stronger feedback in simulations leads to a higher fraction of mass in the outflows and higher outflow velocities. However, Jets-especially the Jet_fgas−4σ model-show a slightly different scaling with radius and halo mass. Thus, the qualitative conclusions from Sect. 4.2.2 remain valid.

However, there are some differences between the two outflow definitions (Eq. 10 vs Eq. D.1). First, at all radii and halo masses, Eq. D.1 identifies a greater amount of mass as outflowing, and at higher velocities. For example, at r₂₀₀ = 5r_{200, GrO}, 40% of the mass is classified as outflowing at M≈10¹³h⁻¹M_⊙ in the fgas−8σ run, compared to 30% using our original criterion. This is expected, as Eq. 10 is more restrictive than Eq. D.1.

Second, the radial scaling of outflows is qualitatively different: at larger distances, more mass is classified as outflowing, and at higher velocities. This unphysical behaviour arises because some gas particles identified as outflows are actually inflowing toward other halos, and not outflowing due to feedback events. A similar pattern would emerge if we applied this method to dark matter particles. In our original criterion, we avoid this issue by requiring that, for gas to be considered outflowing, the dark matter in the same region must be inflowing (the second condition in Eq. 10).

All Figures

Fig. 1. Ratio of density power spectra between hydrodynamic and gravity-only runs for several species. The colours indicate different FLAMINGO models. For comparison, we add the non-radiative run. Left panel: Results at z = 0 for gas (solid lines), dark matter (dashed lines), stars and black holes (dashed double-dotted lines), and total matter (dotted lines). Middle panel: Zoom of the large-scale suppression of four representative FLAMINGO models in the gas (solid lines) and total baryons (dashed-dotted lines) components at z = 0. Right panel: Redshift evolution of the large-scale gas bias, defined as the ratio between gas and gravity-only power spectra at k = 0.01 h/Mpc at each redshift.

In the text

	Fig. 2. Ratio of velocity-divergence power spectra between hydrodynamic and gravity-only runs for several species at z = 0. The different colours indicate different FLAMINGO models. We display the results for gas (solid lines) and dark matter (dashed lines).
In the text

	Fig. 3. Contours showing the correlation between velocity divergence and density contrast in the gravity-only (black lines) and gas in the hydrodynamic runs (coloured lines). Different FLAMINGO models are represented with different colours. The solid lines show the mean of the distribution, while the dotted and dashed lines define the 68.3 and 99.7 percentiles of the distribution, respectively.
In the text

	Fig. 4. Ratio of mean pairwise velocities between hydrodynamic and gravity-only runs for the gaseous (solid lines) and dark matter (dashed lines) components at z = 0. The colours represent different FLAMINGO models.
In the text

Fig. 5. Ratio of mean pairwise velocities between fiducial hydrodynamic and gravity-only runs at z = 0. The gas and dark matter ratios are displayed in both panels in solid and dashed green lines, respectively. Upper panel: Ratio of mean pairwise velocities with dark matter velocity field sampled at gas positions (dashed-dotted lines) and gas velocity field sampled at dark matter positions (dotted lines). Lower panel: Predictions from linear theory assuming different boosts for velocity and density power spectra, as indicated by the legend. The details of the procedure are explained in the main text.

In the text

Fig. 6. Density (left) and radial velocity (right) fields around three halos of mass M200 = 8.7·1012, 9.92·1013, 3.97·1014 h−1 M⊙ and radius r200 = 0.50, 1.12, 1.77 h−1 Mpc, from top to bottom. Figures 6a and 6b display the respective dark matter and gas fields of the fiducial (fgas) run. The panels show a region of ±10r200 around the halo centre, projecting a slice of thickness 4r200 in the z direction. The white circles represent a 5r200 radius. Positive (negative) radial velocities are shown in blue (red).

In the text

Fig. 7. Gravity-only (dotted black lines) and gas (solid lines) spherically averaged radial profiles around isolated halos of three mass bins at z = 0. Top row: Density profiles, where the gas is normalised by the cosmic gas fraction. Middle row: Enclosed mass profiles over the gravity-only enclosed mass, where the gas is normalised by the cosmic gas fraction. Lower row: Radial velocity profiles over v200.

In the text

Fig. 8. Density (left) and radial velocity (right) fields around the same halo of MGrO = 1.00·1013 h−1 M⊙ and r200 = 500 h−1 kpc. The upper panel shows the halo in the gravity-only simulation, while all the others show the gas component of the halo in different feedback scenarios, with increasing strength from top to bottom. The panels show a region of ±10r200 around the halo centre, projecting a slice of thickness of 4r200. White circles represent a 5r200 radius. Positive (negative) radial velocities are shown in blue (red).

In the text

Fig. 9. Radial and halo mass scaling of gas outflows at z = 0. Upper row: Gas mass fraction in outflows, defined as the ratio of the outflowing gas mass over total gas mass in a given radial shell. Lower row: Mean velocity of outflows in a given radial shell. Left column: Halo mass scaling of outflows at fixed radius r = 5r200, GrO. Right column: Radial scaling of outflows at fixed halo mass M200 = 1013 h−1 M⊙. The averages are taken over 100 isolated halos in each mass bin between M200∈(1013, 1015) h−1 M⊙, and the error bars denote the standard deviations. The colours stand for different FLAMINGO models.

In the text

	Fig. 10. Average halo baryon fraction at r200 as a function of halo mass at z = 0. The results are shown for isolated halos with M200>1013 h−1 M⊙. Colours denote different FLAMINGO models.
In the text

Fig. A.1. Split of the gas bias on its components: suppression of the gas power spectra (S(k = 0.01h/Mpc)) and the evolution of the total gas to baryon mass squared ($m_{\mathrm{gas}}^2/m_{\mathrm{b}}^2$). If the gas was not biased with respect to the gravity-only field, the solid lines should lie on top of the dotted lines. Solid lines: Redshift evolution of the suppression of gas power spectrum, as defined by the ratio of gas and gravity-only power spectra at k = 0.01hMpc. The gas power spectrum is normalised to the total baryon mass. In dashed lines, we plot the linear theory prediction for the baryon suppression. Dotted lines: Redshift evolution of the gas to baryon mass fraction squared, defined as the total mass in gas over total baryon mass. The colours denote different FLAMINGO runs.

In the text

Fig. C.1. Gravity-only (black lines) and gas (coloured lines) spherically averaged radial profiles around halos of two mass bins at z = 0. The solid lines correspond to exclusively isolated halos, while in the dashed lines we do not impose any isolation criterion. Top row: Density profiles, where the gas is normalised by the cosmic gas fraction. Middle row: Enclosed mass profiles over the gravity-only enclosed mass, where the gas is normalised by the cosmic gas fraction. Lower row: Radial velocity profiles over v200.

In the text

	Fig. D.1. Convergence of outflow selection criterion with cell-size. As in Fig. 9, but only for the fiducial fgas simulation. In solid and dashed lines, the outflows identified with $l=\frac{30}{128}{r}_{200,\mathrm{GrO}}$ and $l=\frac{30}{256}{r}_{200,\mathrm{GrO}}$ cell-sizes, respectively. The results are well converged within errorbars.
In the text

	Fig. D.2. Dependence on outflow selection criterion. As in Fig. 9, but with the outflow defined as any particle with a positive radial velocity.
In the text