A Hydro-Particle-Mesh Code for Efficient and Rapid Simulations of the Intracluster Medium

We introduce the cosmological HYPER code based on an innovative hydro-particle-mesh (HPM) algorithm for efficient and rapid simulations of gas and dark matter. For the HPM algorithm, we update the approach of Gnedin&Hui (1998) to expand the scope of its application from the lower-density intergalactic medium (IGM) to the higher-density intracluster medium (ICM). While the original algorithm tracks only one effective particle species, the updated version separately tracks the gas and dark matter particles as they do not exactly trace each other on small scales. For the approximate hydrodynamics solver, the pressure term in the gas equations of motion is calculated using robust physical models. In particular, we use a dark matter halo model, ICM pressure profile, and IGM temperature-density relation, all of which can be systematically varied for parameter-space studies. We show that the HYPER simulation results are in good agreement with the halo model expectations for the density, temperature, and pressure radial profiles. Simulated galaxy cluster scaling relations for Sunyaev-Zel'dovich (SZ) and X-ray observables are also in good agreement with mean predictions, with scatter comparable to that found in hydrodynamic simulations. HYPER also produces lightcone catalogs of dark matter halos and full-sky tomographic maps of the lensing convergence, SZ effect, and X-ray emission. These simulation products are useful for testing data analysis pipelines, generating training data for machine learning, understanding selection and systematic effects, and for interpreting astrophysical and cosmological constraints.


Introduction
Cosmological simulations have greatly improved our understanding of the formation and evolution of the cosmic structures throughout our universe and are widely used to interpret observations and design new instruments and surveys. N-body simulations of gravitational dynamics make detailed and reliable predictions for the distributions of dark matter, which forms the backbone of structure formation, and dark energy responsible for the accelerated expansion of the universe (e.g., Springel et al. 2005;Boylan-Kolchin et al. 2009;Klypin et al. 2011;Angulo et al. 2012;Skillman et al. 2014;Klypin et al. 2016;Ishiyama et al. 2021). Hydrodynamic simulations, which model the coupled evolution of dark matter and cosmic gas (e.g., Iannuzzi & Dolag 2012;Vogelsberger et al. 2014;Khandai et al. 2015;Schaye et al. 2015;Dolag et al. 2016;Dubois et al. 2016;Feng et al. 2016b; Barnes et al. 2017;Kaviraj et al. 2017;McCarthy et al. 2017;Tremmel et al. 2017;Pillepich et al. 2018a;Davé et al. 2019), are also able to predict many directly observable properties of cosmic gas and galaxies.
Understanding the large-scale structure of our universe requires two parts: (i) an accurate solution to the equations of motion for the dark matter and (ii) physically reasonable approximations for the behavior of baryonic components of the universe. Together, these two components build the foundation of the modern observational cosmology. Dark matter N-body simulations have achieved significant progress in developing the understanding of the structure of dark matter halos (e.g., Navarro et al. 1996Navarro et al. , 1997Navarro et al. , 2004 and their clustering (e.g., Springel et al. 2006;Tinker et al. 2008;Watson et al. 2013;Bocquet et al. 2016). They were instrumental in establishing the Λ cold dark matter (CDM) cosmological model as the dominant paradigm for the nature of both dark matter and dark energy, but suffer from the fundamental limitation of being incapable of providing any direct prediction for the baryonic component. Thus, the hydrodynamic simulations that also simulate the baryons that form the visible components in our universe are crucial for interpretation and calibration of the results of observational data. A detailed review of the modern hydrodynamic simulations studying the properties, growth, and evolution of galaxies is given in Vogelsberger et al. (2020).
Hydrodynamic simulations are generally preferred for solving the nonlinear physics of structure formation and predicting the survey observable dependence on cosmological parameters. However, the effective volumes of modern surveys keep growing, and achieving the science goals of these surveys requires numerical simulations of exceptionally large volumes -both for correctly capturing the statistics of the rare objects and for computing the covariance matrices between the observables. Simulations in spatial volumes comparable to the surveys in size are generally too expensive to make many large-scale mock observations and explore both astrophysical and cosmological parameter space. In the face of increasing demand for multiple realizations of simulated mock catalogs for comparison with the large-scale structure observations, fast approximate approaches for dark matter simulations based on semi-numerical methods and Lagrangian perturbation theory have been developed. For example, PTHALOS (Scoccimarro & Sheth 2002) has been used for efficiently generating mock galaxy distributions, PINOCCHIO (Monaco et al. 2002) is capable of accurately predicting formation and evolution of individual dark matter halos, and COLA (Tassev et al. 2013) and FastPM (Feng et al. 2016a) can be used for cheaply generating large ensembles of accurate mocks that properly account for nonlinear evolution. These fast approximate methods have shown their ability to reduce computational complexity and required computational resources by orders of magnitude without sacrificing accuracy on large scales.
Though we have seen significant progress in various approaches aiming to speed up dark-matter-only N-body simulations, there is still a notable lack of fast approximate hydro-simulation methods. Previously, Gnedin & Hui (1998) used the particle-mesh (PM) solver for dark matter dynamics and allowed for the additional gas pressure force to approximate hydrodynamics. Their hydro-particle-mesh (HPM) algorithm substantially relies on the existence of a tight temperate-density relation in the intergalactic medium (IGM) and has been successfully used to model the highredshift Lyα forest with moderate precision (McDonald et al. 2002). However, the tight correlation between the gas density and temperature in the low-density IGM breaks down in denser regions. Yet, it is possible to extend the range of validity of HPM-like techniques further. For example, in order to model the intracluster medium (ICM) of galaxy clusters, we can adopt empirical or simulated ICM pressure profiles (e.g., Arnaud et al. 2010;Battaglia et al. 2012a;He et al. 2021) and build a mapping relation between the gas temperature or pressure and some properties of cosmic gas that can be captured by, say, a simple PM solver. Such an approach will allow us to implement a fast approximate method for modeling hydrodynamics in the high-density ICM; the gas physics can then be modeled very efficiently in both the IGM and the ICM regime, which together fill most of the spatial volume in a fast hydrosimulation.
This paper introduces an innovative HPM code for efficient and rapid (HYPER) simulations of gas and dark matter. HPM simulations take approximately two to three times as long to run as PM simulations and are orders of magnitude faster than expensive hydrodynamic simulations. HYPER allows one to systematically vary the ICM and IGM models to study different baryonic physics and effects. We organize this paper as follows: Section 2 discusses the model for radial profiles of dark matter and gas in the ICM. Section 3 briefly reviews the HPM algorithm in the IGM, then describes the implementation of the HPM for the ICM regime, including how to modify the PM code to calculate the designed HPM fields and infer the gas temperature and pressure from the local field information using a pre-computed mapping derived from a given ICM model. In Section 4 we evaluate our new fast hydro-simulation performance by comparing the radial profiles, the integrated halo quantities, and the statistical quantities of the tSZ effect in our simulation to the predictions of the ICM model we use to implement the HPM algorithm. In Section 5 we conclude with our findings for the output of an HYPER simulation, and at the end, we also bring out some perspective of use cases for HYPER and the future extensions of this work.

Models
This section introduces our approach for constructing an analytical model for the radial profiles of different components in galaxy clusters, under the assumption of hydrostatic equilibrium. It is based on thoroughly studied modeling work on dark matter distribution and gas pressure support. We also briefly introduce the analytical relation used by the HPM algorithm (Gnedin & Hui 1998) to simulate gas thermal properties in the IGM.

Halo Model
As a model for the mass distribution in a dark matter halo of a galaxy cluster, we adopt a universal Navarro-Frenk-White (NFW) density profile (Navarro et al. 1997) ( ) (  where R 500c defines the radius of a spherical over-density region with average density 500 times of the critical density ρ crit (z). The total mass of the over-density region can be used to define the halo mass as 3 . The halo concentration c 500c has been calibrated as a function of cluster mass and redshift (e.g., Diemer & Kravtsov 2015). Thus, given the halo mass and redshift, we can specify the halo density profile of dark matter. The total mass enclosed within the radius r then has an analytical expression We use the NFW profile to approximate the total matter distribution, considering that dark matter dominates the mass contribution to the halo. The error on the total mass profile is mainly caused by the difference between the gas density profile and the NFW profile, which is generally considered to be a goof approximation to the dark matter distribution of halos. Though the deviation of the normalized gas density profile from the NFW model may not be negligible, the difference between the total matter density of halos and the NFW profile is much less significant since gas is subdominant in the total matter distribution. The difference between the total matter density profile and the NFW model decreases by a factor of the baryonic fraction f b ∼ 0.15 compared to the gas density profile and should be much less than unity. Thus teams with a lot of analytical work choose to model the total matter density profile with the NFW model, so that ρ m (r) = ρ dm (r) + ρ gas (r) ; ρ NFW (r) (e.g., Olamaie et al. 2012;Shi & Komatsu 2014). X-ray and lensing studies have also shown that NFW profiles can generally provide adequate descriptions of the mass distribution of cluster halos (e.g., Mandelbaum et al. 2006;Schmidt & Allen 2007;Morandi & Limousin 2012).

ICM Model
We construct an analytical model of the ICM based on the assumption of hydrostatic equilibrium and models for describing the gas pressure of galaxy clusters. The gas thermal pressure profile has been widely studied by both hydrodynamic simulations and SZ and X-ray observations. In most of these studies, the thermal pressure of gas generally follows an analytical form where P 500c is the pressure scale derived from the self-similar scaling relation where E(z) = H(z)/H 0 and ( ) = x r R 500c is a generalized NFW (GNFW) model (e.g., Nagai et al. 2007), Parameters P 0 , c 500 , α, β, and γ are normally fitted to the observed or simulated scaled pressure profiles. The ICM gas pressure profile modeled with only the GNFW term assumes a simple self-similar behavior (e.g., Kaiser 1986;Voit 2005). However, several numerical and observational works (e.g., Arnaud et al. 2010;Gupta et al. 2017) claimed to observe a deviation from the self-similar scaling relation in the gas pressure, which we denote by the additional terms f (M), g (z) in Equation (5). We adopt power-law forms for these terms, where M å is a chosen constant referred to as the characteristic cluster mass. However, the dependence on halo mass and redshift of the gas pressure model needs further study, as the studies on these deviation terms with numerical simulations and observations have not reached a final agreement.
In this work, we adopt the debiased pressure profile (DPP) from a recent study of the gas pressure model of the galaxy clusters (He et al. 2021), which adjusted the universal pressure profile (Arnaud et al. 2010) for hydrostatic mass bias in X-ray observation. The DPP model strictly follows the GNFW model and omits the terms describing the deviation from f (M) and g(z) for the gas thermal pressure P th (r). The GNFW parameters P 0 , c 500 , α, β, and γ of the pressure model are [5.048,1.217, 1.192,5.490,0.433].
Further studies of the cluster outskirts with the latest cosmological simulations (e.g., Shaw et al. 2010;Nelson et al. 2012;Battaglia et al. 2012b;Lau et al. 2013;Nelson et al. 2014;Gupta et al. 2017) and the observations of galaxy clusters (e.g., Bautz et al. 2009;George et al. 2009;Reiprich et al. 2009;Hoshino et al. 2010;Kawaharada et al. 2010;Simionescu et al. 2011;Urban et al. 2011) showed existence of the nonthermal pressure support, which is non-negligible particularly in the outskirts of galaxy clusters. The nonthermal pressure is mainly due to the nonthermal gas processes like virialized bulk motions and turbulent gas flows, which are primarily generated by mergers and accretion during the cluster formation. Then the total pressure of gas P tot (r) should consist of both thermal and nonthermal contributions, modeled by where f th (r) and f nth (r) are ratios of the thermal and nonthermal terms to the total pressure. Studying the sample of 65 galaxy clusters from a highresolution hydrodynamic cosmological simulation, Nelson et al. (2014) characterized the nonthermal pressure fraction profile f nth (r), and found it was universal across redshift and cluster mass of the studied 65 galaxy clusters. When scaling the cluster radii with respect to the mean matter density of the universe, the thermal fraction we mentioned above can be then expressed as with the best-fitted parameters A = 0.452, B = 0.841, and γ = 1.628. Then based on the assumption of the cluster dynamical state being in hydrostatic equilibrium, we are able to predict the gas density profile from analytic models for the enclosed mass function of halo M(r), the gas thermal pressure profile P th (r), and the gas thermal fraction f th (r) as The gas temperature profile of the cluster is derived by the ideal gas equation of state, where μ is the mean mass per gas particle, and k B is the Boltzmann constant. We emphasize that the analytical models we choose for both the thermal and nonthermal pressure of the gas in this paper are only presented as an example, not some hard-coded choice. In HYPER, the hydro part that drives the gas particle hydrodynamics can be constructed based on general models for the gas pressure. Thus we can systematically vary the pressure profiles and study different ICM models with HYPER.

IGM Model
In the low-density IGM regime, where the local density is less than the mean density of the universe that  r r 10 m m , shock-heating is not important, and the gas physics is set by different processes than in the ICM regime in Gnedin & Hui (1998). Hui & Gnedin (1997) proposed a semi-analytical method to predict the temperature-density relation for any given cosmology and reionization history, using the Zel'dovich approximation. They find a tight correlation (to better than 10%) between the gas density and temperature (and hence pressure as well), well described by a power law, where T 0 (z) is a function of redshift only and is of the order of 10 4 K, Δ is the relative gas density,r r D = gas gas , and γ is between 1 and 1.62. Both T 0 and γ are found evolving with time. Hui & Gnedin (1997) derived analytical approximations to the temperaturedensity relation in a scenario when the universe reionizes rapidly, which is often considered a suitable approximation for the low-redshift evolution of γ and T 0 .

Methods
In this section we introduce our further development of the HPM algorithm, which has already shown its potential for being an efficient and accurate alternative to full hydrodynamic simulations for simulating the low column density Lyα forest (Gnedin & Hui 1998), by expanding it from the IGM regime to the high-density ICM with the help of the analytical ICM model discussed in Section 2.

Particle-mesh
The PM solver, which laid the foundation of the HPM algorithm, is widely used for tracking the evolution of collisionless dark matter (e.g., Hockney & Eastwood 1981). We start with introducing a PM code developed by Trac & Pen (2004), which demonstrates how we resolve the evolution of collisionless dark matter particles in our HPM algorithm. This PM code has been adapted for solving Newton's equations of motion for dark matter particles in a Friedman-Robertson-Walker (FRW) universe. In the expanding FRW background, we use co-moving coordinates x c = x/a, where the scale factor a is governed by the Friedman equation assuming a spatially flat background. And to preserve the timeinvariant conservation form of the Euler fluid equations, we take a new variable for the time coordinate (e.g., Doroshkevich et al. 1980;Gnedin 1995;Martel & Shapiro 1998;Pen 1998;Trac & Pen 2004). The time-dependence of the cosmological expansion remains only in the gravitational source term, Here the matter density is per co-moving volume d 3 x c = a −3 d 3 x, and is related to the proper mass density by r r = a m 3 m p . The dynamical equations of dark matter particles under the spatial and time coordinate transformations are where code quantities velocity v and gravitational potential Φ are related to the corresponding physical proper quantities by v = av p , Φ = a 2 Φ p . The PM solver solves the equation of motion for the dark matter particles treating the gravitational force as a field and approximating it on a uniform mesh. We first use the "Cloudin-a-Cell" scheme (Hockney & Eastwood 1981) to deposit mass, including both the dark matter and gas, to the mesh to create the density field, then calculate the gravitational potential Φ by solving Poisson's equation, with fast Fourier transform (FFT). In Fourier space, the modified Poisson's equation of the continuous system is expressed as For the PM solver in our HPM algorithm, the Green function is obtained from directly transforming the finite-difference approximation of the Laplacian in the Poisson's equation, where Δl = L/N is the length of unit grid cell, L is the length of simulation box, and N is the number of mesh cells per side. When calculating the force field, differential operators, such as the gradient ∇, are replaced by the finite-difference approximations. Potential and accelerations at particle positions are obtained by interpolating on the array of mesh-defined values. We adopt the kick-drift-kick leapfrog integrator also used for GADGET-2 (Springel 2005) in our simulation, updating the position and velocity for each particle with the time evolution operator ( ) . The kick and drift operators, K and D, are defined as where the acceleration a i is interpolated for each particle from the acceleration field on the grid.

Hydro-particle-mesh
In a hydro-simulation, the main difference between dark matter and gas, dynamically, is that the latter is subject to pressure forces in addition to the gravity. Our HPM method is designed by modifying the dark-matter-only PM algorithm to also solve the dynamical equation for gas particles, gas gas which, in addition to the gravitational force, also includes the pressure force − ∇P/ρ gas . The co-moving gas pressure P is related to the proper pressure P p by P = a 5 P p . The pressure force on gas particles is inferred from HPM variables calculated and saved for each gas particle interpolated to the particle positions on HPM fields and a mapping relation derived from the gas model of ICM and IGM we have already discussed in Sections 2.2 and 2.3. We will show how to select the HPM fields in Section 3.3 and how to construct the mapping relationship in detail in Section 3.4. We calculate the gravitational force − ∇Φ following the same procedure as for the original PM code described above. The pressure force − ∇P/ρ gas is calculated by depositing the particle pressure to the mesh, applying the finite difference, and finally interpolating it back to the particles in the simulation. Specific details of the implementation of the gas pressure calculation with the constructed mapping relation between the gas particle thermal properties and the HPM variables are discussed in Section 3.4.
Notice that we are no longer adopting a single component model to reduce computational cost as in the original HPM. As explained by Gnedin & Hui (1998), in the IGM regime, the gas pressure is dynamically subdominant on large scales; hence, if one is only interested in modeling the baryonic component, one can treat the part of the gravitational force acting on baryons from the dark matter as if the dark matter followed the baryons. However, in the high-density ICM regime, pressure and gravity are comparable quantities, dark matter distribution can not be accurately approximated by the distribution of baryons, and the difference between the two components can not be neglected when solving for the evolution of baryons. Hence, this difference between the IGM and the ICM requires one to track gas and dark matter particles separately.
The results in this paper are drawn from an HPM simulation of a box whose per-side length L = 500h −1 Mpc, with periodic boundary conditions and equal numbers of dark matter and gas particles N dm = N gas = 1024 3 , and the mesh size is set to be N mesh = 4096 3 . We show the visualization of simulation results at redshift z = 0.0 in Figure 1. This is only a test simulation, and we want to eventually run Gpc h −1 sized boxes for production simulations to capture the abundance of clusters (and large groups) and resolve them with sufficient mass and spatial resolutions. Our choice of ratio of particle number to simulation box size guarantees that there are enough simulated particles (N  10 4 ) for the composition of each massive halo (M  10 14 M e ) in the simulation. In standard PM simulations, the ratio of the number of particles to the number of meshes per side length is normally = N N : 1 :2 p 1D mesh

1D
; but in this way, we can only capture halos with at least 1000 particles and do not resolve them. Note that one particle in a mesh cell corresponds to an over-density of 8. In this paper, we choose = N N : for the test simulation so we can capture halos with about 100 particles. Ideally, we want about = N N : for the future production simulations cause then one particle per cell corresponds to an over-density of 512, which is the internal density of clusters. Our choice of = N N : in HYPER also makes the unit length of mesh grid small enough to suppress the region that suffers the numerical effect due to the limited resolution of simulation in the center of the simulated halos. Our choice ensures that the numerical effect is negligible beyond the core region (r  0.1R) of massive halos (M  10 14 M e ), where M and R are the mass and radius of simulated halos. We only want to resolve the region of r > 0.1R because both observed and simulated cluster profiles are uncertain below this scale.
The HYPER code is OpenMP parallelized and allows us to scale the code over multiple cores in a single node. The simulation runs on a single node consisting of four Intel Xeon Platinum 8260M 'Cascade Lake'' CPUs, which contain 96 cores in total, from start point z = 100 to z = 0 while the hydro part for gas particle is turned on at z = 6 takes ∼4500 CPU hours in total. Though we only adopt OpenMP for parallelization in the current development implementation, there is nothing in the method that prevents it from being implemented with Message Passing Interface (MPI) in a straightforward manner for future production simulations. In (Left) Shown is a projection of the dark matter density field. The cosmic large-scale structures like massive collapsed halos, elongated filaments, and near-empty voids in the HYPER simulation can be easily identified in the thin slice. (Right) The projection of the gas density field is shown to have a strong correlation with the dark matter distribution. In this slice, brightness indicates the projected mass density, and color hue visualizes the mean projected temperature (dull-red to brilliant-yellow indicating cold to hot, as shown by the color bar aside). the simulation, we adopt a concordance ΛCDM model (Ω m = 0.3, Ω Λ = 0.7, Ω b = 0.045, h = 0.7, and σ 8 = 0.8).

HPM Variables and Fields
For the ICM model we have discussed in Section 2.2, if the halo mass M and the distance of the gas to the halo center r are specified, one can derive the gas thermal properties from these two variables, where X refers to the gas thermal quantities like temperature T or pressure P. However, in the HPM algorithm, the halo mass and the displacement of gas particles with respective to the halo center are not directly available (without adding on-the-fly halo finding), and evaluating the gas properties in the ICM region after specifying the halo it resides in will also enforce unwanted spherical symmetry. Our goal is to avoid mapping the gas temperature or pressure directly from the halo mass and the gas particle displacement in the halo. Instead of identifying the halo where gas resides and locating the gas particle with respective to the halo center, we choose to use designed HPM variables and fields that could more accurately reflect the local environment information of simulated particles. The HPM variables interpolated to the position of gas particles with the HPM fields should also show enough connection to the halo mass and radius when assuming an ideal spherical symmetry scheme, so that we can use them to build a pre-computed mapping for inference of the gas thermal properties based on the ICM model. Since it requires two variables M and r to calculate the thermal quantities of the gas in the ICM model, the number of HPM fields used for calculating the HPM variables as alternative to halo variables M and r should also be at least two to break the degeneracy. According to the description above, the HPM fields chosen should satisfy two conditions to keep our HPM algorithm efficient enough: (a) the calculation of the HPM fields in the HPM method should be efficient enough, and (b) derivation of HPM variables in the adopted ICM model should be simple, preferably with analytical expressions as a function of halo mass and radius.
In HYPER, we adopt the matter density and a new field, the scalar force, to be our HPM variables for the inference of gas thermal properties. We choose the matter density because this quantity is readily available, since it is used to solve the Poisson's equation and is saved for each particle in the simulation. For the other HPM variable, the scalar force, we create a new variable based on the idea that originates from the gravitational force calculation. Recall that the gravitational force (per unit mass) is defined by where G is the gravitational constant, ⊗denotes the operation of convolution. For the newly designed scalar force (per unit mass) variable, we simply replace | | x x 3 in convolution with | | x 1 2 . The new variable shares the same units as the gravitational force, but is a scalar instead, The scalar force f scalar satisfies the requirements for the HPM variables. Thus we can quickly implement it in the PM solver, since the Fourier transform of | | which can be solved directly analogous to computation of the gravitational potential with the PM solver described in Section 3.2. And in the halo model, we can calculate the radial profile of the scalar force by using the NFW profile for the matter density of the halo convolved with | | x 1 2 . For the scalar force of a halo at radius |x| = r, we can can rewrite the integral in the spherical coordinates, so that Thus, for the given halo mass and radius, we can calculate the scalar force with Equation (29), which will facilitate our construction of the mapping relation between the HPM variables and the gas thermal properties in Section 3.4. Other candidates for the HPM variables that have been considered are gas density or gravitational acceleration, since they are also computed by HPM code. However, if the gas density is one of the fields determining the gas temperature and pressure, then there is no way to prevent artificial numerical fragmentation when the local Jeans' length becomes too small (Truelove et al. 1997) with the fixed spatial resolution of the uniform PM grid. With using the total mass density as the HPM variable, such numerical artifacts are greatly suppressed.
For another candidate, the gravitational acceleration mentioned above is similar to the matter density since dark matter also dominates the source of gravitational force. Thus, this quantity is also very stable against numerical fragmentation. However, the gravitational acceleration is significantly affected by finite mesh resolution and is underestimated in the center regions of simulated halos. The finite resolution effect in the simulation results in a non-monotonic profile of the gravitational acceleration. The non-monotonicity leads to multiple solutions if the gravitational acceleration is used to predict the thermal properties of gas particles in simulated halos, while we find that the scalar force is not subject to the non-monotonicity problem. We have also considered using the gravitational potential as one of the HPM variables, whose radial profiles of simulated halos are also monotonic. However, the gravitational potential suffers from the dynamic range being too narrow, and it has large-scale contributions that bias the local environment. After experimentation, we found that using the gravitational potential is less optimal than the scalar force.

HPM Table
With the HPM variables discussed in Section 3.3, matter density and scalar force, we aim to construct a mapping relation between the gas thermal properties and these designed halo variables. This mapping relation plays a crucial role in efficiently modeling the properties of gas in the HPM method for both the IGM and the ICM, and is referred to as the HPM table.
Construction of the HPM table in the low-density IGM regions where  r r 10 m m and in the high-density ICM region where¯ r r 10 m m in the simulation follows very different rules due to the dissimilar behavior of the gas in these two regimes. Our choice of the over-density threshold is based on the results in the previous study of the HPM simulation (Gnedin & Hui 1998). It shows that in the low-density IGM regime, the gas thermal properties can be very easily characterized through a power-law temperature-density relation mentioned in Section 2.3 with just one of the HPM variables, the matter density ρ m . They find the average error of gas density distribution generally stays within 5% in the HPM simulations for  r r 10 m m when compared with the gas density fields in full hydrodynamic simulations. In the experiments for the HYPER simulation, we adopt a different over-density threshold   r r 5 20 m m since the over-density at far edge of the halos r ∼ 1.5 − 2R is around¯ r r 20 m m . We do not observe a significant effect on the simulation results of the ICM regime. In the high-density regime, the tight correlation between the gas density and the gas pressure breaks down, and we are no longer able to describe the densitytemperature relation by a power law. However, with the help of the ICM model, which introduces a mapping relation between the gas thermal properties and halo information, M and r, we can construct a bridge between the HPM variables and the halo information and a mapping relation from the halo variables to the gas thermal properties.
We construct a probabilistic relation between the gas thermal properties, like temperature or pressure, and the designed HPM variables using a Bayesian analysis. Given the HPM variables for a gas particle ρ m and f scalar , its average temperature or pressure can be estimated using a Monte Carlo approach: and X(M, r) is the gas pressure or temperature defined by the ICM model discussed in Section 2.2.
In the continuous limit, this mapping is deterministic, since p (ρ, f|M, r) should be modeled as where ρ m (M, r) and f scalar (M, r) are obtained from the ICM model with Equations (1) and (29). In HYPER, we consider a discretized sampling method for our temperature estimation scheme, where we use a top-hat function to approximate p(ρ m , f scalar |M, r) where H defines the width of the top-hat function. We note that in the limit of a vanishing width of the top hat, we recover the aforementioned deterministic relation. For the prior distribution for M, r, we decompose the distribution as where the prior for mass p(M) and the conditional p(r|M) are given as Here dn dM is the halo mass function, and ρ NFW (M, r) is the halo density function defined by the NFW model. Since we can not directly sample from p(M, r) without an on-the-fly halo finder, we apply the importance sampling, which is a general technique in statistics for estimating properties of a particular distribution while only having samples generated from a different distribution than the distribution of interest, for estimation of temperature or pressure, where M, r are sampled from a uniform distribution on a logarithmic scale. We denote this covering distribution as p S (M, r). We can then re-weight the sampled points as As mentioned above, we infer temperature and other thermal properties of the gas particles in the simulation from a mapping relation between HPM variables and the thermal properties of gas particles we are interested in by

scalar
For the intermediate region between the IGM and the ICM model, we adopt a log-linear interpolation along the density axis of the HPM table.
In Figure 2 we plot the HPM tables, which show the mapping relation between the HPM variables, the matter density ρ m and the scalar force f scalar , and the gas temperature and pressure. In both panels of the plot, we see a vertical band along the axis of the scalar force atr r~10 m m that separates the HPM table into the IGM and ICM parts. We find the HPM table for the gas pressure generally follows a trend of larger pressure for particles with higher density and larger scalar force. However, the mapping relation for the gas temperature has a more complicated pattern in the top-right region of the HPM table, where particles with high density and large scalar force are found. These particles are more likely to reside in the halo centers, where the temperature profile is not necessarily monotonic.

Smoothing and Clumping Effects
In the HPM simulation, due to the finite resolution of the mesh and the finite number of particles, the simulation results in the highest-density regions like the inner cores of simulated halos suffer from a smoothing effect, which leads to underestimating output quantities of the simulation, like the density and the scalar force. In addition, few of the simulated halos are completely relaxed (especially in the outskirts) or perfectly spherical, which violates the spherical symmetry assumption in the ICM model. We refer to the asymmetry and inhomogeneity in the distribution of particles in the outskirt region of simulated halos as the clumping effect. Both the smoothing effect and the clumping effect will cause the deviation of the simulated profiles for the HPM variables from the ICM model prediction. The smoothing effect will lead to an underestimation of simulated HPM variables calculated for the gas particles in the halo inner region, while the clumping effect in the outer region will lead to an overestimation.
In Figure 3 we show the comparison between the radial profiles of HPM variables, the matter density, and the scalar force, and their ICM model predictions calculated for gas particles of simulated halos within the halo mass bin centered at 3 × 10 14 M e at redshift z = 0. In both panels, the simulation results tend to underestimate the HPM variables from the ICM model in the inner-core region due to the smoothing effect of the finite numerical resolution of the simulation. Outside the core, the volume-weighted matter density profile is in good agreement with the NFW profile. We also find that the massweighted profiles are significantly higher than the volumeweighted ones, which indicates that the clumping effect in the outer regions of the simulated halos is non-negligible.
Notice that since we infer the gas pressure or temperature from the HPM variables calculated for each particle, the offset between the simulation results and the ICM model would introduce bias to the inference of gas temperature or pressure with the HPM table and will result in incorrect dynamics in the simulation. We must treat the smoothing and clumping effects carefully when using the HPM variables calculated for each gas particle to model their thermal properties in the simulation. In Figure 4 we plot the temperature profiles of simulated halos within the mass bin centered at M 200c = 3 × 10 14 M e at redshift z = 0.5, where the temperature of the gas particles is interpolated from the HPM table built with the original ICM model. Ignoring the smoothing and clumping effects on HPM variables results in underestimating the temperature of gas particles in the center of the halo and in overestimating it in the outer region. The deviation from the ICM model prediction on the gas temperature could be up to 20%. Mapping relation from the HPM variables matter density, ρ m , and the scalar force, f scalar , discussed in Section 3.3, to the gas temperature (left) and pressure (right). Both panels show that the HPM variables have a significant correlation with the target mapped quantities, though the mapping relation for the gas temperature has a more complicated pattern in the top-right corner due to the non-monotonicity of the temperature profile in the core region of halos in the ICM model.
The reason for this mismatch is that the smoothing effect in the inner region of simulated halos suppresses magnitude of the HPM variables to the values below the ICM model prediction. According to the HPM table we show in Figure 2, the HPM variables have a significantly positive correlation with the temperature. Thus, underestimating the HPM variables means we are also very likely to underestimate the inferred temperature simultaneously. Then gas particles that reside in the halo center are more likely to be assigned a temperature value that is biased low. In the outer region of simulated halos, the clumping effect leads to HPM variables for a substantial number of particles being higher than the ICM prediction. Due to the same reason that temperature and HPM variables are positively correlated in the HPM table, overestimation of the HPM variables leads to the temperature of the gas being too high in the outer region.
To mitigate this inconsistency in the temperature inference caused by the limited resolution and breakdown of spherical symmetry, we need to take the smoothing and clumping effects into account while building the halo table. We first run a darkmatter-only simulation and fit a calibration function C(r, ζ) to the ratio of simulated profiles of HPM variables to their ICM model predictions, which enables us to capture the differences between the simulation results and the original ICM model. We have discussed in Section 3.3 that the HPM variables are only mildly affected by the gas distribution, so the profiles of HPM variables derived from the dark-matter-only simulation can properly emulate their values in an HPM simulation. We verified that the offsets between simulated profiles of HPM variables and their ICM model predictions measured in the dark-matter-only simulation and HPM simulation were indeed very similar.
The calibration function measured from a dark-matter-only simulation can effectively correct the gas temperature or pressure inference in HYPER. We adjust the calculation with the fitted calibration function when deriving the HPM variables in the halo table aŝ   where ζ refers to any information needed to specify numerical resolution. The calibration function is fitted to the simulation results for two mass bins M 200c = 10 14 M e , 3 × 10 14 M e at redshift z = 0.0, 0.5, and we do not find evident mass or redshift dependence . The HPM table built with the adjusted  halo table more accurately relates the HPM variables of simulated particles to their thermal properties and better approximates the hydrodynamics of gas particles through the ICM model. The profiles of HPM variables adjusted for the calibration functions are also shown in Figure 3. More details about the calibration functions are presented in Appendix B. After adjusting the HPM table for the effects mentioned above, we show in Figure 4 that the temperature inferred from the adjusted HPM table is in good agreement with the ICM prediction. We observe a downturn in the radial temperature profile beyond R 200c , which is due to the existence of a substantial amount of low-temperature IGM gas at the outskirts of the massive halos.

Artificial Viscosity and Pressure Filtering
Shocks are a generic feature of gas flows. When solving the fluid equations in a particle-based simulation like SPH, the entropy generation on shocks is captured by an artificial viscosity term. HYPER is a Lagrangian particle method for the dynamical evolution of gas similar to SPH, and hence needs an artificial velocity term in the gas momentum equation: 3 8 gas gas visc Though we do not need to resolve the conversion between kinetic energy and thermal energy in HYPER as the thermal properties of gas are directly inferred from the HPM table built on the ICM model, we still need to include the artificial viscosity term to prevent particle interpenetration in shocks. The artificial viscosity in HYPER removes the part of the kinetic energy that should be converted into heat. Otherwise, we will not be able to accurately resolve the motion and distribution of the gas in our simulation (Appendix C). Many different forms have been suggested for the artificial viscosity, with the Von Neuman-Richtmyer artificial viscosity being the simplest one: where h is the particle smoothing length, proportional to the local mean inter-particle separation, r µh ; gas 1 3 and c is the speed of sound in gas. In this work, α ∼ 0.1, β ∼ 0.05 are found to be the best-fit values when calibrating the gas radial profiles and scaling relations of different integrated quantities of the simulated halos. The viscosity term we adopt has a similar form as the gas pressure force term ∇P/ρ gas , with Q acting as an excess pressure assigned to gas particles. To integrate the artificial viscosity in HYPER, we only need to solve the dynamical equation for gas particles in Equation (23) with the modified pressure ∇(P + Q)/ρ gas . In our HPM simulation, in addition to adopting the artificial viscosity for gas particles, we also apply a force filter on the pressure force ∇P/ρ gas . This approach aims to sustain a proper hydrostatic equilibrium in the central core region of simulated halos, since the interpolated temperature of gas particles in the center of the halos is not affected as much by the smoothing effect of finite resolution as the gas density, and that in turn leads to less smoothing of the gas pressure compared to the gravitational force. Filtering on the small-scale structure of the pressure field will impose the same degree of smoothing effect on the hydro-force as the gravitational force in the central region of the simulated halos. Filtering also suppresses the numerical noise introduced by the interpolation of gas temperature from a model-based HPM table. We filter the gas pressure field in the Fourier space bỹ f and the filtered pressure force is given by replacing the gas pressure P in ∇P/ρ gas with the filtered pressure P f . For the specific functional form for the filtering, we adopt the Weibull function for the high-frequency filter. Parameters A f L and A f S are tuned to match the simulated results with the ICM model predictions: . A combination of the artificial viscosity and filtering on gas pressure helps to prevent too much gas from being ejected from the halos in our simulation and to ensure we get a reasonable gas fraction for the simulated halos. Tuning of the small-scale filter and artificial viscosity is presented in Appendix C.

Results
In this section we compare the results from an HYPER simulation to the predictions of the ICM model used to implement the HPM mapping relation for gas thermal properties. Good consistency in the properties of the ICM such as the radial profiles, scaling relation of integrated halo properties, and measurements of the tSZ effect implies one can systematically control the ICM physics in the HYPER simulation by varying the ICM model while constructing the HPM mapping relation.

Halo Radial Profiles
We first compare the profiles of matter density and gas thermal properties of simulated halos to the profiles derived with the ICM model. In Figure 5 we show the radial profiles of matter density and their scatter for the simulated halos whose masses fall into three mass bins centered at 10 14 M e , 3 × 10 14 M e , and 10 15 M e at redshifts z = 0 and z = 0.5. We find that the volume-weighted density profiles agree with the NFW model well except in the inner region, where they suffer from the smoothing effect due to the limited resolution of the simulation. The clumping effect in the mass-weighted density profiles mentioned in Section 3.5 leads to an overestimate of the matter density at the outskirt region of the halos.
In Figure 6 we plot the simulated profiles of gas density, temperature, and pressure and their scatter for the halos whose masses fall into the mass bin centered at 3 × 10 14 M e at redshifts z = 0 and z = 0.5. We also plot the prediction for the gas profiles in the ICM model used for implementing the hydro part of the simulation for comparison, and we also show the ratio of simulated profiles to the ICM model predictions in the bottom panels. We find that the radial profiles of the gas density, temperature, and pressure for simulated halos are in about 5% agreement with the ICM model from 0.1R 200c to R 200c and remain within 20% even in the inner core and in halo outskirts up to about 1.5R 200c .
In the inner core, we observe the smoothing effect due to the limited resolution of the simulation in the gas density and pressure profiles. An analogous smoothing effect in the simulated gas temperature profiles in the halo centers is much less significant because we already account for the numerical resolution effect when constructing the HPM table. Notice that we even account for the numerical effect on the HPM variables used to infer the thermal properties of gas and obtain the hydroforce less affected by the numerical effect in HYPER simulation. The smoothing effect will still affect the gravity calculated with the PM scheme and further influence the gas and dark matter density distribution. This explains the smoothing effect we observed in the gas and total matter density profile and gas pressure profile.

Integrated Halo Quantities
X-ray observables, such as luminosity, temperature, mass of the ICM, and SZ flux of galaxy clusters have been proposed and used as proxies for the total cluster mass (e.g., Voit 2005).
Calibrating relations between cluster mass and these observables is important for exploiting the full statistical power of the cluster surveys. In this section we show further comparison between the simulation results and the ICM model predictions by investigating the integrated quantities of identified halos and exploring the scaling relation between the SZ effect signal, X-ray observables, and halo masses in the simulation. For a more detailed comparison between the simulation results and the ICM model expectations, we split the simulated halos into eight groups by mass bins, which uniformly span over the mass range 10 13.7 M e − 10 15.3 M e in logarithmic scale and with 0.2 dex for the width of log bins. We choose this mass range for the analysis because the halos below the low-mass end (5 × 10 13 M e ) significantly suffer a resolution effect, and we only see few massive halos above the high-mass end (2 × 10 15 M e ). We calculate the mean and the uncertainty of the gas fraction of simulated halos falling into each mass bin. For the SZ effect signal, X-ray observables, we calculate the mean and the uncertainty of these quantities of the simulated halos belonging to each mass bin in a logarithmic scale. We ignore the results for mass bins containing less than 10 halos, as they lack statistical reliability.
In Figure 7, we plot the gas fractions of simulated halos enclosed within spheres with the averaged density of 200 and 500 times the critical density, respectively, at different redshifts z = 0, 0.5, and 1.0. In addition, we plot the mean and uncertainties of the gas fraction of simulated halos for different mass bins. For comparison, we plot the gas fraction derived from the ICM model, which equals the ratio of the integral of Figure 5. Comparison between the simulation results and the ICM model (NFW profile) for the radial profiles of matter density for halos in the mass bin centered at 10 14 M e (first column), 3 × 10 14 M e (second column), and 10 15 M e (third column) at redshifts z = 0 (first row) and z = 0.5 (second row). In each plot, the top panel shows that the simulated volume-weighted matter density profile (blue line) agrees well with the NFW profile (black line). In contrast, the mass-weighted matter density profile (red line) overestimates the matter density in the outskirt region compared to the standard ICM model due to the clumping effect. Also shown is the scatter in the simulated matter density profiles (thin blue/red bands); the bottom panel shows the ratio between the simulation results and the ICM model and a ±20% region (dotted black lines). gas density given by Equation (11) over certain over-density regions to the total halo mass. The universal baryonic fraction where Ω b and Ω m are the cosmological parameters set for the simulation, is also shown. The plots also show that, after implementation of gas pressure in solving the motion of baryonic components in the simulation, the gas fraction of the Figure 6. Comparison between the simulation results and the ICM model prediction for the volume-weighted radial profiles of gas density (first column) and pressure (third column) and mass-weighted temperature (second column) for halos in the mass bin centered at 3 × 10 14 M e at redshifts z = 0 (first row) and z = 0.5 (second row). For each plot, the top panel shows the simulated halo profile and its scatter (blue line and band), which is found to be in good agreement with the profile of the ICM model (black line). The bottom panel shows the ratio between the simulation results and the ICM model and a ±20% region (dotted black lines). simulated halos is substantially lower than the universal baryonic fraction f b . If assuming the evolution of gas is only affected by gravity, the gas particles will strictly follow the dark matter particles in simulation in a pair-wise pattern. And the gas fraction within any region in simulation equals the fraction that the baryon mass takes of the mass of all of the matter in the universe, which is also the definition of the universal baryonic fraction f b . The difference between the simulation results and f b indicates a considerable portion of gas gets propelled out of the gravitational potential well of collapsed halos due to the existence of gas pressure. Furthermore, we find the simulation results of gas fraction match the prediction derived from analytical models of gas pressure profiles by assuming hydrostatic equilibrium reasonably well. For each mass bin, the mean gas fraction of simulated halos agrees with the model derivation within 10%. Scatter in the results of simulated halos is not negligible, as the dynamical state of different halos can vary significantly in a real simulation. We observe slightly more scatter inside the radius encompassing the over-density of 500 times the critical density, since the inner region is more sensitive to the dynamic state of the halos. The failure of the simulation to resolve the inner cores of halos may also contribute to the overall scatter.
In Figure 8, we plot the Compton Y-parameter integrated over a spherical volume, on average for different mass bins is found to be slightly larger than the scatter of Y 500c ,  s 0.06 Y log 10 500c , which is due to greater intrinsic scatter in the pressure profile in the outer region of the simulated halos, as shown in Figure 6. Similar over-density dependence of the scatter in the Compton Y-parameter is also observed in the full hydro-simulations (e.g., Battaglia et al. 2012b;Gupta et al. 2017). Scatter of the Y 500c in the simulation is smaller than scatter of the X-ray observable L X that  s 0.13 L log 10 X on average for different mass bins shown in Figure 9; apparently the Compton Y-signal of the simulated halos is less affected by the poorly resolved inner-core region. Our results agree with the conclusion in other studies of the SZ effect (e.g., Komatsu & Seljak 2002) that the Compton Y-signal is less sensitive to gas physics in the cores of clusters.
In Figure 9 we plot the bolometric X-ray luminosity L X , T is the cooling function assuming that Bremsstrahlung emission dominates, and n gas is the gas Another common X-ray observable is the Y X parameter (Kravtsov et al. 2006), proportional to the gas thermal energy as defined by the product of the gas mass and the spectroscopic X-ray temperature, where M gas is the gas mass within the spherical over-density region of radius R 500c , and T X is the spectroscopic X-ray temperature. Finally, we also show in Figure 9 the characteristic temperature T 500c , The scaling relation of the integrated X-ray quantities and halo mass in the simulation results are in good agreement with the ICM gas model predictions from redshift z = 0 to z = 1. We Figure 9. Integrated X-ray quantities: the bolometric X-ray luminosity (top left), the spherical Compton-like Y x -parameter (top right), the mass-weighted average temperature (bottom left), and the emission-weighted temperature (bottom right) of the simulated halos (dots) compared with the ICM model prediction (dashed lines) at different redshifts z = 0 (blue), z = 0.5 (red), and z = 1 (green). The mean and uncertainties in logarithmic scale (square dots with error bars) of the X-ray observables of simulated halos within different mass bins are plotted for more strict comparison with the ICM model derivation. The bottom panel shows the difference between the mean and uncertainty of the X-ray observables of the simulated halos within each mass bin in logarithmic scale and the corresponding ICM for most of the mass bins supports a good agreement between the simulation results and ICM model derivations within 10%, where Y could denote different X-ray observables L X , Y X , T 500c , and T ew . Notice that all of these X-ray observables are exclusively defined by the density and the temperature of the gas in the halos. Since we have shown in Figure 6 that the radial profiles of both the gas density and temperature are in good agreement with the analytical model, this explains the consistency between the integrated X-ray quantities of the simulated halos and the ICM model prediction. Simulation results for the Y X have been shown to have less scatter,  s 0.07 Y log 10 X on average for different mass bins, than the L X , which may be due to the anticorrelation between the deviation of M gas and T X from the scaling relation prediction, which has been found in X-ray observations (Kravtsov et al. 2006).
We also look into the scatter about the scaling relation in the HYPER simulation, as it could reflect the information on important dynamical effects. Gupta et al. (2017) studies the scatter about the Y-M relation in the Magneticum simulation and finds that the scatter s Y ln about the mass-observable relations at over-density 500c is consistent with the lognormal distribution. We also calculate the distribution of scatter about the Y 500c -M 500c relation in HYPER and fit it to a lognormal distribution. We find that the rms scatter s =  0.094 0.002 Y ln is in good agreement with the value s =  0.088 0.006 Y ln from the Magneticum simulation. We show the probability distributions of scatter about the Y 500c -M 500c for HYPER and its best lognormal fit, as well as the best-fit lognormal from the Magneticum simulation in Figure 10. We also measure the scatter about the X-ray observable-mass relation by calculating the rms dispersion 10 10 model 2 10 following Barnes et al. (2017), where Y i is the X-ray observables as i runs over all simulated halos, and Y model is the corresponding ICM model prediction. We present the comparison between our measurement for the scatter about L X -M 500c , Y X -M 500c , and T 500c -M 500c and the MACSIS simulation results in Table 1.
We find the scatter about the scaling relation in HYPER is in general consistent with the state-of-the-art full hydro-simulations, suggesting that HYPER captures the important dynamical effects modeled in full simulations. However, to further verify this point, we need to adopt the ICM model drawn from the full hydro-simulations to construct the HPM table for inferring the gas thermal properties. Studies show that the ICM physics could also affect the scatter about the scaling relation (e.g., Battaglia et al. 2012b). The excellent agreement on the scatter about the Y 500c -M 500c relation between HYPER and Magneticum may result from the ICM gas pressure model we adopt for HYPER being similar to that found in the Magneticum simulation. This implies that HYPER adopting the current ICM model may properly emulate the ICM physics in the Magneticum simulation, and the good match on the scatter is a result of a more fair comparison. Hence, a more detailed comparison for scatter about the observable-mass relation of HYPER and full hydrodynamic simulations is required in the future study.

Thermal SZ Angular Power Spectrum
The tSZ power spectrum C l is a powerful probe of cosmology and can provide promising constraints on cosmological parameters, in particular σ 8 , since s µ -C l 8 7 9 (e.g., Komatsu & Seljak 2002;Shaw et al. 2010;Trac et al. 2011). Because the cluster signal dominates tSZ anisotropies, we can model the analytical prediction of the tSZ power spectrum using the standard halo model (e.g., Cole & Kaiser 1988;Komatsu & Kitayama 1999). The tSZ angular power spectrum at a multipole moment l for the one-halo term is given by is the spectral shape of the tSZ signal. Integration over redshift and mass are carried out from z = 0 to z = 5 and from M = 10 10 M e to M = 10 16 M e , respectively. For the differential halo mass function dn(M, z)/dM, we adopt the fitting function from Tinker et al. (2008) based on N-body simulations.
Following Komatsu & Seljak (2002), the 2D Fourier transform of the projected Compton Y-parameter,˜( where x = r/r s is a scaled dimensionless radius, r s is a characteristic radius for an NFW profile defined as R 500c /c 500c , and we use the average halo concentration c 500c calibrated as a function of the cluster mass and redshift from Diemer & Kravtsov (2015). The corresponding angular  and we use the FFT to calculate the power spectrum of the gas pressure field drawn from the simulation snapshots at redshift z i to directly estimate the thermal SZ angular power spectrum from our simulation output. In Figure 11 we show the tSZ power spectrum measured from the HYPER simulation, which is calculated with the 3D power spectrum of the gas thermal pressure. We compare our result to the analytical prediction evaluated with the ICM pressure model and the halo model assuming the same fiducial cosmological parameters as adopted by simulation. We also plot the measurements of the tSZ power spectrum by Planck, which are in good agreement with the analytical calculation of the tSZ power spectrum derived from DPP in He et al. (2021). The frequency-dependent terms are all scaled to f 2 (ν) = 1 for direct comparison.
The tSZ power spectrum measured with the output of the HYPER simulation is in good agreement with the analytical prediction calculated with the ICM model, within ∼10%. It shows that the HYPER simulation results are consistent with the assumed ICM model used to implement the hydrodynamics of HYPER not only for properties of simulated halos, but also for the most widely used statistical measures of the tSZ effect. According to Equation (49), Compton signals of galaxy clusters dominate the contribution to the tSZ power spectrum, and we show in Section 4.2 that the scaling relation between the Compton Y-parameter and the mass of simulated halos matches the ICM model derivation very well, which also indicates a good agreement on the tSZ power spectrum calculation.

Conclusion and Discussion
In this paper we introduce HYPER, a new implementation of a fast approximate hydro-simulation based on an N-body solver. HYPER applies a power-law density-temperature relation for the gas in the IGM regime of low-density and constructs a mapping relation between two designed HPM variables and the gas temperature and pressure in the highdensity ICM regime (which is based on the adopted ICM gas pressure model) to simulate the evolution of baryonic matter in an efficient way.
We investigate the properties of gas inside the simulated halos by measuring the radial profiles of density, temperature, and pressure of the gas for the identified halos in the simulation. We also present the integrated quantities of observables in the X-ray and SZ survey and calculate the tSZ angular power spectrum from the simulation outputs. We emphasize that one of the crucial strength of HYPER is that one can systematically control the gas physics of simulated halos with the adopted ICM model. We show the following: 1. The radial profiles of density, temperature, and pressure of the gas inside the identified halos in the simulation are in good agreement with the ICM model predictions within 5% for 0.1R 200c -R 200c . The deviation is limited to 20% in the inner core r < 0.1R 200c and the outskirt region R 200c -1.5R 200c . Mild inconsistency found in the inner region might be due to the resolution limit of the HPM solver. 2. The integrated X-ray and SZ observables of the simulated halos are in good agreement with the scaling relations derived from the gas radial profiles of the ICM model at redshifts from z = 0 to z = 1. The scatter in the relation comes from two primary sources: the variations in the dynamical states of different simulated halos and the finite numerical resolution in the inner cores of clusters. The latter may also contribute to the bias of the gas fraction of simulated halos as compared to the model prediction. The scatter in the simulation results of Figure 11. (Top) tSZ angular power spectrum evaluated with the 3D power spectrum of the gas thermal pressure and Equation (53)  different observables is comparable to that in the full hydrodynamic simulation Magneticum and MACSIS. 3. The tSZ angular power spectrum measured for the HYPER simulation, which is calculated using the 3D power spectrum of the gas thermal pressure drawn from the simulation at different redshift snapshots, is in good agreement with the analytical predictions evaluated with the halo model and the ICM model used to implement the HPM algorithm. Good consistency in the simulation output and ICM model derivation for properties of the ICM regime includes the cluster radial profiles, SZ, and X-ray observable-mass relation. Statistical quantities of the tSZ effects indicate that the HYPER simulation allows us to systematically control the ICM physics by varying the ICM model implemented in the HPM mapping relation construction.
We envision three main-use cases for HYPER.
1. HYPER runs orders of magnitude faster than ordinary hydrodynamic simulations. It can be useful for generating a large number of mock catalogs and creating maps of various physical quantities (X-Ray, temperature, SZ effect, etc.) for galaxy clusters. These outputs will help in the development of data reduction and analysis pipelines, for understanding systematics and selection effects, and for interpreting cosmological and astrophysical constraints. One can also envision training a multiband deep learning model with mock observations generated by HYPER for the mass estimation and mass distribution measurement of galaxy clusters. 2. In HYPER, we implement the dynamics of baryons via the HPM mapping relation built on the ICM model for the gas profiles. Modifications of the matter power spectrum due to baryonic physics are one of the major theoretical uncertainties in cosmological weak lensing measurements. The ability of HYPER simulations to efficiently model the joint effects and varied cosmological parameters is a powerful tool for studying mitigation schemes for baryonic effects in weak lensing cosmic shear measurements (e.g., Huang et al. 2019). 3. Baryonic physics such as star formation, energetic feedback, and nonthermal pressure support affect the tSZ angular power spectrum in nontrivial ways (e.g., Trac et al. 2011). The difference in gas physics is imprinted in the cluster gas pressure profile. With HYPER, one can systematically vary the gas pressure model of galaxy clusters. For example, one can implement different ICM models drawn from current state-of-the-art high-resolution hydrodynamic simulations in a large-scale fast hydro-simulation of size up to ∼1-2 Gpc, which would be unrealistically expensive for ordinary hydrodynamic simulations. With this approach, one can systematically study how different gas physics influences, for example, the tSZ angular power spectrum. Moreover, from an inverse perspective, the efficiency of HYPER also enables us to generate a large number of mock observations; when combined with proper statistic techniques (e.g., Gaussian process), this allows us to quickly explore parameter space for the ICM model and use observational data of the SZ effect to put constraints on the ICM model. Furthermore, it can be used to examine the gas physics implemented in cosmological simulations.
In future extensions of this work, we will focus on improving the finite spatial resolution of HYPER. One avenue could be adopting a hybrid scheme combining the multigrid method with FFT (e.g., Kravtsov et al. 1997), which could eliminate the resolution effect in the highdensity ICM regime and sustain the high computational efficiency throughout the rest of the simulation volume. This may help solve the problem of overestimating the gas velocity dispersion in the inner region of halos. We consider using the two-level PM scheme (e.g., Trac & Pen 2006) to achieve higher spatial resolution while significantly reducing memory for future development. We also plan to study the performance of HYPER compared to other large-scale stateof-the-art hydrodynamic simulations like Magneticum and Illustris by replacing the ICM pressure model adopted in this work with a simulated ICM pressure profile drawn from the corresponding hydrodynamic simulation. That would allow for a direct comparison between HYPER and full hydrosimulations, and would enable us to study further the reliability of our new HPM algorithm.
We thank Renyue Cen, Markus Michael Rau, and Qirong Zhu for helpful discussions. H.T. acknowledges support from the NSF AI Institute: Physics of the Future, NSF PHY-2020295. Computer simulations and analysis were supported by Pittsburgh Supercomputing Center under grant AST190014P, PHY200019P and PHY200058P.

Appendix A Computational Performance of HYPER
In Figure 12, we plot the CPU time of HYPER simulations of different sizes, with numbers of simulated particles N p = 2 × 128 3 , 2 × 256 3 , 2 × 512 3 , and 2 × 1024 3 . We show that HYPER simulations take approximately two to three times as long to run as standard PM simulations of the same size. Since the gravity and HPM variables for inference of hydro-force in HYPER are calculated with the fast Fourier transform algorithm, the CPU time of the HYPER simulation is in general scaled by the change of the number of particles with the expected ( ) N N log p p scaling, the same as in the PM simulation. Most of the other particle calculations take a comparable time and scale with the N p . We also compare the HYPER simulation with the current state-of-the-art full hydro-simulations: the CPU time of TNG100 and TNG300 simulations of IllustrisTNG is reported in Pillepich et al. (2018b), the Box2 and Box3 simulations of Magneticum are reported in proposal 6 , and the L100N1504 simulation of EAGLE is reported in Schaye et al. (2015). All of these contain a total number of simulated particles (N p ∼ 10 8 -10 10 ) and are of a similar order of magnitude to the largest run of HYPER simulation (N p ∼ 10 9 ). We find HYPER simulations are orders of magnitude faster than the expensive hydrodynamic simulations when interpolated or extrapolated to the same number of particles as these simulations following ( ) N N log p p scaling. HPM Variables Adjustment We adopt the form of the Weibull function for the calibration function C(r, ζ) to characterize the difference between the simulated radial profiles of HPM variables and their theoretical values in the ICM model: where Y denotes the HPM variables ρ m and f scalar , and ζ = Δl = L/N is the length of the grid cell in our simulation. We fit the calibration function to the ratio of the simulated mass-weighted profiles of ρ m and f scalar to the radial profiles derived in the standard ICM model for two mass bins M 200c = 10 14 M e , 3 × 10 14 M e at redshifts z = 0.0, 0.5. The importance of the data points in fitting is weighted by the uncertainties of the radial profiles. The best-fit results for the parameters A A A , ,  1.30,0.48,0.60] for HPM variables ρ m and f scalar , respectively. These parameters of the calibration functions may be spatialresolution dependent.
In Figure 13 we plot the ratio of the simulated radial profiles of ρ m and f scalar to their respective values in the adopted ICM model and its uncertainty. We also show the best-fit results for our calibration function. The fitted calibration functions agree well with the simulated mean profiles in the inner regions of halos, but show more deviation as the uncertainties of the simulation results grow in the outskirt region.

Appendix C Pressure Filtering and Artificial Viscosity Tuning
As discussed in Section 3.6, we apply a pressure filter on the gas pressure field in Fourier space to smooth the gas pressure in the inner region of the halos and suppress the numerical noise. The pressure filter is in the form of the Weibull function expressed by Equation (42). We tune the parameters A L and A S , which determine the scale range and the degree of smoothing at small scales for the pressure filter, to make the simulated gas density profiles match the ICM profiles.
In Figure 14 we compare the simulated gas density profiles with the analytical derivation in the ICM model after filtering the gas pressure field in Fourier space. We experiment with using the pressure filters with different combinations of A f S and A f L . We find that the simulation results are less sensitive to the filter parameter A f S , and we choose to adopt = A 10 f S for our simulation. We tune the A f L starting from = A 1.0 f L then gradually decrease to adjust the smoothing effect on the gas pressure imposed by the filter until we can balance the hydroforce and the gravity and get a proper density profile when = A 0.3 f L . However, with only the pressure filter, we can obtain a proper density profile but overestimate the velocity dispersion for gas particles in the simulated halos. By applying the Figure 12. The CPU time of HYPER simulations of different sizes, with the number of simulated particles N p = 2 × 128 3 , 2 × 256 3 , 2 × 512 3 , and 2 × 1024 3 (hexagonal dots connected by the solid black line), is found approximately follow the ( ) N N log p p scaling with the total simulated particle number. HYPER simulations take around two to three times as long to run as the standard PM simulations of the same size of total particle number (star dots connected by the dashed black line). The CPU time of full hydro-simulations is also plotted for comparison: TNG100 (red dot), TNG300 (red circle), Magneticum Box3 (blue square), Magneticum Box2 (blue box), and EAGLE (green triangle). All of these are of similar size, with the most extensive realization of HYPER simulation concerning the particle number. HYPER simulations are shown to be orders of magnitude faster than the expensive hydrodynamic simulations. artificial viscosity, which removes the extra kinetic energy of gas particles that should have been converted into heat, we can reduce the discrepancy, although not eliminate it completely. As shown in the right panel of Figure 14, smoothing the gas thermal pressure field avoids an excessive hydro-force in the inner region of the simulated halos and prevents pushing too much gas out of the halo. We also find that the artificial viscosity can effectively reduce the gas velocity dispersion, and it decreases the nonthermal pressure. Since both pressure filter and artificial viscosity can mitigate the gas pressure to balance the hydro-force and gravity, we can tune the pressure filter and the artificial viscosity jointly to keep the simulated gas density profile matching the model derivation and reduce the gas velocity dispersion at the same time. In this case, we find the pressure filter with slightly larger A f L that = A 0.5 f L and α ∼ 0.1, β ∼ 0.05 for the parameters of artificial viscosity in Equation (40) work best for our simulation.
In Figure 15 we plot the ICM gas profiles with (a) both pressure filter and artificial viscosity and (b) the pressure filter only. We show that by tuning the pressure filtering and the . Also shown are the best fits for the calibration functions C(r, ζ) (solid black line). Fitted calibration functions are found to be in good agreement with the simulation results, while deviation appears at large radii where the uncertainties for the simulated radial profiles become much more significant, and data points are less important in the fitting. The bottom panels show the ratio between the simulation results and the best fits for the calibration functions and the ±20% band (dotted black lines). artificial viscosity jointly, we can make the gas density profile match the ICM model and obtain a lower gas velocity dispersion in the halo outskirt region. Even with the artificial viscosity included in our simulation, it is still difficult to match the gas velocity dispersion in the inner region of the simulated halos to the ICM model. This may be because the finite resolution of the HPM algorithm limits the ability to accurately resolve the local artificial viscosity exerted on the gas particles in our simulation, especially in the highdensity halo central region where hundreds, even thousands, of particles can colocate in one grid cell. Further studies are required to explore the reason for the discrepancy and improve the simulation fidelity. One possible solution to increase the resolution in high-density regions is to use the adaptive mesh refinement.