OPTIMIZATION AND QUALITY ASSESSMENT OF BARYON PASTING FOR INTRACLUSTER GAS USING THE BORG CUBE SIMULATION

,


INTRODUCTION
Galaxy clusters play an important part in the current landscape of cosmological analyses (see e.g., Allen et al. 2011, for a review).Their abundance in mass and redshift is predicted by the halo mass function, which is highly sensitive to cosmological parameters -in particular to Ω  and  8 (e.g., Tinker et al. 2008), but also to dark energy properties and neutrino masses (e.g., McClintock et al. 2019;Bocquet et al. 2020;Euclid Collaboration et al. 2023).Thanks to this dependence, one can compare the mass and redshift distribution in a cluster catalog to the halo mass function, and thus infer constraints on cosmological parameters (see e.g., Pratt et al. 2019, for a review).Such cluster counting analyses have been successfully conducted using cluster samples constructed from a large number of datasets, providing tight constraints on cosmological parameters (e.g., Planck Collaboration et al. 2016a;Pacaud et al. 2018;Bocquet et al. 2019;Bocquet et al. 2023;DES Collaboration et al. 2020).
Cluster-based cosmology strongly benefits from the combination of multi-wavelength information.Galaxy clusters can be detected in sky surveys in a wide variety of wavelengths, including optical (e.g., DES Collaboration et al. 2020), X-ray (e.g., Liu et al. 2022), andmillimeter-wave (e.g., Planck Collaboration et al. 2016b;Bleem et al. 2020;Hilton et al. 2021).Each individual wavelength possesses distinct advantages and disadvantages for cluster science.For example, the thermal Sunyaev-Zeldovich effect (tSZ, Sunyaev & Zeldovich 1972) enables the detection of clusters at millimeter wavelengths through their imprint on the cosmic microwave background ★ e-mail:fkeruzore@anl.gov (CMB), and the construction of nearly mass-limited cluster samples (see e.g., Carlstrom et al. 2002;Mroczkowski et al. 2019, for reviews).Complementary optical surveys offer precious insights on cluster redshifts, and enable mass estimation through the gravitational lensing of background galaxies (see e.g., Umetsu 2020, for a review).Multi-wavelength studies, therefore, provide a way to mitigate the different systematics affecting cluster science in each individual frequency domain and to fully harness the constraining power of galaxy cluster cosmology (e.g., Postman et al. 2012;Sereno & Ettori 2015;Mantz et al. 2016;CHEX-MATE Collaboration et al. 2021;Perotto et al. 2022).
Cosmological simulations are of prime importance to investigate the characteristics of cluster observations.Largevolume simulations have been routinely used to study the halo mass function and its dependence on cosmology (e.g., McClintock et al. 2019;Bocquet et al. 2020), as well as general halo properties, such as accuracy of mass estimates (e.g., Becker & Kravtsov 2011;Gianfagna et al. 2021;Debackere et al. 2022;Muñoz-Echeverría et al. 2023), or scaling relation between mass and cluster properties (e.g., Angulo et al. 2012;Sembolini et al. 2013;Child et al. 2018).In addition, large efforts have been made to use simulations to create synthetic catalogs and sky maps, providing accurate mock realizations of cosmological surveys to calibrate cosmological analyses; for example, in optical wavelengths, the cosmoDC2 (Korytov et al. 2019;LSST Dark Energy Science Collaboration et al. 2021), MICE (Fosalba et al. 2015), Euclid Flagship (Potter et al. 2017), andBuzzard (DeRose et al. 2019) datasets; along with the works of Sehgal et al. (2010), and the WebSky (Stein et al. 2020) and Agora (Omori 2022) suites in the millimeter/submillimeter domain.Such datasets provide invaluable insights into cluster cosmology by enabling detailed studies of analysis and modeling systematics, and help pave the way for the exploitation of future cluster surveys.
Cosmological simulations can be broadly separated in two groups.On one hand, gravity-only −body simulations have been widely used to create synthetic datasets, owing to the efficiency of gravitational solvers and to our accurate understanding of gravitational processes.State-of-the-art gravity-only cosmological simulations include the Last Journey (Heitmann et al. 2021), Farpoint (Frontiere et al. 2022), Euclid Flagship (Potter et al. 2017), and Uchuu (Ishiyama et al. 2021) simulations (for a full review of gravity-only cosmological simulations, and a more extensive list of recent products, see Angulo & Hahn 2022).In parallel, hydrodynamic simulations are an increasingly active field of research, and are now reaching volumes large enough to investigate cosmology and largescale structure formation; recent examples including the Magneticum 1, Borg Cube (Emberson et al. 2019), EAGLE (Schaye et al. 2015), BAHAMAS (McCarthy et al. 2017), IllustrisTNG (Pillepich et al. 2018), MillenniumTNG (Hernández-Aguayo et al. 2022), and FLAMINGO (Schaye et al. 2023) simulations (see e.g., Vogelsberger et al. 2020;Crain & van de Voort 2023, for recent reviews of the field).The obvious advantage of hydrodynamic simulations over gravity-only ones is the inclusion of baryonic matter, facilitating the computation of cosmological observables, such as the thermodynamic properties of the gaseous intracluster medium (ICM).Nonetheless, this comes at a cost, as hydrodynamic simulations are much more computationally expensive (Vogelsberger et al. 2020;Angulo & Hahn 2022).Moreover, even in cluster-scale objects, baryonic physics is strongly affected by astrophysical scale processes, such as feedback by active galactic nuclei (AGN) and supernovae (see e.g., Angelinelli et al. 2023).Because these processes typically occur at physical scales smaller than (or on the order of) the mass resolution of cosmological simulations, and because the underlying physics is highly complex and still far from being completely understood, their modelling suffers from large uncertainties and is an active field of research.
In light of the challenges faced by cosmological hydrodynamic simulations, an alternative approach is to use gravityonly simulations and to model baryonic properties in postprocessing.This approach usually consists of using semianalytic models to infer observables from the distribution of dark matter in a gravity-only simulation.Its main advantage resides in the flexibility that is offered; because gravity-only simulations are computationally efficient, and semi-analytic models can be computed in post-processing at relatively low cost, baryon pasting can be used to mock up cosmological observables at a fraction of the cost of a hydrodynamic simulation.Moreover, model exploration with baryon pasting can be performed with a single gravity-only run, thereby eliminating the need to conduct individual hydrodynamic simulations for each model realization (as done in e.g., the CAMELS project Villaescusa-Navarro et al. 2023).
In this context, baryon pasting methods (also referred to as baryon painting) have emerged as a way to estimate ICM gas properties from gravity-only simulations.One of the most common approaches is the use of the gas model introduced by Ostriker et al. (2005) (and further developed by, e.g., Bode et al. 2007Bode et al. , 2009;;Shaw et al. 2010;Flender et al. 2016Flender et al. , 2017;;Osato & Nagai 2023) to model gas densty and pressure as a function of gravitational potential, assuming the gas to un-1 http://www.magneticum.org/dergo a physical rearrangement to follow a polytropic equation of state (see §3.2).Other gas models have also been used, relying on different physical assumptions (e.g., Battaglia 2016;Mead et al. 2020), along with advances in machine learning techniques to draw a mapping between the matter field in gravity-only simulations and baryon properties in hydrodynamic volumes (e.g., Tröster et al. 2019;Thiele et al. 2020;de Andres et al. 2023;Chadayammuri et al. 2023).Baryon pasting methods today constitute a cornerstone of synthetic cosmological datasets generation, having been used on various simulation suites to create multi-wavelength sky maps for cluster science (e.g., Sehgal et al. 2010;Stein et al. 2020;Omori 2022).In parallel, several families of techniques have also been developed to post-process gravity-only simulations to adapt their outputs to account for baryonic physics.Halo models have been used to modify summary statistics (such as matter power spectrum and halo mass function) to account for the presence of baryons (e.g., Mead et al. 2016Mead et al. , 2020)).Similarly, so-called baryonification (or baryon correction) models have been developed, focusing on rearranging particles according to analytical models of halo properties (e.g., Schneider & Teyssier 2015;Schneider et al. 2019;Aricò et al. 2020) to quantify the impact of baryons on various observables related to the total matter density field.
In this work, we investigate the performance of baryon pasting in reproducing the three-dimensional gas density and pressure distribution in cluster-sized objects.We utilize the Borg Cube simulation suite (Emberson et al. 2019), which consists of a pair of cosmological volume simulations.One of the simulations was evolved solely under gravity, while the other incorporated non-radiative hydrodynamics.Since the initial conditions of the two simulations are identical, this approach enables the direct comparison of the gas properties of halos inferred via baryon pasting with those obtained for the same halos evolved consistently with hydrodynamics.In particular, we focus on the aforementioned Ostriker et al. (2005) baryon pasting model, and seek to find the model parameters giving the best agreement between pasted gas and hydrodynamic simulations.Unlike many previous works that have performed baryon pasting using azimuthally symmetric gas profiles, we take advantage of the full information available in the Borg Cube by working with arbitrary potential distributions, thus enabling the creation of detailed realization of cluster gas.
It is worth noting that the Borg Cube hydrodynamic simulation is evolved using solely non-radiative hydrodynamics, and therefore does not include subgrid modelling.As such, feedback processes and radiative cooling are absent from the hydrodynamic data used in this study, and thus from the baryon pasting model we investigate.This makes our results approximate, as feedback -in particular from AGN -and cooling are known from observations to affect the distribution of gas in cluster-scale halos (see e.g., Donahue & Voit 2022, for a recent review), both at the population level (e.g., Battaglia et al. 2010) and at the scale of individual objects (e.g., Hlavacek-Larrondo et al. 2015;Ruppin et al. 2023).Nonetheless, given the aforementioned complexity of subgrid modelling, a solution to replicate non-radiative hydrodynamics from gravityonly simulations provides a useful tool with known caveats.These results will set a baseline on the expected performance of baryon pasting, and will allow us to isolate extensions describing subgrid modeling to address them separately.This work thus marks the first step in a longer series of studies, and the inclusion of more complex physics will be the main subject of future publications (in particular focusing on joint develop-ments of baryon pasting and hydrodynamic simulations).
This article is structured as follows.In §2, we provide an overview of the Borg Cube simulation and the data subsets we are using.We discuss the baryon pasting model and methodology in §3.In §4, we present the method used to evaluate the best-fitting baryon pasting model parameters as a function of redshift.We describe the assessment of baryon pasting performance in §5, focusing on three-dimensional gas properties and integrated tSZ signal.Lastly, we discuss the limitations and outlooks of this work in §6, and conclude in §7.
Notations -Quantities denoted with a 500 (200) subscript refer to the properties of a halo within a characteristic radius  500 ( 200 ), corresponding to the radius enclosing an average density 500 (200) times greater than the critical density of the Universe at the redshift of the halo.Similarly, quantities indexed with a "vir" subscript refer to the properties within the virial radius  vir , computed from the virial overdensity defined in Bryan & Norman (1998).Unless explicitly stated, distances and coordinates are expressed in proper units; comoving units will be noted with a 'c' prefix (e.g., cMpc for comoving Mpc).

DATA 2.1. The Borg Cube simulation
The Borg Cube runs (Emberson et al. 2019, hereafter E19) are a pair of cosmological simulations, each evolving 2×2304 3 particles in a comoving volume of (800 ℎ −1 cMpc) 3 .A set of initial conditions was created according to the WMAP-7 bestfitting cosmology (Komatsu et al. 2011), with parameter values reported in Table 1.In each run, half the particles are flagged as baryons, with particle mass   = 5.21 × 10 8 ℎ −1  ⊙ , while the other half are cold dark matter (CDM), with particle mass   = 2.56 × 10 9 ℎ −1  ⊙ , enforcing that the fraction of total mass made up by baryon particles is the cosmic baryon fraction Ω  /Ω  of the input cosmology (where The first simulation evolves the particles using the HACC code (Habib et al. 2016) and considers only gravitational interactions between all particles.The second run uses the CRK-HACC solver (Frontiere et al. 2023), in which all particles interact through gravitation, and the baryons are in addition subject to hydrodynamic physics.Each of these two volumes uses a softening length of  soft = 14 ℎ −1 ckpc.In the following, we refer to the first run as "gravity-only" (GO), and to the second as "non-radiative" (NR).More details on the Borg Cube simulation can be found in E19.

Massive halos in the Borg Cube
Our baryon pasting algorithm seeks to produce accurate realizations of thermodynamic properties of the ICM, i.e. the gas component of the most massive halos in the Universe.We use high mass objects in the Borg Cube halo catalog, described in E19, and summarized here.Halo finding and characterization is done in two steps within the HACC framework.First, halos are detected in particle snapshots using a friends-of-friends (FOF) algorithm on the CDM particles, with linking length  = 0.168.The position of the most bound CDM particle in each FOF detection is flagged as the center of the halo.Second, the spherical overdensity (SOD) properties of the halo are estimated by integrating outwards the properties of member particles (from both species) in concentric shells, starting at the halo center.Unlike in E19, the SOD properties are estimated out to three characteristic radii for each halo, corresponding to different overdensities:  500 ,  200 , and  vir .Mass-redshift coverage -In order to limit the amount of data being analyzed while still retaining relevant cosmological information, we work on a subset of five fixed-redshift full snapshots of the Borg Cube:  ∈ {0, 0.5, 1, 1.5, 2}.At each redshift, we focus on halos with masses in the gravity-only simulation  GO 500 ⩾ 10 13.5 ℎ −1  ⊙ , roughly corresponding to the mass scale conventionally treated as marking the transition between the galaxy group and galaxy cluster regimes (see e.g., Pratt et al. 2019).Table 2 summarizes the number of halos matching this criterion in each snapshot; distributions of the relevant properties of halos considered in this work are presented in Figure 1.
NR/GO Halo matching -In order to run comparisons of gas properties on a per-halo basis, we need to identify the counterpart to each gravity-only halo in the non-radiative simulation.We perform this matching in each snapshot based on physical separation and mass difference.For each halo in the gravity-only run, we compute the distance to all halos in the non-radiative run, accounting for the periodicity in boundary conditions.The closest NR halo with a mass  500 compatible with that of the GO halo within 20% is accepted as a counterpart if the distance between the two halo centers is smaller than 0.5×  500 .Because the simulations are evolved from the same initial conditions, the halo distribution is expected to be similar.This similarity is illustrated in Figure 2, showing the matter distribution in a subvolume of both simulations.We can see that halos of similar sizes can be found in the same physical locations in both simulations, which illustrates the validity of comparing halo properties between runs.Therefore, we expect the matching algorithm to yield a hydrodynamic counterpart for the vast majority of gravity-only halos, with a few exceptions (due to, e.g., slightly different merging histories in the GO and NR runs).The number of matched halo pairs per snapshot is reported in Table 2.For each redshift snapshot, we find that less than 2% of the gravity-only halos are not matched to a non-radiative counterpart.

ICM GAS PROPERTIES FROM GRAVITY-ONLY
PARTICLE DISTRIBUTION In this section, we describe the methodology used to infer the properties of the intracluster medium from gravity-only simulations.We follow the prescription introduced by Ostriker et al. ( 2005) and further developed in later works (see e.g.,  TABLE 2 Number of halos with  500 ⩾ 10 13.5 ℎ −1  ⊙ in each redshift snapshot of the Borg Cube used in this study.We report the number of halos within the mass range of interest in the GO run (middle column), and the number of halos matched as being present in both GO and NR simulations (last column, see text).Bode et al. 2007Bode et al. , 2009;;Shaw et al. 2010;Sehgal et al. 2010).

Data pre-processing
We start by isolating the particles belonging to each halo in both the GO and NR full particle snapshots.A particle is flagged as belonging to a halo if it is located at a distance smaller than 2× vir from its center 2. Once flagged, the particle properties are saved for processing.For each CDM particle, the properties consist of its mass, 3D position in the comoving volume, and 3D velocity.For each baryon particle in the NR run, the local gas density and thermal energy per unit mass are also saved.
The resulting particle output is projected on a regular grid using a cloud-in-cell (CIC) scheme (see e.g., Hockney & Eastwood 1981; Angulo & Hahn 2022, §5.1.2).For each halo, the cell size is set to  = 0.1 ×  500 , and the total grid size to 4×  vir .For gravity-only halos, we deposit the particle masses and squared velocities on the CIC grid, allowing us to compute the total matter density  tot and average kinetic energy per unit volume  2 tot .The distribution of gravitational potential  is then computed on the grid by solving the Poisson equation in Fourier space, using a non-periodic FFT by zero-padding the grid by twice its original size on each side (as prescribed by Hockney & Eastwood 1981).For halos in the non-radiative run, CDM and baryon particles are projected independently.The former are projected following the same procedure as for the particles in the GO run.For baryon particles, we directly deposit the SPH baryon density at the position of each particle,  bar , and also store the corresponding pressure,  bar , computed as: where  = 5/3 is the adiabatic index of a monoatomic gas, and  bar is the thermal energy per unit mass of each particle.

Gas injection model
To paste gas onto gravity-only simulations, we follow an adapted version of the prescription of Ostriker et al. (2005) in the case of an arbitrary potential distribution.The prescription is based on the assumption that the gas density and pressure are linked through a polytropic equation of state, and depend on the gravitational potential distribution.They can be computed through the following steps.

Initial state
The gas properties are supposed to initially be determined by a constant gas fraction across the halo  g (the value of which will be discussed in §4.1).The gas density and pressure in each cell  of the grid are given by: Using this matter distribution, we compute the gas radius  g , defined as the radius enclosing a gas mass  g =  g  vir .We then derive the total enclosed gas energy, corresponding to the product of the gas fraction and the total matter energy enclosed within the halo virial radius: where the sum runs over all grid cells  within the gas radius  g .
We then compute the initial gas surface pressure: where the sum runs over all   cells within a spherical shell with radii between  g and  g + 0.5 500 3.
3 We chose this value to guarantee that the shell thickness contains over ten times the simulations softening length even for the least massive halos.We confirmed that this parameter selection has negligible impact on the final results.

Polytropic rearrangement
The gas is then assumed to rearrange itself to follow a polytropic equation of state, in which the gas density and pressure in cell  are given by: where Γ is the gas polytropic index;  g,0 and  g,0 are the value of the density and pressure at the potential peak, and where  is the gravitational potential, and  0 its peak value.The gas injection is performed considering a fixed Γ (the value of which will be determined in §4).For a given distribution of gravitational potential, the gas density and pressure are thus fully determined by two parameters:  g,0 and  g,0 .To fix them, we postulate that the rearrangement of the gas conserves the total energy and the surface pressure.First, the gas radius  g is recomputed as the radius enclosing the initially-fixed total gas mass, by integrating the polytropic gas density defined in eq. ( 5).This allows the gas to expand or contract during rearrangement.Energy and surface pressure conservation are then enforced using this new gas radius.
Energy conservation -The rearrangement is set to conserve the initial gas energy  g, computed from eq. ( 3): The change in gas energy due to this contraction or expansion is accounted for through the Δ  term in the right-hand-side of eq. ( 7), computed as: The last term in the right-hand-side of eq. ( 7) is included to allow a fraction  DM of the initial dark matter energy to be transferred to the gas during the polytropic rearrangement.It was first introduced by Bode et al. (2009) as a way to account for energy transfer due to non-gravitational processes.The value of  DM is a free parameter of the model, and will be discussed in detail in §4.
Surface pressure conservation -We also enforce the conservation of the surface pressure of eq. ( 4) during rearrangement, i.e.
where the sum runs over a spherical shell with the same thickness as the one used in eq. ( 4), but located at the new gas radius  g .Fixing Γ and  DM , equations ( 7) and (9) constitute a system of two equations with two unknown parameters,  g,0 and  g,0 .We numerically solve for the two parameters by enforcing gas energy and surface pressure to be conserved to less than 0.1%.

MATCHING BARYON PASTING TO HYDRODYNAMIC HALO PROPERTIES
This section describes the methodology used to optimize baryon pasting to obtain gas properties matching those from hydrodynamic simulations.In particular, we focus on finding the numerical values of baryon pasting parameters that best reproduce radial profiles of gas density and pressure on an individual halo basis.

Model parameter space probed
As detailed in §3.2, for any given halo, numerically solving eqs.( 7) and ( 9) for  g,0 and  g,0 provides a complete description of the gas density and pressure on the grid used to project the total matter properties.This modeling is based on three physical parameters: the gas fraction contained within the virial radius  g ; the gas polytropic index Γ; and the fraction of dark matter energy transferred to the gas during rearrangement,  DM .We seek to find the optimal parameter combination for the model; to that end, we design a grid in the parameter space, and run a comparative study of gas properties for each node on the grid, for each independent redshift snapshot available.
We first choose to fix the value of the gas fraction  g to the cosmic baryon fraction,  g = Ω b /Ω m .This choice is motivated by both observation-and simulation-based studies concluding that the gas fraction enclosed within radius  converges towards the cosmic baryon fraction as  increases (and can thus be used as a cosmological probe, e.g., Mantz et al. 2022).With the presence of feedback, the radius at which this convergence is reached depends on halo mass, and might be well beyond the halo virial radius in the mass range considered here (see, e.g., Figure 2 in Ayromlou et al. 2022, for a comparative study over several simulations with different feedback prescriptions).Nonetheless, the Borg Cube is a non-radiative hydrodynamic simulation, without feedback or cooling; therefore, the gas fraction converges at smaller radii, with halo mass having little impact (see the lower panel of Figure 5 in E19).
For the polytropic index of the ICM we choose to study a range of Γ ∈ [1.13, 1.20].While many previous baryon pasting studies choose to fix Γ to a fiducial value of 1.2, hydrodynamic simulations and observations seem to indicate different values.In the case of the non-radiative run of the Borg Cube, values of Γ ∼ 1.17 are favored for cluster-scaled halos (see Figure 8 of E19), in agreement with analytical derivations for a self-similar, adiabatic collapse (Bertschinger 1985).Observationally, Γ can be estimated through the combination of ICM density and pressure measurements, which points towards slightly higher values (see e.g., Ghirardini et al. 2019, and references therein).Nonetheless, since we are comparing baryon-pasted gas to that of the Borg Cube, we choose a range close to the values measured in E19.Moreover, as we will show in §4.3, Γ values have a relatively small impact on the pasted properties.
Finally, an appropriate range of  DM is hard to predict a priori, as it is a specific parameter of this polytropicrearrangement gas modelling that has no direct measurable counterpart in observations or in hydrodynamic simulations.Previous studies have considered fixed values ranging between 0 (e.g., Osato & Nagai 2023) and 5% (e.g., Shaw et al. 2010).We thus choose to study that range, and will compute baryonpasted gas properties for  DM ∈ [0, 0.05].
It is worth emphasizing again at this point that some of the components of the model as introduced by Ostriker et al. (2005) and developed by e.g., Bode et al. (2009); Shaw et al. (2010) are not included here, such as the contribution of gas mass dropout towards stellar formation, and energy injection through feedback processes.This is due to our comparison of baryon pasting to non-radiative hydrdodynamics, which do not include these effects, preventing us from calibrating their modeling in baryon pasting.We address this point further in §6.We also choose not to focus on non-thermal pressure support, which will be the focus of a subsequent paper.

Gas properties comparison
With baryon-pasted gas properties computed for the desired set of parameters, we can compare them to the same properties in the hydrodynamic run of the Borg Cube.Since we used the same grid for the CIC projection of particle properties for both the GO and NR runs, we can perform the comparison at the cell level.
For each halo, and each set of model parameters, we compute the grids of relative differences, , as: From these grids, we compute radial profiles of the relative differences, ().The results at redshift  = 0 are shown in Figure 3.The impact of changing the two model parameters are clear.On one hand, an increase in Γ mostly reduces the central gas density and pressure (as can be seen by sliding from left to right in the panels of Figure 3).On the other hand, increasing  DM (i.e.sliding from top to bottom in Figure 3) increases gas concentration, creating more peaked profiles, and thus, slanted relative differences profiles.
Overall bias correction -The comparison of pasted gas properties with those from the hydrodynamic simulation reveals an overall multiplicative bias, which can be seen as a systematic offset in the relative differences profiles (e.g., in Figure 3).As an extension of the modeling presented above, we introduce two new parameters, Δ  and Δ  , which correspond respectively to an overall offset of the gas density and pressure.In this parametrization, the final pasted gas density and pressure are written as: where  g,pol.and  g,pol.are the polytropic gas density and pressure computed from eq. ( 5).The inclusion of these two parameters allows us to correct for the overall bias between baryon pasting and hydrodynamic simulations.For a given set of parameters (Γ,  DM ), the values of Δ  and Δ  are determined from the radial profile of relative differences   () and   ().First, the profiles of relative differences are computed for each halo by averaging   and   , defined in eq. ( 10), in concentric shells of width 0.1 ×  GO 500 .Then, we compute the median profiles of relative differences across all halos in the redshift snapshot.Finally, the bias values are computed as the average value of the resulting median profiles, on radii comprised within (0.25, 1.25) ×  GO 500 .The innermost radii are excluded as there is no set of parameters that provides a good description of the gas properties in the halo cores at low redshift (as can be seen in Figure 3).Moreover, this region poses inherent modelling challenges.One -Relative difference between gas properties from baryon pasting and from the Borg Cube NR run as a function of cluster-centric radius.We show the profiles at redshift  = 0 for gas density (blue) and pressure (orange).The solid line indicates the median over all halos, while the envelopes show the interval between the 16 th and 84 th percentile.Each subfigure shows the results for one set of baryon pasting parameters (Γ,  DM ), indicated at the top of the panel: each column corresponds to a fixed value of Γ (increasing from left to right), while each row shows results for a fixed value of  DM (increasing from top to bottom).concern is that very small scales are affected by the softening length of the simulation.In addition, clusters are most strongly affected by star formation and subgrid processes in their cores; the Borg Cube being a non-radiative simulation, this region is expected to be the least accurate representation of observed clusters.Therefore, we choose not to focus on cluster cores, and to save the reproduction of halo cores for a future study, that will include further modelling of baryonic effects and possibly feature higher resolution simulations (see discussion in §6).We also exclude the outskirts from the analysis because of their low particle counts: for low-mass halos, gas density and pressure profiles become very noisy beyond ∼ 1.5×  500 in the NR run.This is particularly problematic at high-redshift, where high-mass halos are not yet formed (see Figure 1).Therefore, in order for the comparison between baryon pasting and hydrodynamic gas to be valid for the full mass range of halos in this study, we choose to exclude radii beyond ∼ 1.5 ×  500 from the analysis.
The resulting values of Δ  and Δ  are then used to compute the final, unbiased gas density and pressure distributions for each halo through eq. ( 11).While the inclusion of these overall bias parameters enables achieving high-accuracy realizations of ICM gas from baryon pasting (as quantified in the next section), the origin of the discrepancy remains unclear.It could be partially explained by the chosen value for the injected gas fraction being too high, which would result in overestimated gas density and pressure.As discussed in §4.1, this effect is unlikely to be large enough to explain a 20% bias, as the gas fraction in clusters usually converges towards the cosmic baryon fraction at large radii.In the specific case of the nonradiative Borg Cube, the virial gas fraction of our halo sample is found to be at least 95% of the cosmic baryon fraction, corresponding to a possible global offset by 5% at most between baryon pasting and hydrodynamic densities.Another partial source of explanation could be that the fictitious polytropic rearrangement artificially pushes gas towards the cluster center in a way that cannot be captured by changes in (Γ,  DM ), which would point towards the need to modify the polytropic modelling of the gas.Finally, in extended versions of the Ostriker et al. (2005) model (e.g., Bode et al. 2007Bode et al. , 2009;;Shaw et al. 2010), a fixed fraction of the gas is allowed to cool down and be turned into stars.Although the Borg Cube hydrodynamic simulation does not include star formation, we note that this fraction would be degenerate with the bias parameters.3.

Best-fitting baryon pasting parameter selection
The performance of each parameter set is judged by the mean squared value of the relative pressure difference between pasted and hydrodynamic gas pressure.We choose to focus on pressure reconstruction as we are mostly interested in reproducing the tSZ signal, which is sourced by electron pressure in the ICM; this point is addressed further in §6.First, for each set of (Γ,  DM ), we compute the relative difference between the unbiased baryon pasting pressure (computed from eq. 11) and that of the gas from the NR run, using eq.( 10).For each halo, we then infer the corresponding radial profile in concentric shells of width 0.1 ×  GO 500 .We then compute the median profile across halos   (), and then the mean squared relative difference, where the average is again taken on the radii comprised within (0.25, 1.25) ×  GO 500 .Results -For each redshift snapshot, the set of (Γ,  DM ) that minimizes the value of  2  is then chosen as the set of parameters providing the best reproduction of gas properties using baryon pasting.Results are shown in Figure 4, with best-fit values reported in Table 3.We can identify a clear redshift trend for both Γ and  DM : from redshift  = 0 to 2, Γ increases from 1.15 to 1.18, and  DM from 0.5% to 3%.As for the overall biases on pasted gas density and pressure, we find that Δ  increases from ∼ 15% from  = 0 to ∼ 20% at  = 2, while Δ  is approximately constant over this range, with an average value of ∼ 20%.

QUALITY ASSESSMENT
We now focus on the quantitative assessment of the accuracy and precision of gas properties estimated via baryon pasting.While this assessment is often neglected in baryon pasting studies, it is necessary to quantify the reliability of the algorithm, and to understand how the procedure can be applied to cosmological studies.With this in mind, we measure the accuracy and precision of baryon pasting in its ability to reconstruct gas properties of cosmological interest, focusing on gas radial profiles and the scaling relation between the integrated Compton parameter,  500 , and halo mass,  500 .

Density and pressure profiles
The procedure described in §4 delivers an estimate of the baryon pasting parameters that give the best reproduction of gas pressure for each redshift snapshot.We now focus on the radial profiles of relative differences between baryon-pasted gas properties and their hydrodynamic counterparts for both density and pressure,   () and   (), defined in eq. ( 10). Figure 5 shows the corresponding radial profiles for every redshift snapshot.
To quantify bias, we first compute the median of   () and   () across all halos (shown as the solid blue and red lines in Figure 5).We then define the mean bias as the average absolute value of the resulting profiles across a fixed radial range.We again only consider the radial range  ∈ [0.25, 1.25]  GO 500 , excluding cluster cores because of their inherent inaccuracy in non-radiative hydrodynamic simulations.Over that range, we find that the average bias on density (pressure) is smaller than 3% (2%) for all redshift snapshots.
To judge the precision, we estimate the standard deviation of the profiles across halos.Over the radial range considered, we find that the density profiles are scattered, with an average standard deviation of the order of 15%.As for the gas pressure, we obtain higher values, with an average scatter ranging between 18% and 22% for the different redshift snapshots.
Note that part of this scatter is inherently due to the halo masses being slightly different in the GO and NR simulations.For example, in the case of a halo that has a higher mass in the NR run than in the GO run, we expect hydrodynamic gas density and pressure to be higher than their pasted counterparts.Since there is scatter in the masses of halos in both runs (as illustrated on the right panel of Figure 1), this directly implies scatter in the comparison between hydrdodynamic and pasted gas.

𝑌 500𝑐 |𝑀 500𝑐 scaling relation
A fundamental property of galaxy clusters that is highly important in cosmological analyses is the statistical relation between the integrated Compton parameter and cluster mass.The (spherically-) integrated Compton parameter,  500 , is defined as the three-dimensional integral of ICM electron pressure within  500 , i.e. -Relative difference between gas properties from baryon pasting and from non-radiative hydrodynamic measurements as a function of clustercentric radius.We show the profiles for gas density (blue) and pressure (orange).Each subfigure shows the results for a redshift snapshot, from  = 0 (top) to  = 2 (bottom), using the corresponding best set of baryon pasting parameters.The hatched grey regions correspond to the radial ranges not considered in the parameter search.
where  t is the Thomson interaction cross-section,  e  2 is the electron rest-frame energy, and  e is the thermal pressure of free electrons in the ICM. 500 is expected to scale as a power law of cluster mass (e.g., Kaiser 1991;Motl et al. 2005;Kravtsov et al. 2006), often written as (see e.g., Arnaud et al. 2010;Planck Collaboration et al. 2011 where Ỹ500 ≡  −2/3 () 500 , with  () the reduced Hubble parameter at the cluster redshift  ()/ 0 .Given the complexity of physical processes in cluster-scale halos, scatter is expected around this relation.It is therefore common to write it as the probability distribution followed by the integrated Compton parameter of a cluster given its mass, where  and  are the log-scaled integrated Compton parameter and mass, and  is the log-normally distributed intrinsic scatter around the scaling relation.
The spherically-integrated Compton parameter  500 as defined in eq. ( 13) is tightly linked to the integrated tSZ signal measurable in millimeter-wave cluster surveys (see Arnaud  et al. 2010;Hadzhiyska et al. 2023).This gives  500 practical use as a mass proxy, as knowing its scaling relation with mass enables its use to calibrate cluster masses from detections in a tSZ survey (which has been the basis of many cluster count cosmological constraints, as e.g., Planck Collaboration et al. 2016a).Moreover,  500 is known to be a relatively low-scatter mass proxy (i.e. with a tight correlation with mass and a low intrinsic scatter, Motl et al. 2005;Kravtsov et al. 2006), meaning it can be used to achieve precise mass estimates.Measuring the scaling relation between  500 and  500 has therefore been the driving science goal of many studies based on both observations (e.g., Arnaud et al. 2010;Planck Collaboration et al. 2011;Sereno et al. 2015;Kéruzoré et al. 2022) and simulations (e.g., Angulo et al. 2012;Wadekar et al. 2022;Pakmor et al. 2022;Hadzhiyska et al. 2023;Baxter et al. 2023).As our goal is to create datasets that provide accurate cluster properties for multi-wavelength cluster cosmology studies, ensuring the accuracy of the scaling relation between  500 and  500 is crucial.
We compute the value of  500 for each halo in our sample, using both the ICM pressure distribution in the NR run, and that reconstructed using baryon pasting.For each redshift snapshot, we then independently fit the scaling relation for both simulations using the BCES algorithm (Akritas & Bershady 1996), and use the generalization by Pratt et al. (2009) to estimate the intrinsic scatter around the power law.Results are shown in Figure 6.We note that the values we find for the slope of the relation  are a lot closer to the expected selfsimilar value of  = 5/3 than what is commonly measured in observations; for example, Arnaud et al. (2010) andPlanck Collaboration et al. (2011) respectively find slopes of  = 1.78 and 1.79.This is likely the consequence of the Borg Cube hydrodynamic simulation being non-radiative; feedback and cooling are known to impact the gas distribution in halos, in particular at low mass (e.g., Ayromlou et al. 2022), thus modifying the slope of the relation (as observed in, e.g., Yang et al. 2022).
Nonetheless, for each snapshot, we find closely-agreeing scaling relations between the hydrodynamic simulation and baryon pasting.Figure 7 shows the resulting scaling relation parameters, including confidence intervals estimated via bootstrap resampling -although they are very small at low redshift due to the sample size.We see that the reconstructed scaling relation parameters for pasted gas are consistently close to those obtained from the Borg Cube hydrodynamic simulation.We do note that the scatter around the relation  is slightly -but consistently -higher in the baryon pasting run, which could be due to cluster cores not being accurately reproduced by the model.Across all five redshift snapshots, we find this excess scatter to range between 1.5% and 4.9% of the value measured in the NR run, which is still far below typical uncertainties for current measurements of the  500 | 500 scaling relation (e.g., Arnaud et al. 2010;Planck Collaboration et al. 2011;Kéruzoré et al. 2022).

First look at tSZ projections
Finally, the major goal of baryon pasting is the creation of realistic sky maps from cosmological simulations; in particular, in the context of this work, the realization of accurate maps of the tSZ signal.The aim of this study was not to provide sky maps for the Borg Cube simulations, but to focus on exploiting the complementarity of the two runs to optimize a baryon pasting model that will be usable for other gravity-only simulations with larger volumes and better mass resolution.Nonetheless, for illustrative purposes, we provide example tSZ thumbnails of some of the halos in the Borg Cube simulation in Figure 8.The maps show the spatial distribution of the Compton− parameter, proportional to the line of sight-integrated ICM electron pressure (Sunyaev & Zeldovich 1972): where the line of sight (LoS) is taken along the  axis of the simulation.
For each redshift snapshot, we show maps of three halos covering the entire mass range of the sample: one of the most massive halos, one near the median mass, and one near the low-mass end.For each halo, we show a map obtained from the ICM pressure distribution in the NR run as well as from the pressure obtained from the baryon pasting model.We see that the baryon pasting model is able to reproduce the overall halo shape.On each map, we add a gray circle with a diameter corresponding to an on-sky angular size of one arc minute, roughly corresponding to the angular resolution of telescopes producing millimeter-wave sky surveys, such as the South Pole Telescope (SPT, Carlstrom et al. 2011;Bleem et al. 2020), the Atacama Cosmology Telescope (ACT, Swetz et al. 2011;Hilton et al. 2021), the Simons Observatory (Ade et al. 2019), and CMB-S4 (Abazajian et al. 2016).We see that only the most massive and low-redshift clusters are well resolved with an arc minute resolution, highlighting the fact that the small-scale discrepancies between baryon pasting and hydrodynamic simulations will be smoothed out in sky maps created to mimic current-and next-generation millimeter-wave cluster surveys.
6. DISCUSSION This work constitutes the first step in a larger effort to implement a baryon pasting pipeline to the suite of post-processing data products of HACC simulations.We begin this effort by focusing on reproducing the distribution of gas in non-radiative hydrodynamic simulations.Because more complete hydrodynamic simulations are complex and require some level of fine-tuning, the simpler non-radiative treatment of hydrodynamics provides a helpful baseline.In this context, ensuring that baryon pasting is able to reproduce the results of such important anchor points is an important step towards the building of a complete framework to populate gravity-only simulations with synthetic gas.Moreover, this effort evolves in parallel with recent developments of the CRK-HACC hydro solver (Frontiere et al. 2023), which include the implementation of complex subgrid modeling of baryonic physics.Future work will make use of results from these new simulations to extend the baryon pasting model, and ensure it can recreate gas properties similar to hydrodynamic simulations beyond the non-radiative treatment.
One of the main limitations of this study is the absence of complex baryonic processes such as star formation, cooling and feedback.This comes as a direct consequence of the dataset being used: since our goal is to reproduce the thermodynamic properties of intracluster gas in a hydrodynamic sim-ulation without subgrid physics modelling, our baryon pasting model cannot account for these effects.These effects are well known to play an important role in gas distribution, even at cluster scales (e.g., Battaglia et al. 2010;Donahue & Voit 2022), which limits the realism of gas properties obtained in our work.We emphasize again that this is to be addressed with future extensions of this work, using the same approach on full-hydrodynamic simulations, as well as information from observational data.In particular, we plan on exploring subgrid models simultaneously in baryon pasting and hydrodynamic simulations, which will allow us to explore a wide variety of models and to understand the impact of subgrid modeling uncertainties on baryon pasting and resulting maps.
We also note that our analysis has focused on accurately reproducing ICM gas pressure at the expense of gas density, meaning different results might have been achieved if we had focused on the density reproduction.As previously stated, this choice was made to prioritize the realistic reproduction of tSZ effect signal, as the amplitude of this effect is sourced by electron pressure in the ICM.A similar study oriented towards the optimal reconstruction of gas density might deliver bettersuited results for X-ray emission modelling (e.g., Lau et al. 2023), or of the kinetic Sunyaev-Zeldovich effect (e.g., Flender et al. 2016).Nonetheless, we find the accuracy and precision of density reconstruction to be on par with that of the pressure, as assessed in §5.1.
Finally, one might have concerns regarding the applicability of these results to other gravity-only simulations is the fact that both runs of the Borg Cube simulation comprised two species, with different particle masses.This makes the Borg Cube GO run atypical, as most gravity-only simulations consist in the evolution of a single species representing the combined fluid.The impact of this duality has been studied in detail in E19, showing evidence of a small-scale bias in the power spectra of the individual species due to artificial mass segregation.Nonetheless, the same study also demonstrates that the impact is negligible when considering the total matter density field.In particular, it was shown that the total matter field agrees with that of a standard single-species gravity-only run down to roughly twice the softening length.Hence, we do not expect baryon pasting to be significantly affected by the number of species used to run the input gravity-only simulation.

CONCLUSIONS
The realization of realistic sky maps from cosmological simulations is a crucial issue for observational cosmology.In particular, as millimeter wave-detected galaxy clusters surveys are poised to deliver high precision constraints on cosmological parameters, synthetic datasets provide a precious and necessary tool for the calibration of cosmological analyses.In this context, in parallel with efforts to create high-accuracy hydrodynamical simulations, it is useful to implement frameworks to take full advantage of gravity-only simulations, such as baryon pasting, which can be realized in post-processing at a fraction of the computational cost of complete simulations.
We have presented a study detailing the implementation of a baryon pasting algorithm and the investigation of its model parameters.In particular, we used the simple model introduced by Ostriker et al. (2005), and focused on comparing the gas thermodynamic properties produced by baryon pasting with those obtained from a non-radiative hydrodynamic simulation evolved from the same conditions.In doing so, we found the optimal model parameters to reproduce the pressure distribution in the ICM gas of cluster-scale halos, and assessed the performance achieved by the model when using these parameters.
Our main conclusions are as follow.
• Within the framework of the Ostriker et al. (2005) baryon pasting model, we find that the parameters achieving the best reproduction of pressure in non-radiative hydrodynamic simulations vary smoothly with redshift.From  = 0 to 2, the gas polytropic index Γ increases from 1.15 to 1.18, while the fraction of dark matter energy transferred to the gas during the polytropic rearrangement  DM increases from 0.5% to 3%.Values at intermediate redshifts are reported in Table 3 and illustrated in Figure 4.
• Using these sets of parameters and fixing the virial gas fraction to the cosmic baryon fraction, we find that baryon pasting produces systematically overestimated gas density and pressure compared to hydrodynamic simulations.We estimate these mean biases Δ  and Δ  , and find that the bias on density varies between ∼ 15% and ∼ 20% across the considered redshift range, while the pressure bias remains constant at ∼ 20%.Values are also reported in Table 3 and illustrated in Figure 4.As discussed in §4.2, the need for bias parameters may originate from a combination of an overestimated virial gas fraction, an incorrect assumption in the polytropic rearrangement at the heart of the pasting method, or due to the absence of star formation which would otherwise be captured in degenerate bias parameters.
• Correcting for these biases, we find that baryon pasting produces median gas pressure (density) profiles biased by less than 2% (3%) across all redshifts, on a radial range comprised between 0.25 and 1.25 500 compared to hydrodynamic simulations, with a scatter of ∼ 20% (∼ 15%).The bias and scatter of baryon pasting profiles compared to hydrodynamic ones is shown in Figure 5.
• We compare the scaling relation between the integrated Compton parameter  500 and halo mass  500 using baryon pasted-gas, and compare it to the scaling relation obtained from the same halos in the hydrodynamic run of the Borg Cube.We find consistent scaling relations for all redshifts, highlighting the usefulness of baryon pasting for calibrating the analysis of tSZ surveys.
The implementation of the model presented in this work marks the first step towards a systematic estimation of ICM gas properties in HACC simulations.This will result in the production of accurate millimeter-wave sky maps from the whole suite of HACC gravity-only simulations, adding yet another component to the synthetic datasets created from these products.In light of the importance of these simulations in the current cosmological landscape, this will open the door to many multi-wavelength cosmological studies to prepare the exploitation of next-generation large-scale surveys.

Fig. 1 .
Fig. 1.-Distribution of halo properties for the Borg Cube halos used in this work.From left to right, we show the distribution of (a) GO halo masses  500 and (b) concentrations  500 , (c) the offset distance between the FOF centers of matched halos in the GO and NR runs, and (d) their relative mass difference.Each color represents the distribution for one of the redshift snapshots.
Fig. 2.-Left: Comparison between halo distributions in a (100 × 50 × 800) (ℎ −1 cMpc) 3 subvolume of the GO (top) and NR (bottom) runs of the Borg Cube at  = 0. Right: Zoomed-in surface density maps of the different matter components in the region highlighted with an orange square on the left panel, for the GO (top) and NR (bottom) runs.In every panel, each open circle marks the position in the ( , ) plane of a halo with  GO 500 ⩾ 10 13.5 ℎ −1  ⊙ , and its radius is the halo characteristic radius  500 .Only halos matched in both runs are shown.The same structures appear in the halo distribution and the total matter distribution of both runs, showing the relevance of comparing the properties of matched counterparts.Surface density maps were computed from particle data using SDTFE (Rangel et al. 2016).Redshift  GO halos  both halos  = 0 20458 20120  = 0.5 11880 11703  = 1 4698 4644  = 1.5 1315 1311  = 2 263 260 Total 38614 38038 Fig.3.-Relativedifference between gas properties from baryon pasting and from the Borg Cube NR run as a function of cluster-centric radius.We show the profiles at redshift  = 0 for gas density (blue) and pressure (orange).The solid line indicates the median over all halos, while the envelopes show the interval between the 16 th and 84 th percentile.Each subfigure shows the results for one set of baryon pasting parameters (Γ,  DM ), indicated at the top of the panel: each column corresponds to a fixed value of Γ (increasing from left to right), while each row shows results for a fixed value of  DM (increasing from top to bottom).
Fig. 4.-Redshift evolution of the baryon pasting parameters Γ,  DM , and biases on density Δ  and pressure Δ  which achieve the best results in reproducing the ICM pressure distribution found in the NR run of the Borg Cube.Numerical values are reported in Table3.
Fig.5.-Relative difference between gas properties from baryon pasting and from non-radiative hydrodynamic measurements as a function of clustercentric radius.We show the profiles for gas density (blue) and pressure (orange).Each subfigure shows the results for a redshift snapshot, from  = 0 (top) to  = 2 (bottom), using the corresponding best set of baryon pasting parameters.The hatched grey regions correspond to the radial ranges not considered in the parameter search.

Fig. 6 .
Fig. 6.- 500 |  500 scaling relation obtained on the full halo samples, using baryon pasting (BP, left) and the NR run (right).Each row represents a redshift snapshot, indicated on the right of the figure.In each panel, the orange line shows the BCES estimator of the power-law relation between Ỹ500 and  500 of eq.(14).The corresponding values of parameters ( , , ) are reported in the lower right corner of each panel.The grey hatched region shows the mass cut on gravity-only  500 used to select the sample.

Fig. 7 .
Fig. 7.-Best-fit parameters of the  500 |  500 scaling relation from eq. (14).Each color represent a redshift snapshot, and shows the results using baryon pasting (filled circles, BP) and the non-radiative run (open circles, NR).Error bars are computed by bootstrap resampling of the sample, and show the 16 th and 84 th percentile of the resulting parameter distribution.The dashed vertical black line shows the self-similar expectation of  = 5/3.