A new method to quantify the effects of baryons on the matter power spectrum

Future large-scale galaxy surveys have the potential to become leading probes for cosmology provided the influence of baryons on the total mass distribution is understood well enough. As hydrodynamical simulations strongly depend on details in the feedback implementations, no unique and robust predictions for baryonic effects currently exist. In this paper we propose a baryonic correction model that modifies the density field of dark-matter-only $N$-body simulations to mimic the effects of baryons from any underlying adopted feedback recipe. The model assumes haloes to consist of 4 components: 1- hot gas in hydrostatical equilibrium, 2- ejected gas from feedback processes, 3- central galaxy stars, and 4- adiabatically relaxed dark matter, which all modify the initial dark-matter-only density profiles. These altered profiles allow to define a displacement field for particles in $N$-body simulations and to modify the total density field accordingly. The main advantage of the baryonic correction model is to connect the total matter density field to the observable distribution of gas and stars in haloes, making it possible to parametrise baryonic effects on the matter power spectrum. We show that the most crucial quantities are the mass fraction of ejected gas and its corresponding ejection radius. The former controls how strongly baryons suppress the power spectrum, while the latter provides a measure of the scale where baryonic effects become important. A comparison with X-ray and Sunyaev-Zel'dovich cluster observations suggests that baryons suppress wave modes above $k\sim0.5$ h/Mpc with a maximum suppression of 10-25 percent around $k\sim 2$ h/Mpc. More detailed observations of the gas in the outskirts of groups and clusters are required to decrease the large uncertainties of these numbers.

cosmology provided the influence of baryons on the total mass distribution is understood well enough. As hydrodynamical simulations strongly depend on details in the feedback implementations, no unique and robust predictions for baryonic effects currently exist. In this paper we propose a baryonic correction model that modifies the density field of dark-matter-only N -body simulations to mimic the effects of baryons from any underlying adopted feedback recipe. The model assumes haloes to consist of 4 components: 1-hot gas in hydrostatical equilibrium, 2-ejected gas from feedback processes, 3-central galaxy stars, and 4-adiabatically relaxed dark matter, which all modify the initial dark-matter-only density profiles. These altered profiles allow to define a displacement field for particles in N -body simulations and to modify the total density field accordingly.
The main advantage of the baryonic correction model is to connect the total matter density field to the observable distribution of gas and stars in haloes, making it possible to parametrise baryonic effects on the matter power spectrum. We show that the most crucial quantities are the mass fraction of ejected gas and its corresponding ejection radius. The former controls how strongly baryons suppress the power spectrum, while the latter provides a measure of the scale where baryonic effects become important. A comparison with Xray and Sunyaev-Zel'dovich cluster observations suggests that baryons suppress wave modes above k ∼ 0.5 h/Mpc with a maximum suppression of 10-25 percent around k ∼ 2 h/Mpc. More detailed observations of the gas in the outskirts of groups and clusters are required to decrease the large uncertainties of these numbers.

Introduction
Cosmological structure formation is driven by the dominant yet invisible dark matter (DM) component acting as a self-gravitating collisionless fluid. The subdominant baryonic component is collisional and dissipative, what makes it considerably more complicated to model. For a long time it was considered sufficient to only predict the DM clustering while assuming galaxies to trace the centres of matter over-densities. In recent years, however, results from hydrodynamical simulations suggested that baryons affect the total matter distribution up to scales large enough to become relevant for cosmology [1,2].
Upcoming galaxy and weak-lensing surveys such as DES 1 , LSST 2 , and Euclid 3 have the potential to become the next leading cosmological probes [3][4][5]. They are complementary to the highly successful cosmic microwave background (CMB) measurements [6], potentially providing cosmological information from smaller scales at much later times. The success of weak lensing surveys, however, crucially depends on our understanding of how baryons influence the matter distribution, making the study of baryons a critical task for cosmology.
Numerical simulations with full hydrodynamic and radiative treatment have become increasingly popular in the last decade and major progresses have been made [7]. The fundamental difficulty of hydrodynamical simulations, with respect to the much simpler dark-matteronly simulations, is that the large-scale gas distribution is affected by small-scale events not resolved by the simulation. Active galactic nuclei (AGN) driven by black hole accretion or single supernovae explosions provide energy feedback mechanisms altering the formation of galaxies. As they are far too small to be properly included in simulations, they are usually implemented phenomenologically as sub-grid processes. These feedback mechanisms are poorly constrained by observations and severely restrict the predicability of hydrodynamical simulations for cosmology.
Recently, several analytical approaches have been proposed to complement results from hydrodynamical simulations and to enable a more direct understanding of how baryons affect large-scale structure formation. Examples are baryonic extensions of perturbation theory [8,9] or modifications of the halo model [10][11][12][13][14]. The downsides of these two approaches are that perturbation theory is only accurate at large scales, while the halo model has a rather limited overall accuracy.
In this paper, we propose yet another approach to include the effects of baryons on the matter density field. The baryonic correction model consists of an analytically motivated modification of dark-matter-only simulations. It assumes every halo to be made of four different matter components -dark matter, stars, bound gas, and ejected gas -with specific fractions and profiles that can be constrained with observations. Together, the components sum up to a total halo profile which differs from the initial DM-only profile. This allows us to define a Lagrangian mapping and displace particles in a DM-only simulation in order to recover the modified profile for every halo.
The baryonic correction model has the advantage of providing a simple and adaptable recipe and a direct link between observations of the gaseous and stellar profiles on the one hand and the total matter density field on the other hand. For example, fraction and radius of the ejected gas can be arbitrarily varied, illustrating how the efficiency of energy feedback affects the matter distribution. The model is not supposed to replace hydrodynamical simulations but to provide a fast and controllable tool to quantify the effects of baryons on the large-scale matter distribution.
Concerning the analysis of the matter distribution, we focus on the power spectrum, the prime statistical measure of the large-scale structure. As a requirement from upcoming weaklensing surveys, such as Euclid, the theoretical knowledge of the matter power spectrum has to be at percent precision at wave modes below k ∼ 10 hMpc −1 . In Ref. [15] we showed that it is challenging to obtain this level of accuracy with DM-only simulations. As for baryons, an additional modification at the ten percent level above k ∼ 0.1 − 1 hMpc −1 is expected. It is therefore of prime importance to further quantify or at least parametrise the effects of baryons on the matter power spectrum.
The paper is structured as follows. In Sec. 2 we introduce the model, providing details about the different model components and how they are constrained observationally. Sec. 3 is devoted to the details of modifying N -body simulations and shows an example case of a single cluster. We then discuss the effects on the power spectrum and define a simple fitting function for the baryonic effects in Sec. 4, before concluding in Sec 5.

Mimicking the effects of baryons
The baryonic correction model aims to mimic the effects of gas and stars on the halo profiles of cosmological N -body simulations. The idea is to replace the initial dark matter NFW profile with a combination of an adiabatically relaxed dark matter profile, a stellar profile, a hot gas profile, plus a component of ejected gas, where profiles and abundances are determined by observations. We start with an overview of the model, before providing details about the shapes of the profiles and the parametrisation of the different abundance fractions.

The model
The outcome of N -body simulations can be approximatively described by haloes (identified with a halo finder) in a smooth background component. Around each halo centre the density is given by the dark-matter-only (dmo) profile where ρ nfw is a truncated NFW-profile (with finite total mass M tot ) and ρ bg is the background density of non-collapsed matter, both defined in Sec. 2.2. The simplified description of Eq. (2.1) can be modified to account for the effects of baryons on the total density field. A model with baryons may consist of an adiabatically relaxed DM profile (y rdm ), a profile of bound gas in hydrostatic equilibrium (y bgas ), expelled gas due to feedback ejection (y egas ), and a stellar component from the central galaxy (y cgal ). All these components can be combined to obtain the final profile of the baryonic correction model (ρ bcm ), which has the form ρ bcm (r) = f rdm y rdm (r) + f bgas (M )y bgas (r) + f egas (M )y egas (r) + f cgal (M )y cgal (r) +ρ bg , (2.2) where the normalised individual profiles y χ and the corresponding abundance fractions f χ are described in Sec. 2.3 and 2.4. The mass dependences of the abundance fractions comes from the fact that heating and ejection of gas by the AGN and stellar feedback depends on the total halo mass. Gas is easily blown out of smaller haloes, while it stays trapped in the deep potentials of galaxy clusters.
The different density profiles can be integrated to obtain the enclosed mass where the subscript χ stands for both individual components and the total initial and final profiles (dmo and bcm, respectively). By construction the mass profiles Y χ are normalised so that where M ≡ 4π∆ 200 ρ c r 3 /3. This defines a Lagrangian mapping from the initial (dmo) to the final (bcm) density profile and allows to directly shift particles in the simulation output as in order to recover ρ bcm starting from ρ dmo . Here r i and r f refer to the initial and final radii from the halo centre. In Fig. 1 we illustrate the method by showing the density profile, mass profile, and displacement field of a schematic example with exaggerated baryonic effects. The left panel shows that the presence of baryons leads to a steepening of the inner and a flattening of the outer density profile, the former due to the central galaxy and the latter due to the bound and ejected gas components. The integrated mass profile (middle panel) is then used to define a displacement function (right panel) providing a Lagrangian mapping between the initial and final particle distribution around each halo. The displacement function is obtained by inverting the final mass profile and subtracting it from the reference radius as described in Eq. (2.6) and illustrated by the blue arrows in the middle panel.

Initial profiles
This section contains the definitions of the DM-only profiles introduced in Eq. (2.1) of the last section. We use a truncated version of the NFW-profile ρ nfw [16] to avoid mass divergence. This is important for a proper normalisation of the different mass components, as shown in Eq. (2.4). We follow Refs. [17,18] and use the profile where x = r/r s and τ = r tr /r s and where r s < r 200 < r tr are the scale radius, virial radius, and truncation radius, respectively 4 . We tested the profile with haloes from DM-only simulations and obtain best results for τ = 8c, where c = r 200 /r s is the concentration parameter. More about the truncated NFW profile as well as the integrated mass profile can be found in Appendix A.1.
The constant background componentρ bg represents all non-collapsed matter. In a simulation of a volume V and with the total amount of N p haloes,ρ bg is given as where the sum goes over all haloes in a simulation and M (i) tot is the total halo mass obtained by integrating over Eq. (2.8).

Final profiles
In this section we provide the definitions for the bcm profiles introduced in Eq. (2.2) of Sec. 2.1. We start with the gas and stellar components and then discuss their back-reaction effect on the dark matter profile.
The bound gas in galaxy clusters mainly consists of X-ray emitting ionised hydrogen. Assuming hydrostatic equilibrium as well as a polytropic form for the gas pressure, i.e. P = T ρ ∝ ρ Γ , the bound gas profile can be described as [19][20][21] y bgas (x) = y 0 where x = r/r s . The polytropic index Γ is fixed by assuming the slope of the gas profile to equal the one from the NFW profile at r 200 / √ 5, yielding [19] Γ eff (c) In our model the bound-gas profile is described by Eq. (2.10) until r < r 200 / √ 5 and follows the NFW profile (i.e. Eq. 2.8) further out. This is physically motivated as the gas is believed to be in hydrostatical equilibrium in the centre and to act like a collisionless fluid in the outskirts of the halo. Furthermore, the profile has been shown to agree well with both simulations [13] and observational data [20].
The ejected gas profile is supposed to capture all the remaining gas which has neither been transformed into stars nor is part of the bound hot gas component, but has been expelled from the galaxy due to strong energy feedback from the AGN. The shape of the profile and especially how far it extends into inter-galactic space influences the power spectrum at rather large scales, a good estimate of the density profile is therefore important. We assume that the AGN energy induces velocity kicks to all gas particles following a Maxwell-Boltzmann Furthermore, we assume the unbound gas to freely expand into space reaching the approximate distance r = vt 0 after a certain time t 0 6 . As a result, the Maxwell-Boltzmann distribu-tion naturally transforms into the ejected-gas profile where r ej = v cr t 0 determines how far the gas with critical velocity v cr is pushed into intergalactic space.
We will see later on that the choice of r ej is crucial for studying the large-scale density distribution, as it determines up to what maximum scales baryonic effects are visible. In Sec. 2.5 two models for the halo mass dependence of r ej are presented, one being heuristic and very simple and the other being physically motivated via the expected AGN energy input.
For the stellar profile of the central galaxy we follow Ref. [13] and use which behaves as a power law in the inner part and drops off exponentially well beyond the half-light radius R h . The relation R h = 0.015 r 200 is based on observations from Ref. [22] and has been shown to agree with simulations of galaxy clusters [13]. So far, we have defined specific profiles to model the baryonic fraction of the total matter distribution. There is, however, an additional dynamical modification of the DM component due to a back-reaction effect from the baryons which has to be taken into account. This effect is expected to contract the DM profile in the inner and expand it in the outer regions of the halo, the former because of additional mass from the central galaxy and the latter because of the missing gas in the halo outskirts.
The baryonic back-reaction effect can be included in the model by allowing for adiabatic relaxation (i.e. contraction and expansion) of the DM profile. Early models by the Refs. [23,24] assumed shells of DM to contract due to the presence of the central galaxy while conserving angular momentum, i.e. r i M i = r f M f (M f and M i being the total mass with and without central galaxy). More recent work by Refs. [25,26] showed that a better agreement with numerical simulations is obtained if the slightly modified relation is used instead. The authors of Refs. [25,26] argue that reducing the constant a to 0.68 (with respect to a = 1 for angular momentum conservation) corrects for the fact that the growth of the central galaxy is not an instantaneous process. In Ref. [27] it is furthermore shown that Eq. (2.15) is not only a good model for the adiabatic contraction induced by the central galaxy but can also account for the back-reaction due to missing gas in the halo outskirts. We follow Ref. [27] and write for the mass terms

Abundance fractions
By construction, the sum over all component fractions equals one, i.e.
While the relaxed DM fraction is a constant given by f rdm = 1 − Ω b /Ω m (where Ω b /Ω m is the cosmic baryon fraction), the stellar fraction and the gas fractions depend on the total mass of the halo. For the bound and ejected gas components the mass dependence comes from energy injection of AGN feedback. While in large clusters the gas may only be heated up by the AGN activity, it is ejected outside of the halo in the case of Milky-Way sized galaxies. The hot gas inside galaxy clusters can be observed via X-ray radiation. The signal mainly comes from within r 500 and is therefore probing the bound gas component (bgas). Based on X-ray observations, the bound gas fraction can be parametrised as follows where M c and β are free model parameters describing the mass scale and the slope of the hot-gas suppression towards small halo masses [13].
In the left-hand-side panel of Fig. 2 we plot observational X-ray data from the Refs. [28][29][30] together with the fit from Eq. (2.19). The bulk of the observations lie below the baryon fraction given by Planck (horizontal dotted line), showing that a significant fraction of the gas has been pushed out beyond r 500 . Furthermore, the clear trend in the data suggests that smaller haloes loose more of their gas due feedback effects. The trend is well described by Eq. (2.19) with the best fitting parameters M c = 1.2 × 10 14 M /h and β = 0.6. However, as the data is subject to significant statistical and potentially systematical errors, we probe wide ranges of values for M c and β as well as the possibility of scatter in this paper.
The fraction of stars in the central galaxy can be parametrised assuming abundance matching. Ref. [22] proposed a functional form based on [31] with updated parameters = 0.023, M 1 = 1.526 × 10 11 M /h, α = −1.779, δ = 4.394 and γ = 0.547. In the right-hand-side panel of Fig. 2 we plot the stellar fraction of the central galaxy from X-ray [22,28] and weak lensing data [32] given by star symbols and the pink shaded area. The abundance-matching curve from Eq. (2.20) is added as black solid line and agrees well with the direct observations as shown by Ref. [22].
Based on the observations of the bound gas and stars, it is now straight forward to estimate the ejected gas fraction by subtracting the former components from the total baryon abundance, i.e., The ejected gas is therefore assumed to make up for the remaining unobserved baryons expected from the Planck CMB measurement [6].

Ejected gas radius
We model the ejected gas component with a Gaussian profile defined in Sec. 2.3. This is motivated by assuming that the velocity-kicks from the AGN follow a Maxwell-Boltzmann distribution and escape the haloes with constant speed. What we have not done so far is to estimate the critical ejected gas radius (r ej ) defined in Eq. (2.13). Later on in Sec. 4, we will show that the value of r ej is important as it sets the minimum wave mode above which the power spectrum is suppressed by baryons. A natural assumption is to link the critical radius to the typical distance travelled by a gas particle with the halo escape velocity v esc . The escape radius is given by where we approximate the typical time-scale t 0 by half the Hubble time. The most simple physically motivated model is to assume with the free model parameter η a , where the ejection radius is simply proportional to the virial radius. We will call this heuristic approach model (A) throughout the paper. A more evolved model can be constructed by using the efficiency of gas ejection given by the observationally constrained expelled gas fraction (f egas ). The value of f egas determines how many gas particles obtain velocity-kicks beyond the escape velocity, which sets the width of the Maxwell-Boltzmann distribution. Integrating Eq. (2.12) and equating it with the gas fraction yields  which can be solved iteratively to determine r ej . This consists of a more evolved approach referred to as model (B) in this paper. It also has one free parameter η b to account for uncertainties in the approximative calculation of r esc . The free model parameter η χ (with χ = a, b depending on the model) does not only incorporate uncertainties related to the time-span t 0 of AGN activity, but it also accounts for the fact that gas does not escape with constant speed (as assumed in Sec. 2.3). The velocity decrease of gas escaping a spherical NFW potential is indeed proportional to the viral radius and can therefore be absorbed into the η χ -parameter.

Correcting N -body simulations
The baryonic correction model presented in the last section can be applied to outputs of Nbody simulations. The displacement function defined by Eq. (2.6) allows to shift particles around haloes, thereby modifying the initial NFW profiles to account for the effects of baryons.
In this section we present our simulations and we show how the matter distribution is modified by the particle displacement.

Setup of simulations
We use N -body simulations from previous work [15,33] (performed with Pkdgrav3  [36,37], assuming an overdensity criterion of ∆ vir = 200 times the critical density ρ c . In principle, the model is independent of the halo finding technique and also works for other common definitions of the halo mass. However, it crucially relies on halo concentrations and therefore requires an accurate measurement of the halo profiles. Haloes with fewer than 500 particles are not included in the analysis, as their concentration cannot be measured accurately. We have verified that this choice does not affect our final results significantly.

Case study of a cluster
We now illustrate the effects of the baryonic correction model with the example of a randomly selected cluster of mass M 200 = 10 14 M /h, viral radius r 200 = 0.77 Mpc/h, and concentration c = 3.2. The cluster profile of the initial N -body output is well described by an NFW profile plus a constant background component 8 .
We study three hypothetical cases for the cluster with different fractions for the bound and the ejected gas components: (a) nearly all the gas is bound, f bgas = 0.151, f egas = 0.005; (b) half of the gas is bound and half is ejected, f egas = f egas = 0.078; (c) nearly all gas is ejected, f bgas = 0.005, f egas = 0.151. For simplicity the dark matter and the stellar fraction are kept constant at f rdm = Ω c /Ω m = 0.839 and f cgal = 0.005.
The effects due to the different choices for the bound and ejected gas fractions are illustrated in Fig. 4, where the density profile, the mass profile, the displacement function, and the density map are plotted from top to bottom. The initial (dmo) and final (bcm) model profiles are shown as black and red lines, while the corresponding density and mass measurements from the simulations are given as black and red symbols. Faint non-continues red lines show profiles from the individual components. The displacement function is plotted as dashed line for negative and a solid line for positive values. Fig. 4 provides similar information than Fig. 1, however with realistic baryon fractions, making it more difficult to see the baryonic effects on the matter profile. The characteristic density increase in the inner and flattening in the outer part of the profile -due to gas being condensed into stars and pushed away by feedback effects -are nonetheless visible. The more gas is in the ejected component, the more the outer density profile is flattened, resulting in an increase of the displacement function. The cluster density maps at the bottom of Fig. 4 give a visual indication of the changes imposed on the N -body outputs. Close inspection allows to observe a slight displacement of satellite haloes in the outer parts of the cluster.
The bottom-line of this section is that the displacement function most crucially depends on the ejected gas fraction. For a cluster with only bound gas the required particle displacement is indeed more than an order of magnitude smaller than for a cluster where all the gas is ejected.

Convergence
It has been shown by previous studies [15,[38][39][40][41]   with very high particle numbers. In Ref. [15] we find minimum values M p = 10 9 M /h and L = 500 Mpc/h for particle number and box size to avoid nonlinear (non Gaussian) deviations from percent convergence. However, these numbers hold for an absolute convergence criterion and do not apply to relative convergence where the difference between a specific model and its corresponding baseline model are examined. It turns out that relative convergence can be achieved with smaller boxes and particle numbers, arguably because box-size and resolution effects are divided out. We illustrate the relative convergence properties of the baryonic correction model in Fig. 5, where ratios of power spectra from simulations with increasing box-size and particle numbers are shown at fixed physical resolution. While the very low resolution outputs with N = 64 (cyan) and N = 128 (magenta) deviate at five and two percent, the N = 256 (red) result is converged at sub-percent level over all scales with respect to the N = 512 result (blue). This rapid convergence of the relative power spectrum is very encouraging. It suggests that once a baseline model is simulated to high accuracy with large particle numbers, alternative models can then be simulated at much lower resolution in order to obtain relative deviations to the baseline solution. This is not only true for baryonic effects but should hold for any deviation from the standard model, such as nonlinear effects of massive neutrinos or modifications of gravity. vations finally allows us to set constraints on the baryonic power suppression and to come up with a simple, physically motivated fitting function.

Single component analysis
There are different baryonic effects shaping the matter power spectrum at different scales. While the central galaxy and the subsequent adiabatic contraction of the dark matter increases power at small scales, feedback effects such as the AGN activity tend to push the remaining gas away reducing power at larger scales. With the baryonic correction model it is straight forward to study these effects individually, which allows to gain physical understanding of how baryons shape the matter power spectrum. We do this by treating all but one component as passive non-modified matter components with no influence on the final result.
The outcome of the single component analysis is illustrated in Fig. 6, where the solid blue line shows the total (bcm) power spectrum and where the individual components correspond to the non-continues coloured lines. In general and as expected, the gas components lead to a power depletion at rather large scales, while the central galaxy increases power at the smallest scales.
The dominant effect on the power spectrum comes from the ejected gas component in form of a significant suppression at small wave numbers. The suppression is further amplified by the adiabatic relaxation of the dark matter which is mainly governed by the missing gas in the halo outskirts. The hot-gas component, on the other hand, only has a very mild effect on the power spectrum yielding an additional suppression of less than three precent. The central galaxy strongly affects the smallest physical scales inducing a significant boost of power at k > 10 h/Mpc. Future galaxy and weak-lensing surveys will probe the matter power spectrum up to k ∼ 1 − 10 h/Mpc. It is therefore of prior importance to fully understand the physics of the ejected gas component since it has the strongest effect on the power spectrum at the smallest wave numbers. In the next sections we will show that the ejected gas fraction and the typical ejection radius are key quantities capturing at what scales and how strongly the power spectrum is suppressed.

Hot versus ejected: the influence of the gas
In the last section we have established that the most stringent modification of the power spectrum comes from the ejected gas. Unfortunately, this is also the baryonic component with the weakest observational constraints, since the ejected gas has a shallow and extended density profile which is very difficult to observe directly. In principle observations of X-ray radiation and the SZ effect can be used to measure gas profiles in clusters and groups. While the X-ray observations predominantly come form the centres of clusters and are ideal to measure the gas fraction inside M 500 , the SZ signal probes the cluster outskirts and therefore the ejected gas profile.
In Sec. 2.4 we have briefly discussed how X-ray observations can be used to determine the fractions of the bound and ejected gas components (the latter indirectly via a subtraction from the cosmic baryon fraction). We now investigate the influence of these components on the matter power spectrum. We do this by following the parametrisation of Eqs. (2.19) and (2.21) with the model parameters M c and β describing the mass scale and slope below which (above which) the hot gas fraction (ejected gas fraction) is suppressed. Although the two parameters can be directly constrained with X-ray measurements of the gas fraction at M 500 , there is significant uncertainty due to statistical and potentially systematic errors. We therefore allow the parameters to freely vary within a reasonable parameter range (enclosing the observations) and study their effects on the power spectrum.
In the top panels of Fig. 7  The results presented in Fig. 7 are encouraging because they show that the suppression of the power spectrum mainly depends on the critical mass and not so much on the slope of the ejected/hot gas-mass relation (i.e. Eq. 2.19). This makes it easier to find a physically motivated fitting function connecting key values of the gas distribution to the power spectrum suppression (see Sec. 4.6).

Ejection radius: effects at the largest scales
The profile of the ejected gas component can only be determined with observations from the outskirts of haloes. Recently it has become possible to constrain the gas fraction beyond the viral radius with SZ measurements from Planck [42]. These results allow to tentatively distinguish between different values for the ejected gas radius (r ej ) at cluster scales. For smaller haloes, r ej can be determined with either model (A) or (B) defined in Sec. 2.5. While model (A) simply assumes a constant fraction r ej /r 200 over all mass ranges, model (B) accounts for the fact that gas is more efficiently ejected out of smaller haloes. We mainly focus on model (B) which is physically motivated and constrain the corresponding model parameter η b . However, further observations of the gas fraction around galaxy groups and small clusters are crucial in order to confirm the validity of model (B) and its predictions with respect to the power spectrum.
In the upper panels of Fig. 8, we plot the gas fraction from the baryonic correction model together with observational results based on Planck SZ data. The left panel shows measurements from the Planck collaboration [42], who combined SZ pressure profiles with extrapolated X-ray data from XMM-Newton. Different extrapolation schemes yield different results (labelled H1 and H2 in the figure, see Ref. [42] for more information) which were claimed to roughly bound the true profile [43]. The results are in general agreement with our model, however with too large error bars to show clear preferences for a certain value of η b . While model H1 prefers smaller ejection radii (i.e smaller η b ), model H2 seems very  inconclusive.
The right panel of Fig. 8 shows results from Ref. [44] based on SZ measurements from Planck and ROSITA data, using a different mass estimate based on the assumption of hydrostatic equilibrium. This leads to a generally higher gas fraction around the viral radius, especially for non-cool core clusters, pointing towards rather small ejection radii with η b < 0.5. Very recent results solely based on X-ray observations from Chandra confirm these rather high gas fractions [45]. However, Ref. [43] argue that assuming hydrostatic equilibrium plus not resolving small gas clumps may artificially boost the observed gas fraction.
In Appendix A.2 we compare the gas fraction profile of our model to hydrodynamical simulations from Ref. [43] and find good agreement at all investigated mass scales. The best match is obtained for η b ∼ 0.6 pointing towards efficient feedback in the simulation in slight tension with the observational results presented above. The good agreement between the simulations and the baryonic correction model further confirms the validity of the assumed gas profiles.
The bottom panels of Fig. 8 illustrate the ejection radius for a given value of η b (left) and the corresponding ratios of the power spectrum (right). The mass dependence of the ejection radius results from the fact that small haloes have shallower potentials allowing the ejected gas to expand further as described by model (B) in Sec. 2.5. The power spectra shown on the right-hand-side of Fig. 8 are sensitive to the choice of the ejected gas radius. Increasing the value of η b leads to a shift of the power suppression towards smaller k-modes (larger physical scales) while there is only a minor effect on the maximum power suppression. Hence, there is a simple trend between the ejection radius of gas around haloes and the largest scale where baryonic effects become important 10 .
Based on the established connection between power suppression and ejection radius, it is possible to estimate the largest scale where baryons affect the power spectrum. The SZ based measurements from Ref. [44] prefer values of η b 0.4 resulting in baryonic effects being negligible for modes k < 0.5 h/Mpc. Alternative modelling by Ref. [42] confirms this result for one (but not the other) approximation scheme. Further confirmation comes from direct X-ray measurements by Ref. [45] which also point towards small small ejection radii for clusters.
Although it is encouraging that current observations are able to constrain the model parameters of the ejected gas profile, the limits on the power suppression have to be taken with a pinch of salt. Indeed, there are two main caveats potentially affecting the result: First, there could be a bias because of unresolved gas clumping and the assumption of hydrostatic equilibrium, artificially boosting the observed values of the gas fraction [43]. Second, all current observations come from large galaxy clusters while the power spectrum is mainly influenced by group-size haloes. We use a physically motivated model to estimate the ejection radius-mass relation based on the cluster results, which ultimately requires observational confirmation. For the future, it is therefore crucial to find ways to measure gas profiles of Milky-Way to group sized objects in order to better predict the mass clustering of the universe.

The role of scatter
One of the advantages of the baryonic correction model is its adaptability to different approximations or new model assumptions. For example, it is straight forward to introduce a scatter term for the bound and ejected gas fractions. The observational data of Fig. 2 indeed shows a significant amount of scatter. Some of it should come from uncertainties in the data, but there should also be a physical contribution due to the fact that haloes have individual formation histories affecting the gas components.
In order to test the influence of scatter in the bound and ejected gas components, we add Gaussian scatter to the bound gas fraction (f bgas ). This means that, for every halo in the simulation output, f bgas is drawn randomly from a Gaussian distribution with standard deviation σ = 0.1, 0.2 around the mean given by Eq. (2.19).
The influence of the scatter on the power spectrum is shown in Fig. 9. While the left panel illustrates the 2-σ scatter-contours around the best fitting function for the bound gas fraction, the right panel shows how the power spectra is affected. Despite the significant scatter of σ = 0.1 (green) and σ = 0.2 (light green), the power spectrum stays nearly identical to the scenario without scatter (blue). This is very encouraging, as it shows that individual gas fractions from the specific halo formation histories do not influence the matter power spectrum significantly, strengthening the results presented here.

Redshift dependence
The effects of baryons on the matter power spectrum are intrinsically redshift dependent. At high redshift, prior to the typical formation time of group-sized haloes, the effects are negligible, but later on, they are expected to grow steadily due to repeated gas ejection via feedback mechanisms. Quantifying the redshift dependence is especially important for predicting the cosmic shear because the signal depends on the power spectrum between roughly redshift zero and two depending on the distribution of sources (see e.g. [4]). Implementing the full redshift dependence into the baryonic correction model might appear challenging, as parameters like M c and η b could vary with redshift. Pinning down the redshift dependence of these parameters would require high-redshift observations of the halo gas fraction, which currently do not exist. However, AGN activity is observed to peak around z = 2.15 [46], suggesting the the model parameters M c and η b might not evolve much below redshift two.
It is important to realise that the baryonic correction model leads to intrinsically redshift dependent power spectra even for the case of constant model parameters. This is due to the fact that at different redshifts the power spectrum measurement is dominated by haloes of different mass scales. For example, small clusters strongly contribute to the signal at redshift zero, but they are too rare to have an influence above redshift one. In this paper we study the redshift dependence (of the range z = 0 − 2) due to the growth of haloes while assuming redshift independent model parameters. This rather strong assumption ultimately requires testing with simulations and, ideally, observations, an endeavour we postpone to future work.
The bottom-right panel of Fig. 10 shows the redshift dependence of the power spectrum, illustrating how the baryonic suppression grows substantially between z ∼ 2 − 1, slows down at z = 1, and stops to evolve around z = 0.5. This qualitative evolution is in agreement with the redshift dependence observed in hydrodynamical simulations of Ref. [47] (but see Fig. 1 in Ref. [48]), where the growth of suppression is shown to be significant at z = 2 − 1 and small at z = 1 − 0. In both our model and the simulations the total baryonic suppression roughly doubles in scale between redshift two and zero, an encouraging agreement supporting our initial assumption of redshift independent model parameters.

Fits for the power spectrum with baryons
In the last sections we have established that the presence of baryons leads to both a power suppression at medium and a power increase at small scales, the former due to gas ejection and the latter due to the central galaxy. The total amount of suppression is driven by the typical halo mass scale (M c ) below which most of the gas is ejected, while the maximum scale (minimal k-range) where the suppression becomes visible depends on the parameter η b (fixing the r ej -mass relation). Furthermore, the suppression is growing with decreasing redshift as the signal is dominated by larger and larger haloes. These simple trends make it possible to determine a fitting function based on three parameters, M c , η b and z, and relating key observables of the gas distribution in haloes to the shape of the power spectrum. The fitting function has the following functional form: where G describes the suppression due to the gas and S the small-scale increase due to the central galaxy stars. The gas suppression is best captured by the function The influence of the stellar component S on the power spectrum is negligible below k ∼ 5 h/Mpc and stays subdominant until k ∼ 20 h/Mpc. As these small scales are not the prime target of this paper, we leave investigations of the dependence between S and the functional form of f cgal (M ) to future work. In Fig. 10 we show how the fitting function (grey lines) performs with respect to a selection of power spectra with varying M c , η b , and z. While the fit does not account for subtleties like the slight change in the total strength of the suppression for varying η b , it nevertheless provides a good match to all the different power spectra with errors staying well below five percent. This is a remarkable agreement, showing that the baryonic correction of the power spectrum is indeed mainly governed by the two parameters M c and η b as well as the redshift.

Comparison to previous work
First analytical estimates of how baryons affect the matter power spectrum have been carried out by Refs. [49,50] more than a decade ago. These authors investigated the influence of baryonic cooling [49] and hot inter cluster gas [50] finding percent effects on the power spectrum. Similar conclusions were obtained a few years later by Refs. [1,51,52] based on hydrodynamical simulations which included radiative cooling and star formation but no AGN feedback. Ref. [51] found a power suppression of about two percent beyond k ∼ 0.5 h/Mpc. We approximately recover their result using the fitting function of Eq. 4.1 with the parameters M c ∼ 2 × 10 12 M /h and η b ∼ 1.0. Such a low value for M c corresponds to a case with basically no ejected gas component, which is in agreement with the lack of AGN feedback in the simulations of Ref. [51].
The first predictions of the matter power spectrum including AGN feedback were performed by Ref. [2] about five years ago (see also [48,53]). Based on simulations capable of reproducing both X-ray and optical observations [47,54], they found a suppression of up to 25 percent with a deviation starting at k ∼ 0.3 h/Mpc. This is a dramatic increase of the effect with respect to earlier investigations and has important implications for the interpretation of current and future weak lensing surveys.
With the baryonic correction model and using standard model parameters matched to observations, we find a total power suppression roughly similar to the results from Ref. [2]. Concerning the fitting function of Eq. (4.1), the best match is obtained with parameters M c ∼ 5 × 10 14 M /h and η b ∼ 0.4 pointing towards a strong AGN activity able to push matter out of large clusters. We do not recover the exact shape of the suppression reported by [2]. In principle, closer agreement could be found by implementing a somewhat shallower ejected gas profile into the baryonic correction model. However, we want to stress that the shape of the ejected gas profile (and therefore the shape of the power spectrum) crucially depends on the interplay between gas and the AGN which is currently poorly understood 11 .
Further attempts to parametrise the baryonic effects have been carried out very recently by the Refs. [9,10,[12][13][14] applying analytical methods. While the works of [9,10,12,14] are all based on results from [2], Ref. [13] exhibits an independent investigation with a parametrisation of baryonic components similar to the one we use. The main differences are that Ref. [13] built upon the halo model and did neither assume an ejected gas profile (ejected gas is simply added to the linear background perturbations) nor did it include adiabatic DM expansion. They found a total power suppression which is similar to the one we obtain if we neglect adiabatic DM expansion. On the other hand, the largest scale (smallest mode) of suppression found by [13] does not agree with our predictions, most probably because of their simplified modelling of the ejected gas component.

Conclusions
The influence of baryons on structure formation needs to be quantified in order to extract cosmological information from upcoming full-sky surveys. In this paper we present a new approach to investigate how baryons affect the matter power spectrum based on a modification of N -body simulations. The baryonic correction model assumes every halo to consist of an adiabatically relaxed dark matter profile, a stellar component, as well as a bound and an ejected gas component, which sum up to a modified total halo profile. This profile can be recovered from an initial NFW profile by applying displacement functions to all particles around haloes in a N -body simulation.
The main advantage of the baryonic correction model is that observational measurements of the stellar and gas distributions in haloes can be used to estimate the influence of baryons on the total matter distribution without requiring further input about feedback mechanisms. The effects of single components, such as the ejected gas or central galactic stars, can be studied individually, clarifying the importance of different baryonic components on the total matter density field. Furthermore, the model simply perturbs the output of N -body simulations which are accurate and well controlled calculations of the gravity-only clustering in the nonlinear regime [see e.g 15,38,56]. The baryonic correction model is much more accurate than, for example, the halo model, where both the dark matter power spectrum and the baryonic corrections are estimated analytically. Perturbing N -body simulations is a promising approach to produce mock skies for future weak lensing surveys. This is not only true for baryonic corrections, the method is potentially applicable to various deviations from the standard cosmological model, such as the nonlinear effects of massive neutrinos or modified gravity.
In this paper we focus on the power spectrum, which is the prime statistical measure for the large scale structure and contains crucial information about cosmology. Future galaxy and weak-lensing surveys will require precise predictions of the power spectrum in order to be used as cosmological probes. It is therefore crucial to quantify what effects on what scales are expected from the presence of baryons. The baryonic correction model gives important insights to these problems which we summarise in the following: • The most important effect on the power spectrum comes from the gas ejected out of haloes. The ejected gas fraction per halo mass, f egas (M ), regulates how strongly the power spectrum is maximally suppressed with respect to the DM-only case. Based on the baryonic correction model plus X-ray observations of f egas (M ), we predict the power suppression to be 10-25 percent at scales beyond k ∼ 1 h/Mpc.
• Equally important than the ejected gas fraction is the ejected gas radius, i.e. the typical distance the gas is pushed out of the halo centres. It controls the maximum physical scale (minimum wave mode) where baryons start to affect the power spectrum. Recent measurements of the Sunyaev-Zeldovich signal in large clusters suggest that ejected gas might be around or just outside of the viral radius. Based on these observations the baryonic correction model predicts the power spectrum to be unaffected by baryons up to wave modes of k ∼ 0.5 h/Mpc. However, observational uncertainties are large and more data from galaxy groups is required to confirm this constraint.
• The remaining baryonic components have smaller effects on the power spectrum at medium scales. While the bound X-ray emitting gas only contributes weakly by no more than ∼ 3 percent, the stellar components leads to a strong power increase but only at very small scales beyond k ∼ 5 h/Mpc.
Due to the simple connection between galaxy properties and the power spectrum, it is possible to identify a physically motivated fitting function for the baryon-induced power suppression. We suggest a two-parameter fit (with additional redshift dependence), one parameter regulating the total amount of suppression and one the maximum scale (minimal wave number) where the suppression exceeds the percent level. The first parameter is linked to the typical mass above (below) which most gas is bound (ejected), while the second is related to the typical ejection radius of the gas. The two parameters have clear physical meaning and capture scale and shape of the power suppression to good accuracy. Finally, another interesting result of this paper is the rapid convergence rate for ratios of baryonic corrected versus DM-only power spectra. While absolute percent convergence requires simulations with a minimum box-size of L ∼ 1000 Mpc/h and a minimum number of N ∼ 4096 particles per dimension [15,38], relative ratios are converged for L ∼ 128 Mpc/h and N ∼ 256. We argue that resolution and box-size effects appear in both the corrected and the DM-only run and are therefore divided out for relative measures. The rapid convergence rate of ratios is very encouraging, because it means that one high-resolution baseline simulation plus different bias factors based on lower resolution might be sufficient to predict signatures of baryonic physics or alternative cosmologies for next generation surveys.  Figure 11. Upper left: hot gas fraction at M 500 from the simulations of Ref. [43] (green symbols) and from observations (grey symbols, same as in Fig. 7). Other panels: Comparison between the gas profile from the baryonic correction model and stacked gas profiles from the simulations (thick coloured lines) for different cluster masses.
profiles inside haloes (given by Eqs. 2.10 and 2.14) are able to reproduce results from cluster zoom simulations. The model for the gas outside of the viral radius (i.e the ejected gas profile of Eq. 2.13), however, has not been compared to simulations so far.
In this appendix we use results from hydrodynamical simulations of Ref. [43] to compare gas fractions in and outside the viral radius. This is an important consistency test for the baryonic correction model as the gas fraction directly depends on the underlying bound and ejected gas profiles. The simulations of Ref. [43] were run to study galaxy groups and clusters as well as the intra-cluster medium. They used an AGN recipe where the black hole accretion is coupled to the star formation rate. This simple model has been shown to produce identical results in terms of gas distribution than more elaborated AGN models [57] and is ideal for simulating large boxes with a large clusters sample.
In the top-left panel of Fig. 11 we show the simulated gas fraction at M 500 (green symbols in top left panel) which is in good agreement with observations (grey symbols, same as in Fig. 2) albeit with a somewhat shallower mass dependence (green dashed line) than the best fit to the observations (grey solid line, see Sec. 2.4). In the remaining panels we illustrate the mean gas profiles of clusters from the simulations (thick coloured lines) together with the predictions from the baryonic correction model depending on the value of η b (see Eq. 2.24). Different panels refer to different cluster masses from 10 14 to 5 × 10 15 M /h. Fig. 11 shows that the simulated gas profiles from Ref. [43] agree well with the model predictions, which validates our assumptions regarding the density profiles of the ejected gas component. The best agreement is found with the parameter choice of η b ∼ 0.6. This is in some tension with the observations illustrated in Fig. 8 of the main text, which seem to favour smaller values of η b 0.4.