# The Impacts of Modeling Choices on the Inference of Circumgalactic Medium Properties from Sunyaev-Zeldovich Observations

Emily Moser<fnref n="1" />, Stefania Amodeo<fnref n="1" />,<fnref n="2" />, Nicholas Battaglia<fnref n="1" />, Marcelo A. Alvarez<fnref n="3" />,<fnref n="4" />, Simone Ferraro<fnref n="3" />,<fnref n="4" />, and Emmanuel Schaan<fnref n="3" />,<fnref n="4" />

## Abstract

As the signal-to-noise of Sunyaev-Zeldovich (SZ) cross-correlation measurements of galaxies improves our ability to infer properties about the circumgalactic medium (CGM), we will transition from being limited by statistical uncertainties to systematic uncertainties. Using thermodynamic profiles of the CGM created from the IllustrisTNG (The Next Generation) simulations we investigate the importance of specific choices in modeling the galaxy sample. These choices include different sample selections in the simulation (stellar versus halo mass, color selections) and different fitting models (matching by the shape of the mass distribution, inclusion of a two-halo term). We forward model a mock galaxy sample into projected SZ observable profiles and fit these profiles to a generalized Navarro-Frenk-White profile using forecasted errors of the upcoming Simons Observatory experiment. We test the number of free parameters in the fits and show that this is another modeling choice that yields different results. Finally, we show how different fitting models can reproduce parameters of a fiducial profile, and show that the addition of a two-halo term and matching by the mass distribution of the sample are extremely important modeling choices to consider.

Unified Astronomy Thesaurus concepts: Circumgalactic medium (1879); Sunyaev-Zeldovich effect (1654)

## 1. Introduction

The circumgalactic medium (CGM) consists of a large reservoir of gas surrounding the outermost regions of galaxies. It is an important component of the galactic structure, and is believed to govern the cyclical flow of material in and out of galaxies. This cycle includes material falling onto the disks of galaxies from the CGM and intergalactic medium (IGM) and being recycled back out through various feedback mechanisms.

The CGM is known to be complex and multiphase in temperature, ionization state, and kinematics (Tumlinson et al. 2011; Werk et al. 2013; Nielsen et al. 2015; Tumlinson et al. 2017; Oppenheimer et al. 2018) and has been studied through observations (e.g., Lanzetta et al. 1995; Chen et al. 1998; Steidel et al. 2010; Tumlinson et al. 2011; Rudie et al. 2012; Tumlinson et al. 2013; Werk et al. 2014) and hydrodynamic simulations (e.g., Oppenheimer & Davé 2008; Ford et al. 2013; Hummels et al. 2017; Suresh et al. 2017; Oppenheimer et al. 2018; Hummels et al. 2019; Peebles et al. 2019). Several questions regarding the physical processes that govern the different states of the CGM remain unanswered, commonly including what are the properties of quenched galaxies and how do the galaxies become quenched (Baldry et al. 2004; Noeske et al. 2007; Peng et al. 2010; Geha et al. 2012), and the missing baryon problem in which baryons predicted by the current cosmological model are not observed (Persic & Salucci 1992; Fukugita et al. 1998; Fukugita & Peebles 2004; Cen & Ostriker 2006; Bregman 2007; Davé 2009). Star formation in galaxies is known to be inefficient (Fukugita & Peebles 2004; Federrath 2015), which implies that there is some process occurring within the CGM and its interaction with the host galaxy that is preventing the gas from cooling and forming stars efficiently. This process could be related to the various feedback mechanisms injecting energy into the system such as supernovae and active galactic nuclei (AGN;

Di Matteo et al. 2005; Springel et al. 2005; Scannapieco et al. 2008; Somerville et al. 2008; Marasco et al. 2015; Somerville & Davé 2015; Agertz & Kravtsov 2016) that push central material to the far parts of the galaxy while also affecting its thermodynamic properties. These problems and many more discussed in the literature are important to understand for galaxy evolution theory, and can be studied theoretically through cosmological simulations and observationally through galaxy surveys at multiple redshifts and wavelengths.

An emergent method to observe the CGM is observing the cosmic microwave background (CMB). While initially measured to be homogeneous and isotropic, fluctuations (called anisotropies) in temperature have been detected, and these anisotropies can give insight on the seeds of structure formation and growth. Secondary anisotropies due to the Sunyaev-Zeldovich (SZ) effect (Sunyaev & Zeldovich 1970) have also recently been measured, which can be further broken down into thermal (tSZ) and kinetic (kSZ) effects (Sunyaev & Zeldovich 1972, 1980).

The tSZ effect describes the increase in energy of CMB photons due to scattering off ionized electrons in galaxies and galaxy clusters, and produces distortions in the nearly perfect blackbody spectrum of the CMB. It is a function of frequency, $\nu$, and the Compton-y parameter, shown in Equation (1). The amplitude of this distortion is proportional to the line-of-sight (LOS) integral of the electron pressure, so it can be used as a probe of the gas pressure within galactic halos.

$$ \begin{aligned} \frac{\Delta T(\nu)}{T_{\text{CMB}}} &= f(\nu)y(\theta), \\ y(\theta) &= \frac{\sigma_T}{m_e c^2} \int_{\text{LOS}} P_e(\sqrt{l^2 + d_A^2(z)|\theta|^2}) dl, \end{aligned} \quad (1) $$

where $\Delta T(\nu)$ is the shift in temperature measured as the tSZ signal, $T_{\text{CMB}}$ is the temperature of the CMB, spectral function

<footnote n="1">¹Department of Astronomy, Cornell University, Ithaca, NY 14853, USA</footnote>

<footnote n="2">²Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote n="3">³Lawrence Berkeley National Laboratory, One Cyclotron Road, Berkeley, CA 94720, USA</footnote>

<footnote n="4">⁴Berkeley Center for Cosmological Physics, UC Berkeley, CA 94720, USA</footnote>

<footnote>Received 2021 March 26; revised 2021 May 28; accepted 2021 June 17; published 2021 September 16</footnote>
$f(\nu) = x\coth(x/2) - 4$, with $x = \frac{h\nu}{k_B T_{CMB}}$, $h$ is the Planck constant, $k_B$ is the Boltzmann constant, $y(\theta)$ is the Compton-y parameter measured within $\theta$, $\sigma_T$ is the Thomson scattering cross section, $m_e$ is the electron mass, $c$ is the speed of light, $d_A(z)$ is the angular diameter distance at redshift $z$, and $P_e$ is the electron pressure.

The kSZ effect is the Doppler shift of CMB photons scattering off free electrons in galaxies and clusters with peculiar velocities, causing Doppler shifts in the CMB temperature that are directly related to the peculiar momentum. These shifts in temperature are proportional to the LOS integral of peculiar velocity multiplied by electron number density; therefore, we can use the kSZ to measure the density of the CGM in extragalactic halos (Battaglia et al. 2017).

$$ \frac{\Delta T}{T_{\text{CMB}}} = \frac{\sigma_T}{c} \int_{\text{LOS}} e^{-\tau(\theta)} n_e v_p dl, \\ \tau(\theta) = \sigma_T \int_{\text{LOS}} n_e (\sqrt{l^2 + d_A^2(z) |\theta|^2}) dl, \quad (2) $$

where $n_e$ is the electron number density, $v_p$ is the peculiar velocity, and $\tau(\theta)$ is the optical depth.

Combining the tSZ and kSZ effects provides complete thermodynamic information of the CGM (Battaglia et al. 2017) and can provide constraints on the physical processes that govern star formation and feedback, and thus galaxy evolution. The combination of SZ measurements is a relatively new method to study the CGM. While the tSZ effect has been observable in galaxies and galaxy clusters for several years, detections of the kSZ effect have only recently been possible through cross-correlation of CMB observations with galaxy catalogs (e.g., first detected in Hand et al. 2012). Among the kSZ estimators that exist in the literature, a recent powerful method uses a velocity-weighted stack, in which the signal measured for each halo on the CMB map is weighted by an estimate of its LOS velocity, reconstructed from the density field (e.g., Planck Collaboration et al. 2016a; Schaan et al. 2016, 2021). Using this method, recent results from the Atacama Cosmology Telescope (ACT) DR5 and Planck (Amodeo et al. 2021; Schaan et al. 2021) have achieved high signal-to-noise measurements of the kSZ and tSZ effects, and in turn, of the electron density, temperature, and pressure distribution around the CMASS (constant stellar mass) galaxies from the Baryon Oscillation Spectroscopic Survey (BOSS; Ahn et al. 2014), constraining the effects of feedback and finding tensions with cosmological simulations. The signal-to-noise ratio (S/N) in measurements like those of Schaan et al. (2021) and Amodeo et al. (2021) will increase rapidly, as forecasted in Battaglia et al. (2017), and will make the tSZ and kSZ effects even more powerful probes of astrophysical processes related to galaxy formation and evolution (Ade et al. 2019; Abazajian et al. 2019; Battaglia et al. 2019). With such a high S/N, correctly modeling details of the sample becomes very important, and this is the crucial question that this study addresses.

Since it is not possible to directly observe galaxy evolution over cosmological timescales, simulations offer the only method of tracking and characterizing the physical properties and locations of baryons in the galaxies (including the CGM and even farther into the IGM) in controlled environments at multiple redshifts, which can act as predictors for how we believe real galaxies behave. In this study we have developed methods using cosmological simulations to study the

thermodynamic properties of galaxy halos, thus constraining the baryonic processes that are responsible for making star formation inefficient and affecting the galaxy's evolutionary track.

In the current era of high signal-to-noise SZ observations, modeling the halo samples becomes a major and important uncertainty to understand. Not taking certain properties of the sample into account in the model could bias the way in which the observations are interpreted.

In this study we use simulated halo samples from the IllustrisTNG (The Next Generation) simulation (Marinacci et al. 2018; Naiman et al. 2018; Pillepich et al. 2018; Springel et al. 2018; Nelson et al. 2018, 2019) to test various components of modeling an example of an observed halo population. We model a subsample of the CMASS sample from the BOSS survey, Data Release 10 (Ahn et al. 2014), of which SZ profile measurements have already been made (Amodeo et al. 2021; Schaan et al. 2021). This subsample is chosen to contain galaxies in the region covered by ACT, that are mainly central galaxies of group-sized halos, selected with a halo mass lower than $10^{14} M_{\odot}$. The sample is spectroscopic, with a range of redshifts $0.4 < z < 0.7$ and median redshift $z = 0.55$. This is just one example of a possible observed halo sample, but the results are expected to be applicable to modeling other samples as well.

In this paper we summarize the simulations used in Section 2.1. We discuss how we extract the simulated halo information and construct thermodynamic profiles in Section 2.2, the modeling uncertainties explored in Section 2.3, and the fitting procedure in Section 2.4. We discuss our methods of modeling the CMASS sample in Section 3.1, and projecting the simulated three-dimensional profiles into a two-dimensional, observing space in Section 3.2. Finally, we discuss the results of fitting in 3D in Section 4.1 and 2D in Section 4.2, and how the modeling uncertainties affect the shape of the profiles and fitting parameters in each.

## 2. Theoretical Profiles

We begin by describing our methods for producing and analyzing simulated three-dimensional profiles, which will also be used for computing and analyzing the two-dimensional profiles discussed in Section 3.

### 2.1. Simulations

We have developed methods using the Illustris and IllustrisTNG simulations to study halo properties. Illustris (Vogelsberger et al. 2014) is the original large-scale hydrodynamical simulation of galaxy formation, and provides the foundation on which its successor TNG (Springel et al. 2018) was built. Illustris includes models for gas cooling, stellar evolution, various forms of feedback (e.g., from supernovae and AGN), and many other physical processes to produce galaxy and cluster samples that matched observed trends and scaling relations well (Vogelsberger et al. 2014). However, Illustris did struggle to reproduce other trends, such as low gas mass fractions in halos, discussed in detail in Genel et al. (2014); Vogelsberger et al. (2014); Nelson et al. (2015). Thus, the TNG simulations were created to address these discrepancies, along with including numerical advancements to the overall simulation code. TNG includes the same basic physical processes as the Illustris simulations, and additionally includes a new implementation of growth and feedback from supermassive black
<figure><img src="image_1.png" /></figure>

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Mean profiles of 3D gas density (left) and thermal pressure (right) of Illustris (blue), TNG100-3 (red), and TNG100-2 (black) as functions of radius at redshift $z = 0.55$ for halo masses $12 \le \log_{10}(M_{\text{h}}/\text{M}_\odot) \le 13$. The bands show the $\pm 1\sigma$ distribution of the values in each radial bin.</figcaption></figure>

holes (SMBHs), galactic winds, and magnetic fields based on ideal magnetohydrodynamics (Weinberger et al. 2017; Marinacci et al. 2018; Pillepich et al. 2018).<fnref n="5" />

Both Illustris and TNG simulations have varying box sizes and resolutions. For Illustris, the smallest volume and lowest resolution run is Illustris-3, which has a box size of $106.5^3$ Mpc³ and tracks $2 \times 455^3$ gas and dark matter particles. For TNG the smallest volume and lowest resolution run currently available is TNG100-3, with a box size of $110.7^3$ Mpc³ and $2 \times 455^3$ tracked particles. The TNG simulation with the next lowest level of resolution for the same box size is TNG100-2, and has $2 \times 910^3$ resolution elements. In the following section, we compare the thermodynamic profiles of halos generated from these sets of simulations and we focus on TNG100-3 for the remainder of the paper.

## 2.2. Generating 3D Profiles

We have identified simulated halos with certain properties of interest and created thermodynamic radial profiles using the publicly available stacking code repository *Illstack*.<fnref n="6" /> An already existing public TNG repository *illustris-python*<fnref n="7" /> allows for the extraction of particle information from each snapshot, and offers access to the group catalogs of the simulations which include positions, masses, radii, and other quantities of the halos defined by the “friends-of-friends” algorithms. In theory, one could study the simulated halo properties only using the information provided in the group catalog, however the study would be limited by the predetermined cuts and boundaries defined by the TNG halo finding algorithms in post-processing and would lack particles beyond a certain radius. *Illstack* expands on the framework of *illustris-python* by storing the quantities of interest of every particle in the entire snapshot along with the halo positions of the group catalog, then linking the particles to the halos out to arbitrarily high radii. We do not perform cuts in the phase space of the simulated gas particles; rather, we select the halos for which we compute profiles by their group properties (mass and/or color) discussed in Section 2.3. The halos are then split into several radial bins, and the quantities are stacked

in each of these bins to create 3D radial profiles. We define the radial bins to extend from $1 \times 10^{-4} - 1 \times 10^2$ Mpc split into 25 bins in logspace, resulting in radial bin size $\sim 0.55$ Mpc. *Illstack* also accounts for periodic boundary conditions to ensure that all halo and particle information is utilized even if a halo is located near the edge of the simulated box.

Examples of the output from *Illstack* are shown in Figure 1. In particular, the profiles display the properties of pressure and gas density to be directly related to the observable tSZ and kSZ effects, respectively. This figure shows the density and pressure profiles derived from halos of the same mass range ($12 \le \log_{10}(M/\text{M}_\odot) \le 13$, with M as the halo mass of the group catalog) and same redshift ($z=0.55$) of the three different simulations described previously in Section 2.1. The solid lines show the median of each radial bin, and the bands show the $\pm 1\sigma$ distribution of the values within each radial bin. Differences among all of these simulations can be seen in the profiles, such as higher density values at inner regions and lower density values from middle regions from Illustris compared to TNG, and similar differences in the pressure profiles, which could be a result of the different prescription of AGN feedback among the simulations. Since the shortcomings of Illustris are well known and documented, we focus our analyses on TNG for the remainder of the paper.

Furthermore, given the mass ranges in which we are interested we do not see significant differences between the resolution levels of the TNG simulations, and thus do not require the higher resolution. Therefore we can be efficient with our computing resources and use the smallest box currently available with the lowest resolution, TNG100-3.

## 2.3. Modeling Uncertainties

Once the simulated profiles are computed, we can further divide them to study various halo sample modeling uncertainties using *Illstack*. It has the ability to select and make profiles for halos of different types (e.g., within a desired mass range, different redshifts, and colors). We test different options for modeling halo populations including splitting by mass and color, adding a two-halo term to the fitting model, and matching the shape of an observed mass distribution.

### 2.3.1. *Mass-splitting Definitions*

Relating the stellar contents of a galaxy to the dark matter halo in which it resides has been the subject of many studies in

<footnote n="5">⁵ We note that TNG uses the Planck cosmology from Planck Collaboration et al. (2016b), while Illustris adopted the Wilkinson Microwave Anisotropy Probe cosmology from Hinshaw et al. (2013).</footnote>

<footnote n="6">⁶ https://github.com/marcelo-alvarez/illstack</footnote>

<footnote n="7">⁷ https://github.com/illustristng/illustris_python</footnote>
the past decade (e.g., Leauthaud et al. 2011, 2012; Behroozi et al. 2010, 2013), and an understanding of the topic is essential for understanding how dark matter plays a role in galaxy formation and evolution. This widely studied relation is called the stellar–halo mass relation (SHMR), and it attempts to map the stellar mass of a galaxy, usually estimated through photometric or spectroscopic surveys, to the underlying dark matter halo mass, usually calculated through lensing, kinematics, or abundance matching. The functional form of the relation is debated, as well as its redshift dependence (Behroozi et al. 2013).

As described in Behroozi et al. (2010), measurements of the stellar mass of galaxies introduce significantly more systematic uncertainties than measurements of the halo mass, which could be a reason the computed relations vary among studies using different surveys and methods (Shankar et al. 2014). While the uncertainty in stellar mass function is more troublesome, there are still systematics in the derivation of halo masses, described in more detail by Behroozi et al. (2010). With all of the uncertainties in both sides of the relation, the resulting relation is therefore quite uncertain.

We use the TNG simulations in this study so that we have values for both stellar and halo masses. Illstack has the ability to select for halos by either mass option, so we derive stacked profiles for two separate populations based on mass type. This allows us to explore the differences in profiles of these samples, along with effects of stacking on samples with different mass distributions discussed further in Section 2.3.4, we show how important it is to correctly model the sample, including the SHMR. Our exploration is specific to TNG, but it does provide a range of uncertainty that could arise from the SHMR.

Aside from looking at differences in the samples based on the mass type selection, we also perform a study on how the profiles change as a function of mass (both stellar and halo). We derive the TNG samples to match the limits of the CMASS sample analyzed in Amodeo et al. (2021) and discussed further in Section 3.1). We further explore any trends in mass by splitting this sample into four smaller mass bins, discussed further in Section 4.1.2.

### 2.3.2. Color-splitting Definitions

One of the previously mentioned enigmas of CGM theory is the cause of galaxy quenching. The theoretical picture for quenched galaxies is that they are mostly red in color, as a result of little or no star formation. With Illstack we can separate the simulated halos into a red population and study the differences in the equivalent profiles of the entire sample that also includes active, bluer halos. Since the group catalogs of the Illustris simulations do not provide color information for the main halo, we use the colors from the most massive subhalo within the host halo of the TNG group catalog as representative of the whole halo. We derive subpopulations within each of the mass bins based on color using the TNG color cut in Nelson et al. (2018), which is defined as the difference in magnitudes $g-r > 0.6$. Using this color cut, we split the populations into one without a color cut (“tot” for total), and one with a red color cut (“red”).

### 2.3.3. Two-halo Term

Since Illstack can make profiles out to arbitrarily high radii, it can be used to study other theoretical aspects of the interaction of the CGM and IGM, i.e., the two-halo term. As

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Density profile (top) and pressure profile (bottom) of a halo with mass $2 \times 10^{13} M_{\odot}$ and redshift $z = 0.57$ as functions of scaled radius. The profiles are split into the components of the profiles as a function of scaled radius for the one-halo (red, dashed) and two-halo (blue, dashed) terms, using the two-halo term model of Vikram et al. (2017). The profile combining the terms is shown in purple.</figcaption></figure>

can be seen in Figure 1 the profiles are generally decreasing functions of radius but there is a point of inflection near the outer part of the profile. This part of the profile has been called the *two-halo term* and is due to contributions from neighboring halos and the overall enhancement of density and pressure from gravitational clustering of gas that has not yet virialized (Cooray & Sheth 2002; Vikram et al. 2017; Hill et al. 2018).

The part of the profile associated with virialized gas in the main halo is called the *one-halo term* and represents the signal only from the main halo. A model for the two-halo term was first described by Vikram et al. (2017), but it is largely neglected in SZ measurement literature. It has been shown that the two-halo term contribution to measured halo signals can be significant, especially around lower mass halos (Vikram et al. 2017; Hill et al. 2018), and is therefore important to include when modeling observed signals. We show the contributions of the different terms as a function of radius in Figure 2. The figure shows a density profile (top panel) and pressure profile (bottom panel) for a halo of mass $2 \times 10^{13} M_{\odot}$ at redshift $z = 0.57$ as functions of radius scaled by the virial radius, $r_{200c}$, defined as the radius of the sphere whose mean density is 200 times the critical density of the universe at redshift $z$. The red and blue dashed lines correspond to the one-halo and two-halo
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Left panel: Comparison of the SHMR used for calculating corresponding halo masses for CMASS stellar mass sample (Kravtsov et al. 2018) and SHMR of TNG. The SHMR relation compared to the simulations sample is clearly different, showing the importance of this modeling uncertainty, which we explore through modeling the CMASS sample two ways, via halo masses and stellar masses. Right panel: Histograms of the CMASS sample analyzed in Amodeo et al. (2021) shown in gray and the TNG mh sample in red. The blue vertical lines show the weights derived by Amodeo et al. (2021). This comparison shows that while the limits of the TNG sample are chosen to match those of the CMASS sample, the distributions are not the same shapes, resulting in the importance of considering distribution matching.</figcaption></figure>

terms, respectively, and the solid purple line shows the total profile combining the terms. It can be seen that for this halo, the two-halo term begins to significantly contribute to the profile around 1 virial radius. The amplitude and radius at which the two-halo term begins to contribute varies according to the mass and redshift of a halo, discussed further in Vikram et al. (2017).

Here we include the model for the two-halo term for both density and pressure profiles derived in Vikram et al. (2017), multiplied by an amplitude $A_{2h}$ that we allow as a free parameter when fitting the profiles as in Amodeo et al. (2021), discussed further in Section 2.4.

### 2.3.4. Matching the Shape of the Mass Distribution

The distribution of the halos in the simulation does not match the distribution of the data as seen in the right panel of Figure 3, therefore a straight average of the halos in the mass range could be biased. We test the importance of matching the shape of the mass distribution by weighting each halo appropriately to ensure that the contribution of each mass bin of the simulated sample matches the contribution of the corresponding mass bin of the observed sample.

We use the same weights as those used in Amodeo et al. (2021). In more detail, the weights are computed so the area of the histogram of the CMASS sample (shown in the right panel of Figure 3) is equal to 1 (i.e., $\sum_i w_i \Delta m = 1$, where $w_i$ are the weights of each mass bin and $\Delta m$ is the width of the bin). The bin edges are shown by the blue vertical lines in the right panel of Figure 3, so a different weight is applied to the halos falling in each of these bins, according to their mass. This is especially important for the TNG halos that fall into the first two mass bins, since the shape of the distribution within these bins significantly differs from the observed CMASS distribution. We show this can cause considerable bias in the inferences of the profiles' parameters.

The specific method for the weighting is we assign each halo a weight according to its mass, then do a weighted average for all halos within each radial bin to get the *matched* stack. The equation for the matched profile value in each radial bin would be $val_m = \frac{\sum val_m^* w}{\sum w}$. In contrast, the unmatched (unweighted) average

would be $val_{um} = \frac{vals}{N}$ ($N$ is number of halos in the bin), with all halos contributing an equal amount. This process is only relevant for the calculation of the raw 3D radial profiles; we are not weighting the observed 2D profiles, rather we do the weighing and stacking of the 3D profiles then project the average (weighted and unweighted) into the 2D SZ observables.

As the tSZ signal is proportional to $M^{5/3}$ and the kSZ signal is linearly dependent on mass (see Equations (1) and (2)), we expect the distribution matching to have more of an effect on the pressure profiles than density. However, we provide both matched and unmatched results for both quantities. This kind of weighting could be done for the redshift distribution as well, but as described in Amodeo et al. (2021), the CMASS redshift distribution is peaked around the median and the resulting profiles with redshift weighting do not significantly differ from profiles without redshift weighting. Therefore, we simply use the median redshift of the CMASS distribution.

We also explore the effects of the different mass type-selected samples to further study the importance of correctly modeling the SHMR. In deriving the mass bins within each type to perform the distribution matching, we chose limits to align with the sample analyzed in Amodeo et al. (2021). Within each of the limits, we further split into smaller bins and assigned weights based on a normalized distribution (same procedure and weights as in Amodeo et al. 2021), shown in the right panel of Figure 3. Since the CMASS sample was chosen spectroscopically by stellar mass, we need to convert these stellar mass bins into corresponding halo mass bins to calculate the quantities that require halo mass information. Following the same procedure as Amodeo et al. (2021), we use the SHMR of Kravtsov et al. (2018) to perform this conversion. Since TNG has not been tuned to match these sorts of observed relations, the Kravtsov relation does not match the TNG SHMR (see the left panel of Figure 3). By using the Kravtsov relation as the basis of converting our stellar mass weighting bins to halo masses, we can see the effects of using the same weights calculated for a stellar mass distribution (like the CMASS sample) on a sample that has been selected for halo mass. In
other words, we are using the same weights (computed using the SHMR) on halos that are selected by different mass types.

The TNG sample shown in the left panel of Figure 3 is halo mass selected for halos with masses $11 \le \log_{10}(M/M_c) \le 14$. As can be seen in the figure, the slope of the SHMR for TNG is much steeper than the relation of Kravtsov et al. (2018).

## 2.4. Fitting Procedure of Simulated Profiles

We fit the profiles computed for the different samples (red versus total, matched versus unmatched) using the simulated data ($\vec{d}$) and errors ($\vec{\Sigma}$), given by the $\pm 1\sigma$ distribution of profiles in each radial bin shown in Figure 1. We assume the likelihood, $L(\vec{d})$, to have the form

$$ \ln L(\vec{d}) = -\frac{1}{2}[\vec{d} - \vec{\mu}(\vec{\theta})]^T \vec{\Sigma}^{-1} [\vec{d} - \vec{\mu}(\vec{\theta})], \quad (3) $$

where $\vec{\mu}(\vec{\theta})$ is the model evaluated with parameters $\vec{\theta}$. The posterior, $P(\vec{\theta})$, is written in terms of the likelihood as

$$ P(\vec{\theta}) \propto L(\vec{d}) \Pr(\vec{\theta}), \quad (4) $$

where $\Pr(\vec{\theta})$ are the priors on the parameters $\vec{\theta}$. We maximize the posterior probability functions using the Markov chain Monte Carlo (MCMC) calculation package `emcee` (Foreman-Mackey et al. 2013) to quantify the different effects and biases of model selection (one-halo versus two-halo term GNFW model, matched versus unmatched) and sample selection (red versus total, stellar mass versus halo mass).

For the models $\vec{\mu}(\vec{\theta})$ in Equations (3) and (4) we use a generalized Navarro–Frenk–White (GNFW) profile (Zhao 1996, see also Hernquist 1990; Navarro et al. 1997), as it is a simple parametric model that captures the shape of the density and pressure profiles (Nagai et al. 2007).

The GNFW gas density profile has the functional form

$$ \begin{aligned} \frac{\rho(x)}{\rho_{cr}(z)} &= f_b \rho_0(x/x_{c,k})^{\gamma_k} [1 + (x/x_{c,k})^{\alpha_k}]^{-\frac{\beta_k - \gamma_k}{\alpha_k}}, \\ \rho_{cr}(z) &= \frac{3H_0^2}{8\pi G} [\Omega_m(1+z)^3 + \Omega_{\Lambda}], \end{aligned} \quad (5) $$

where scaled radius $x \equiv r/r_{200c}$, $f_b$ is the baryon fraction $\Omega_b/\Omega_m$, $\rho_{cr}(z)$ is the critical density of the universe at redshift $z$, $H_0$ is the Hubble constant, and $G$ is the gravitational constant.

The GNFW thermal pressure profile has the functional form

$$ \begin{aligned} \frac{P(x)}{P_{200c}} &= P_0(x/x_{c,t})^{\gamma_t} [1 + (x/x_{c,t})^{\alpha_t}]^{-\beta_t}, \\ P_{200c} &= \frac{200GM_{200c} \rho_{cr}(z) f_b}{2r_{200c}}, \end{aligned} \quad (6) $$

where $M_{200c}$ is the mass of the sphere whose mean density is 200 times the critical density of the universe at a redshift $z$. In both Equations (5) and (6), $x_c$ is the core radius, $\alpha$ is the slope of the GNFW profile at $x \sim 1$, $\beta$ is the outer slope of the profile at $x \gg 1$, and $\gamma$ is the inner slope of the profile at $x \ll 1$ (Zhao 1996).

For the gas density profiles (Equation (5)) we perform fits for parameters $\rho_0$, $\alpha_k$, and $\beta_k$, keeping parameters $x_{c,k}$ and $\gamma_k$ fixed due to degeneracies, described in Battaglia (2016). The values we adopt for these fixed parameters along with the flat priors on the free parameters ($\Pr(\vec{\theta})$ in Equation (4)) are shown in Table 1. The fits are not sensitive to choice of initial states, but

<caption>**Table 1**</caption>

Values for the Fixed Parameters and Free Parameters in the GNFW Profile Fits

<table><thead><tr><th>Pressure</th><th>Density</th></tr></thead><tbody><tr><td colspan="2">Fixed</td></tr><tr><td>$\alpha_t = 1.0$</td><td>$x_{c,k} = 0.5$</td></tr><tr><td>$\gamma_t = -0.3$</td><td>$\gamma_k = -0.2$</td></tr><tr><td colspan="2">Free</td></tr><tr><td>$P_0, x_{c,t}, \beta_t, A_{t2h}$</td><td>$\rho_0, \alpha_k, \beta_k, A_{k2h}$</td></tr><tr><td>$P_0$: (0.1,30)</td><td>$\log_{10} \rho_0$: (1.0, 6.0)</td></tr><tr><td>$x_{c,t}$: (0.01,1.0)</td><td>$\alpha_k$: (0.1,6.0)</td></tr><tr><td>$\beta_t$: (1.0,10.0)</td><td>$\beta_k$: (1.0,10.0)</td></tr><tr><td>$A_{t2h}$: (0.01,5.0)</td><td>$A_{k2h}$: (0.01,5.0)</td></tr></tbody></table>

Note. $A_{2h}$ is the amplitude of the two-halo term. Priors for the fits are given for the free parameters.

we use the fitting formulas derived in Battaglia (2016) as the starting points for each of the MCMC chains.

For the thermal pressure profiles (Equation (6)) we perform fits for parameters $P_0$, $x_{c,t}$, and $\beta_t$, keeping parameters $\alpha_t$ and $\gamma_t$ fixed as described in Battaglia et al. (2012a). Similarly to the treatment of the density fits, we use the fitting formulas derived in Battaglia et al. (2012a) as the starting points for each of the MCMC chains.

As in Amodeo et al. (2021), we add in the contributions of neighboring halos with the addition of a two-halo term multiplied by an amplitude, $A_{2h}$, which is a free parameter in the fits, such that the final forms of the fitting equations are

$$ \begin{aligned} P_{\text{GNFW1h}} &= P_{1h}, \\ P_{\text{GNFW}} &= P_{1h} + A_{t2h} P_{2h}, \end{aligned} \quad (7) $$

where $P_{\text{GNFW1h}}$ is the model used when just taking the one-halo term into account, and $P_{\text{GNFW}}$ is the model used for taking both one-halo and two-halo terms into the fit. The fitting equations for the different models of density similarly take this form.

We run several independent MCMC chains for each sample until the Gelman–Rubin statistic (Gelman & Rubin 1992) reaches values of $\lesssim 1.1$. We combine the independent chains and derive marginalized estimates for each parameter using the median $\pm 1\sigma$ as the error estimates. We select the model resulting in the minimum $\chi^2$ as our best fit, shown as the fits in Figure 4.

## 3. Observational Profiles

In order to understand how the above modeling choices described in Section 2.3 impact the CGM properties that we will infer from observations we need to forward model the simulated profiles from TNG into an observable space, e.g., projected SZ profiles. Here we choose to model a subsample of CMASS, as recent SZ profile measurements of this sample have been made (Amodeo et al. 2021; Schaan et al. 2021). The general outline of this procedure is to make the modeling selections within the 3D simulated TNG halos, project the computed profiles to how they would actually be observed in the 2D observing space of SZ effect measurements, and fit the profiles with the same process as described in Section 2.4.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Simulated profiles of density (left) and pressure (right) showing the populations fit with different mass selections, color selections, distribution-matching options, and fitting models. The solid lines are the TNG profiles for the population specified by the legend, the dashed lines are the best fit from the GNFW model, and the dotted lines are the best fit from the GNFW1h model. In the bottom panels, the profiles and fits are normalized by a fiducial TNG profile for sample ms-red, m (solid red curve) denoted as $\rho_{\text{gas}}^*$ and $P_{\text{th}}^*$. This figure shows that selections by color and mass type (seen in the left panels) and the distribution matching (seen in the right panels) have varying levels of significance on the profiles. Furthermore, this figure shows the inclusion of a two-halo term to the fitting model provides a closer match to the profile than a model that only includes the one-halo term.</figcaption></figure>

### 3.1. Modeling the Sample Selection

For our fiducial model, we design our simulated halo sample and fitting model to best reflect the properties of the observed CMASS halo sample.

We use the halos of a single TNG snapshot at $z=0.55$ to match the median redshift of the CMASS sample, $z=0.55$ (Ahn et al. 2014). We use Illstack to select for stellar masses to match the sample analyzed in Amodeo et al. (2021) shown in the right panel of Figure 3, and further separate into finer mass bins within this range to weight the TNG sample to match the observed mass distribution. More specifically, the stellar mass limits of the sample are $10.71 \le M \le 11.72$, and the halo mass limits of the sample are $12.12 \le M \le 13.98$, where $M = \log_{10}(M_x/M_{\odot})$ and $M_x$ is either stellar or halo mass. CMASS galaxies are large red galaxies (LRGs; Padmanabhan et al. 2007), so we use Illstack to select for red halos as defined by the TNG color cut described in Section 2.3.2. Lastly, we include a two-halo term to the model as described in Section 2.4. Table 2 summarizes the labels for each modeling uncertainty.

### 3.2. Generating 2D Profiles

The signals of the kSZ and tSZ effects are integrated along the LOS, so we need to project our simulated 3D profiles in the same way. We use the repository Mop-c-GT,<fnref n="8" /> (Model-to-observable projection code for Galaxy Thermodynamics), introduced in Amodeo et al. (2021). This repository inputs 3D gas density and pressure radial profiles and outputs 2D

<table><caption>**Table 2**
Explanation of the Labels Used to Identify the Choices Included for the Halo Sample and Fitting Model, See Section 3.1</caption><thead><tr><th colspan="2">Sample Selection</th></tr></thead><tbody><tr><td>"ms" = Stellar mass-selected</td><td>"red" = Color-selected, g - r > 0.6</td></tr><tr><td>"mh" = Halo mass-selected</td><td>"tot" = No color selection</td></tr></tbody><thead><tr><th colspan="2">Fitting Model</th></tr></thead><tbody><tr><td>"m" = CMASS mass distribution matched</td><td>"GNFW" = With two-halo term</td></tr><tr><td>"um" = Unmatched</td><td>"GNFW1h" = No two-halo term</td></tr></tbody></table>

profiles of observable quantities through LOS projection. Specifically, it produces profiles of the temperature shifts in CMB signal due to the kSZ and tSZ effects, which we derive from the simulated gas density and pressure profiles, respectively. For the kSZ signal, we assume that the rms of the peculiar velocities projected along the LOS is 313 km s⁻¹, as predicted by the linear theory at $z=0.55$ and adopted in Schaan et al. (2021). This value only impacts the amplitude of the signal and is a systematic uncertainty in the velocity reconstruction, which is beyond the scope of this analysis. Mop-c-GT also allows for forecasting of signals received by individual instruments with instrumental beam convolution, and an aperture photometry filter is also used, described in more detail in Amodeo et al. (2021) and Schaan et al. (2021). Here we convolve the profiles with a Gaussian beam of 1.4' at frequency 150 GHz, corresponding to the forecasted Simons Observatory (SO) experimental setup described in Ade et al. (2019). We are using a method to create stacked profiles of

<footnote n="8">⁸ https://github.com/samodeo/Mop-c-GT</footnote>
multiple hundreds of halos. The halos could have differing
individual orientations, but by stacking together we are able to
assume a spherically symmetric profile. Thus, changing the
direction of projection should not yield different results (see
Battaglia et al. 2012b). To further demonstrate that additional
contributions along the LOS are completely subdominant, we
performed a test extending the projected profiles out to ~50
Mpc to see whether the contributions to the profile from larger
radii affect the results of the projection process. We find that
the difference in the projected SZ profiles extending to 10 Mpc
and extending to 50 Mpc is completely negligible, <<0.1%,
which highlights the benefit of using the aperture photometry
filter.

Similar analyses and profiles can be made for other CMB
experiments, such as CMB-S4, and the subsequent error bars
will change with the different sensitivities, frequency channels,
and beams (e.g., see Section 1.4.2.1 in the CMB-S4 Science
Case and Reference Design document, Abazajian et al. 2019).

3.3. Fitting Procedures for Observational Profiles

We compute the 2D observational profile using the parameters $\vec{\theta}$ from the best fit of the 3D theoretical profile using MOP-c-GT, with the same GNFW models $\vec{\mu}(\vec{\theta})$ as defined in Section 2.4. We use this profile as the data, $\vec{d}$, to calculate the likelihood of the same form as Equation (3). However, instead of using the distribution of the 3D profiles as the errors $\vec{\Sigma}$ we use the covariance matrix of forecasted errors for the SO experiment of Battaglia et al. (2017), in which the authors use a semi-analytical foreground model that includes contributions from the cosmic infrared background, primary CMB fluctuations, extragalactic radio emission, and galactic cirrus.

We fit for the same parameters ($\vec{\theta}$) for each kind of profile
and use the same priors, Pr($\vec{\theta}$), previously described in
Section 2.4 and shown in Table 1.

4. Results

4.1. 3D Simulations

Here we show the results of fitting the 3D simulated profiles,
focusing on the different modeling choices and any trends in
mass. The parameters of the fits for each sample and model are
shown in Table 3.

4.1.1. Modeling Choices

As a qualitative example to demonstrate the effects of
different modeling procedures, we show fits of the 3D profiles
in Figure 4. The profiles computed from TNG are shown in
solid lines, while the best GNFW and GNFW1h fits are shown
by dashed and dotted lines, respectively. To be able to see
differences more clearly, all of the profiles and fits have been
normalized by the fiducial TNG profile for the ms-red, m
sample, denoted as $\rho_{\text{gas}}^*$ and $P_{\text{th}}^*$ in the bottom panels of the
figure. We explore different modeling choices in each column,
though the trends are similar for both pressure and density.
Both columns show profiles of differing mass selections;
additionally, in the left column we show profiles differing in
color selection, and in the right column we show profiles
differing in distribution-matching options.

First, the left column of Figure 4 shows the differences in
sample selection, including mass type and color, for the
distribution-matched density profiles. While the inner regions

of the profiles for the ms-tot and mh-tot samples have higher
values than the other profiles, it can be seen that for the rest of
the profiles, the ms samples tend to have higher values than the
mh samples. Beyond the inner part of the profiles, the color
selection does not appear to have a significant effect: for the mh
samples the red color-selected sample has higher values than
the corresponding tot sample (compare solid purple curve with
solid turquoise curve), and within the ms samples the profiles
of the red and tot samples are nearly identical (compare the
solid red and solid black curves).

Next, the right column of Figure 4 shows the differences in distribution-matching options and mass selection of the pressure profiles. A separation between the matched and unmatched profiles can be seen in both panels, with the unmatched profiles having lower values (≤0.7 times the fiducial profile $P_{\text{th}}^*$, which is matched) than their matched counterparts. This is expected since the matched and unmatched samples are compiled from different halo distributions (Figure 3). Mass distribution matching is expected to have a larger effect on the pressure profiles than the density profiles, due to the pressure's dependence on $M^{5/3}$. Differences are seen among the density profiles as well, although to a lesser extent. Similarly to the density profiles of the left column, within each of the weighting groupings the profiles of the ms samples tend to have higher values than those of the mh samples.

Lastly, Figure 4 shows how the fits vary for the different
fitting models, GNFW versus GNFW1h. An obvious feature in
both bottom panels of the figure is the spike in values around
1 Mpc, although more significant for the GNFW1h model
(dotted) and more significant for pressure than density. This
indicates the fits' inabilities to accurately reflect the shape of
the outer profiles, and in this case the best GNFW1h models
return higher values than the actual profiles at these radii by
factors of 2 or 3. It can be seen by eye that the GNFW model
provides better fits to the profiles, as the spikes are less
significant and they more closely align with the actual profiles
at other radii.

4.1.2. Mass Dependence

As described in Section 2.3.1, we split the simulated halos
into mass bins to observe any trends in our models. In Figure 5
we show the marginalized estimates for parameters $\rho_0$ and $\beta_t$,
resulting from the MCMC chains. The points in each mass bin
show the median of all the samples with the error bars showing
the $1\sigma$ range. The top and bottom axes show the stellar mass
and halo mass bins used to select the halos, respectively. We
note that the top and bottom axes are not related by an SHMR
but they just show the mass bins further divided, therefore one
should not compare the ms versus mh (red versus blue) values.

Some of the parameters have flat trends with mass, but others
have relations like those shown in Figure 5. The top panel
shows that the normalization factor for the GNFW density
profile has a negative trend scaling with mass, and the bottom
panel shows that $\beta_t$ for the GNFW1h pressure profile has a
positive trend with mass. The CMASS sample spans a large
range in masses and such a range is probably not unique to
CMASS. Figure 5 shows that some parameters show mass
dependencies while others do not. Capturing such mass trends
is important when modeling observable samples for cross
correlations, and is a motivating reason for why we looked into
matching the sample to the observed mass distribution.
<table><caption>**Table 3**

Parameters from 3D Fits for Each Sample and Model</caption><thead><tr><th></th><th colspan="4">Density</th><th colspan="4">Model GNFW-matched</th></tr><tr><th></th><th>mstar-tot</th><th>mstar-red</th><th>mh-tot</th><th>mh-red</th><th></th><th>mstar-tot</th><th>mstar-red</th><th>mh-tot</th><th>mh-red</th></tr></thead><tbody><tr><td>log<sub>10</sub> ρ<sub>0</sub></td><td>3.67,<br>3.04<sup>+0.60</sup><sub>-0.28</sub></td><td>3.28,<br>2.92<sup>+0.41</sup><sub>-0.21</sub></td><td>4.34,<br>3.22<sup>+0.84</sup><sub>-0.42</sub></td><td>3.38,<br>2.88<sup>+0.46</sup><sub>-0.23</sub></td><td>P<sub>0</sub></td><td>4.18, 8.86<sup>+8.25</sup><sub>-4.17</sub></td><td>4.03, 8.66<sup>+7.89</sup><sub>-4.02</sub></td><td>2.74, 6.94<sup>+8.50</sup><sub>-3.62</sub></td><td>2.86, 7.45<sup>+8.77</sup><sub>-3.96</sub></td></tr><tr><td>α<sub>k</sub></td><td>0.64,<br>1.49<sup>+1.66</sup><sub>-0.79</sub></td><td>0.80,<br>1.65<sup>+1.62</sup><sub>-0.84</sub></td><td>0.43,<br>1.10<sup>+1.19</sup><sub>-0.56</sub></td><td>0.68,<br>1.67<sup>+1.78</sup><sub>-0.88</sub></td><td>x<sub>c,t</sub></td><td>0.97, 0.59<sup>+0.27</sup><sub>-0.30</sub></td><td>1.00, 0.60<sup>+0.27</sup><sub>-0.30</sub></td><td>0.98, 0.51<sup>+0.31</sup><sub>-0.28</sub></td><td>0.93, 0.53<sup>+0.30</sup><sub>-0.29</sub></td></tr><tr><td>β<sub>k</sub></td><td>3.37,<br>4.31<sup>+2.01</sup><sub>-1.06</sub></td><td>3.20,<br>4.10<sup>+1.99</sup><sub>-1.01</sub></td><td>3.44,<br>4.76<sup>+2.42</sup><sub>-1.31</sub></td><td>3.30,<br>4.90<sup>+2.67</sup><sub>-1.55</sub></td><td>β<sub>t</sub></td><td>6.09, 5.46<sup>+2.42</sup><sub>-2.07</sub></td><td>6.14, 5.43<sup>+2.36</sup><sub>-2.00</sub></td><td>6.31, 5.65<sup>+2.65</sup><sub>-2.40</sub></td><td>5.80, 5.67<sup>+2.59</sup><sub>-2.31</sub></td></tr><tr><td>A<sub>k2h</sub></td><td>1.40,<br>1.60<sup>+0.43</sup><sub>-0.36</sub></td><td>1.37,<br>1.57<sup>+0.43</sup><sub>-0.36</sub></td><td>1.31,<br>1.50<sup>+0.43</sup><sub>-0.34</sub></td><td>1.34,<br>1.53<sup>+0.41</sup><sub>-0.34</sub></td><td>A<sub>t2h</sub></td><td>0.57, 0.95<sup>+0.82</sup><sub>-0.47</sub></td><td>0.53, 0.87<sup>+0.72</sup><sub>-0.42</sub></td><td>0.51, 0.98<sup>+1.01</sup><sub>-0.52</sub></td><td>0.51, 0.91<sup>+0.82</sup><sub>-0.46</sub></td></tr></tbody></table><table><thead><tr><th></th><th colspan="4">Density</th><th colspan="4">Model GNFW-unmatched</th></tr><tr><th></th><th>mstar-tot</th><th>mstar-red</th><th>mh-tot</th><th>mh-red</th><th></th><th>mstar-tot</th><th>mstar-red</th><th>mh-tot</th><th>mh-red</th></tr></thead><tbody><tr><td>log<sub>10</sub> ρ<sub>0</sub></td><td>4.11,<br>3.25<sup>+0.89</sup><sub>-0.42</sub></td><td>3.38,<br>2.98<sup>+0.59</sup><sub>-0.26</sub></td><td>5.68,<br>3.62<sup>+1.04</sup><sub>-0.58</sub></td><td>3.48,<br>2.91<sup>+0.57</sup><sub>-0.26</sub></td><td>P<sub>0</sub></td><td>3.95, 9.31<sup>+8.70</sup><sub>-4.53</sub></td><td>4.91, 10.70<sup>+8.88</sup><sub>-5.05</sub></td><td>3.84, 10.16<sup>+9.56</sup><sub>-5.25</sub></td><td>4.47, 10.85<sup>+9.59</sup><sub>-5.56</sub></td></tr><tr><td>α<sub>k</sub></td><td>0.49,<br>1.02<sup>+1.09</sup><sub>-0.51</sub></td><td>0.68,<br>1.24<sup>+1.16</sup><sub>-0.62</sub></td><td>0.30,<br>0.83<sup>+0.77</sup><sub>-0.38</sub></td><td>0.61,<br>1.42<sup>+1.48</sup><sub>-0.75</sub></td><td>x<sub>c,t</sub></td><td>0.98, 0.59<sup>+0.28</sup><sub>-0.31</sub></td><td>0.98, 0.62<sup>+0.26</sup><sub>-0.31</sub></td><td>0.98, 0.59<sup>+0.28</sup><sub>-0.30</sub></td><td>0.95, 0.60<sup>+0.27</sup><sub>-0.30</sub></td></tr><tr><td>β<sub>k</sub></td><td>3.30,<br>3.85<sup>+1.32</sup><sub>-0.79</sub></td><td>2.91,<br>3.41<sup>+1.29</sup><sub>-0.71</sub></td><td>3.54,<br>4.31<sup>+1.54</sup><sub>-0.93</sub></td><td>3.03,<br>4.01<sup>+2.20</sup><sub>-1.10</sub></td><td>β<sub>t</sub></td><td>5.00, 4.75<sup>+2.42</sup><sub>-1.80</sub></td><td>4.65, 4.48<sup>+2.26</sup><sub>-1.63</sub></td><td>4.61, 4.80<sup>+2.72</sup><sub>-2.02</sub></td><td>4.55, 4.79<sup>+2.61</sup><sub>-1.88</sub></td></tr><tr><td>A<sub>k2h</sub></td><td>1.31,<br>1.47<sup>+0.40</sup><sub>-0.33</sub></td><td>1.28,<br>1.47<sup>+0.41</sup><sub>-0.34</sub></td><td>1.20,<br>1.36<sup>+0.38</sup><sub>-0.31</sub></td><td>1.30,<br>1.50<sup>+0.41</sup><sub>-0.34</sub></td><td>A<sub>t2h</sub></td><td>0.49, 0.84<sup>+0.71</sup><sub>-0.41</sub></td><td>0.47, 0.80<sup>+0.64</sup><sub>-0.39</sub></td><td>0.44, 0.86<sup>+0.86</sup><sub>-0.46</sub></td><td>0.47, 0.91<sup>+0.82</sup><sub>-0.46</sub></td></tr></tbody></table><table><thead><tr><th></th><th colspan="4">Density</th><th colspan="4">Model GNFW1h-matched</th></tr><tr><th></th><th>mstar-tot</th><th>mstar-red</th><th>mh-tot</th><th>mh-red</th><th></th><th>mstar-tot</th><th>mstar-red</th><th>mh-tot</th><th>mh-red</th></tr></thead><tbody><tr><td>log<sub>10</sub> ρ<sub>0</sub></td><td>6.00,<br>4.68<sup>+0.89</sup><sub>-1.01</sub></td><td>6.00,<br>4.53<sup>+0.95</sup><sub>-0.96</sub></td><td>6.00,<br>4.76<sup>+0.84</sup><sub>-1.00</sub></td><td>6.00,<br>4.45<sup>+1.00</sup><sub>-1.00</sub></td><td>P<sub>0</sub></td><td>18.33,<br>10.02<sup>+10.94</sup><sub>-5.93</sub></td><td>17.37,<br>10.24<sup>+10.81</sup><sub>-6.01</sub></td><td>13.26,<br>6.51<sup>+10.91</sup><sub>-4.07</sub></td><td>13.80,<br>7.32<sup>+11.28</sup><sub>-4.71</sub></td></tr><tr><td>α<sub>k</sub></td><td>0.17,<br>0.25<sup>+0.14</sup><sub>-0.06</sub></td><td>0.16,<br>0.25<sup>+0.14</sup><sub>-0.07</sub></td><td>0.17,<br>0.24<sup>+0.12</sup><sub>-0.05</sub></td><td>0.16,<br>0.25<sup>+0.15</sup><sub>-0.07</sub></td><td>x<sub>c,t</sub></td><td>0.06, 0.14<sup>+0.21</sup><sub>-0.08</sub></td><td>0.06, 0.14<sup>+0.19</sup><sub>-0.07</sub></td><td>0.05, 0.15<sup>+0.31</sup><sub>-0.10</sub></td><td>0.05, 0.14<sup>+0.26</sup><sub>-0.09</sub></td></tr><tr><td>β<sub>k</sub></td><td>2.06,<br>1.95<sup>+0.16</sup><sub>-0.18</sub></td><td>1.99,<br>1.90<sup>+0.15</sup><sub>-0.16</sub></td><td>2.03,<br>1.94<sup>+0.16</sup><sub>-0.17</sub></td><td>1.92,<br>1.83<sup>+0.16</sup><sub>-0.17</sub></td><td>β<sub>t</sub></td><td>1.53, 1.80<sup>+0.44</sup><sub>-0.33</sub></td><td>1.54, 1.81<sup>+0.43</sup><sub>-0.33</sub></td><td>1.42, 1.73<sup>+0.49</sup><sub>-0.35</sub></td><td>1.41, 1.72<sup>+0.47</sup><sub>-0.35</sub></td></tr></tbody></table><table><thead><tr><th></th><th colspan="4">Density</th><th colspan="4">Model GNFW1h-unmatched</th></tr><tr><th></th><th>mstar-tot</th><th>mstar-red</th><th>mh-tot</th><th>mh-red</th><th></th><th>mstar-tot</th><th>mstar-red</th><th>mh-tot</th><th>mh-red</th></tr></thead><tbody><tr><td>log<sub>10</sub> ρ<sub>0</sub></td><td>6.00,<br>4.80<sup>+0.82</sup><sub>-1.01</sub></td><td>6.00,<br>4.55<sup>+0.96</sup><sub>-1.01</sub></td><td>6.00,<br>4.92<sup>+0.75</sup><sub>-0.99</sub></td><td>6.00,<br>4.40<sup>+1.04</sup><sub>-1.04</sub></td><td>P<sub>0</sub></td><td>19.54, 9.58<sup>+11.05</sup><sub>-5.70</sub></td><td>19.23,<br>10.53<sup>+10.73</sup><sub>-5.99</sub></td><td>16.64,<br>7.79<sup>+11.06</sup><sub>-4.59</sub></td><td>16.50,<br>8.88<sup>+11.09</sup><sub>-5.27</sub></td></tr><tr><td>α<sub>k</sub></td><td>0.16,<br>0.23<sup>+0.11</sup><sub>-0.05</sub></td><td>0.16,<br>0.24<sup>+0.14</sup><sub>-0.06</sub></td><td>0.16,<br>0.22<sup>+0.10</sup><sub>-0.04</sub></td><td>0.15,<br>0.24<sup>+0.16</sup><sub>-0.07</sub></td><td>x<sub>c,t</sub></td><td>0.06, 0.16<sup>+0.23</sup><sub>-0.09</sub></td><td>0.08, 0.19<sup>+0.25</sup><sub>-0.10</sub></td><td>0.07, 0.22<sup>+0.34</sup><sub>-0.13</sub></td><td>0.08, 0.21<sup>+0.31</sup><sub>-0.13</sub></td></tr><tr><td>β<sub>k</sub></td><td>1.97,<br>1.87<sup>+0.16</sup><sub>-0.17</sub></td><td>1.88,<br>1.78<sup>+0.15</sup><sub>-0.16</sub></td><td>1.94,<br>1.84<sup>+0.15</sup><sub>-0.17</sub></td><td>1.84,<br>1.73<sup>+0.16</sup><sub>-0.17</sub></td><td>β<sub>t</sub></td><td>1.39, 1.64<sup>+0.38</sup><sub>-0.30</sub></td><td>1.42, 1.67<sup>+0.38</sup><sub>-0.31</sub></td><td>1.29, 1.57<sup>+0.38</sup><sub>-0.30</sub></td><td>1.39, 1.69<sup>+0.41</sup><sub>-0.33</sub></td></tr></tbody></table>

Note. The best values (minimum $\chi^2$) for each fit are listed first, followed by the marginalized parameter listed as the median $\pm 1\sigma$.

## 4.2. Modeled Observations

Here we show the results of projecting and fitting the 2D simulated profiles. Figure 6 shows the projection process using Mop-c-GT; the left columns show the fits of the 3D simulated profiles from TNG for different halo selections and fitting models, which are input into the projection code, and the right columns show these different populations projected into tSZ and kSZ temperature shifts described in Section 3.2. The different line styles correspond to the different kind of fit, as specified in the legend of the figure, and all of the fits are selected by minimum $\chi^2$. The error bars are the forecasted errors for the SO experiment.

The top panels of Figure 6 show fits for the density profiles. While not as prominent as the differences in pressure, the inclusion of mass-distribution matching does have an effect on the observed density. In general, a profile with higher 3D values will have higher measurable signals, which is what we see in this figure. The GNFW and GNFW1h-matched models have higher values than the unmatched fits in 3D, and thus have higher kSZ signals. Within each of the fitting models the stellar mass-selected sample has higher kSZ signals than the halo mass-selected sample, which is also seen in the 3D profiles. This is the result of SHMR of IllustrisTNG being lower than the Kravtsov et al. (2018) relation, so at a fixed stellar mass
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Mass dependence of marginalized parameters for each sample fit by the GNFW model (top) and GNFW1h model (bottom). The point in each bin shows the median of the MCMC chains with error bars ±1σ showing the distribution within the bin. The red and blue lines (showing different mass selections) have been slightly offset in each mass bin to more clearly show the different error bars, along with different cap sizes to more clearly show the error bars of the red vs. tot samples.</figcaption></figure>

TNG will select a higher halo mass. Furthermore, the GNFW1h models overpredict the values of the profiles at high radii (seen previously in Figure 4), so the amplitudes of the profiles only taking into account the one-halo term are higher than the models including a two-halo term.

The lower panels of Figure 6 shows fits for the pressure profiles. As shown previously in Figure 4, the inclusion of mass-distribution matching has a significant effect on the 3D pressure profiles due to its dependence on $M^{5/3}$. We can see this in more detail in both bottom panels, along with similar trends to the top panel density profiles with the matched models and the stellar mass-selected sample resulting in higher values and amplitudes. It is clear from Figure 6 that the modeling choices for the cross-correlation signal are important simply by comparing the differences in the observed profiles to the size of the forecasted errors.

### 4.3. Testing the Number of Free Parameters

Certain parameters of the GNFW profile are degenerate, as discussed in Section 2.4. We test how changing the number of

free parameters in the fits of the observed profiles affects the results by fixing another parameter for each type of profile, shown in Figure 7. For density we fix $\alpha_k$ and for pressure we fix $x_{c,t}$, as each of these parameters have degeneracies with $\log_{10}\rho_0$ and $\beta_t$, respectively, as seen by the contours in this figure.

Figure 7 shows an example of contours of fits for the GNFW-$m$ model of the ms-red sample using a different number of free parameters. We find that some profile information is lost in the projection process, resulting in the input 3D parameters not being recovered as well in the 2D fit with all parameters free. The solid lines show the values of the best 3D fit parameters, which serve as inputs for the 2D fits. The red contours show the fits with all of the GNFW parameters free (see Table 1), and the blue contours show fits holding one extra parameter constant. The small panels show the corresponding values calculated as a GNFW profile to show the differences in the 3D profiles. It can clearly be seen by the contours that fixing an extra parameter in the fits results in better alignment with the true values, as the peaks of the blue contours closely match the solid lines. Similarly, the corresponding 3D profile computed using the parameters from the blue contours more closely aligns with the 3D profile computed using the input parameters than the red.

This result indicates that with this experiment we cannot resolve the inner/intermediate parts of the profile for neither density or pressure. Perhaps with a new experiment, such as CMB-HD (Sehgal et al. 2019) we would be able to resolve the profile further into the interior at the redshift of the CMASS sample. Furthermore, if we were to choose a different sample to model at lower redshift, we expect that our ability to resolve the inner parts of these profile would improve. Hereafter, we show the fits fixing an additional parameter to be able to constrain whether differences are due to the number of free parameters or the modeling choices.

### 4.4. Biases from Modeling Choices

In Figure 8 we show contours for pressure and density of different samples, all for the GNFW-$m$ model. These fits use different input data ($\bar{d}$) for the 2D observable profiles corresponding to each color and mass selection to show how the shapes of the profiles differ among differently-selected samples. We note the differences in two-halo term amplitudes in this work and Amodeo et al. (2021) are due to the TNG simulations having a larger two-halo term amplitude compared to the observations presented in Schaan et al. (2021).

For pressure, there is a clear separation among the $P_0$ values favored by the halo- and stellar-mass-selected samples, with the ms samples returning higher $P_0$ values than the mh samples. The contours show more of a spread in the remaining parameters, with more significant differences in $\beta_t$ than $A_{2h}$. For density, the normalization parameter $\log_{10}\rho_0$ shows the most significant spread among samples, with the mh-tot sample favoring the highest values and ms-red sample favoring the lowest values. For $\beta_k$ we see a spread similar to $\beta_t$ for pressure, and for $A_{2h}$ the fits for all of the samples return nearly the same values. The differences among the contours for these samples show again that the modeling choices yield different results.

Figure 9 shows contours only for the ms-red sample, but with different fitting models and distribution-matching options. The solid lines show the corresponding value of the 3D best fit for each parameter, which is used as input for the 2D fits. This
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Projection process of Mop-c-GT. Input 3D theoretical profiles (in this case, from TNG) and output the tSZ and kSZ observable temperature shifts as a function of radius in arcminutes. The different colors correspond to different mass selections, and different line styles correspond to different fitting models and inclusion of distribution matching. The error bars correspond to the forecasted SO error bars at 150 GHz, discussed further in Section 3.3. These profiles show the effects of 1. selecting by different mass types, with the ms samples having higher values and amplitudes than the mh samples, 2. mass-distribution matching, with more significant differences for pressure but still relevant for density, and 3. inclusion of a two-halo term, which results in lower amplitudes of the observed profiles as the 1h-only model overpredicts the values.</figcaption></figure>

figure differs from Figure 8 in that the fits all use the same input profile $\bar{d}$ for the 2D observable profiles, calculated by the GNFW-$m$ model, but use different fitting models ($\bar{\mu}$), such as without a two-halo term or without mass-distribution matching. By treating the GNFW-$m$ profile as our truth, we can show how well the other models recover the given parameters. Simply by looking at the contours we can see that assuming an incorrect model results in estimates of the parameters not being very close to the actual values. For density, the GNFW-um model returns nearly correct values for $\log_{10}\rho_0$ and $A_{k2h}$, although slightly lower values for $\beta_k$, indicating that inclusion of distribution matching is not as important to kSZ model fitting as the inclusion of a two-halo term. The estimates from the fits with GNFW1h models are much lower than the true values for $\log_{10}\rho_0$ and $\beta_k$, with no information on a two-halo term.

For pressure we see more of a spread in returned values of different fitting models. The model that returns the closest value to the truth for $P_0$ is the GNFW1h-um model, for $\beta_t$ the

GNFW1h-$m$ model, and for $A_{k2h}$ the GNFW-um model returns nearly the correct value. This indicates that the inclusion of mass-distribution matching is more important for tSZ model fitting than kSZ, and is equally as important as the inclusion of a two-halo term (see the $\sigma$ values of Table 4). For both of the profiles, incorrectly assuming the model induces systematic effects in the resulting parameters.

The resulting values of the 2D fits with the fiducial ms-red, GNFW-$m$ model as the data are shown in Table 4 (contours shown in Figure 9), allowing for comparisons of the importance of different modeling choices. Next to the estimates for each parameter in the table is the number of $\sigma$ away from the corresponding fiducial value of the GNFW-$m$ model.

As a specific example, in the right panel of Figure 9 we show contours of the ms-red sample for density. Taking the GNFW-$m$ model as our truth, the input values for $\log_{10}\rho_0$, $\beta_k$, and $A_{k2h}$ were $3.28 \pm 0.04$, $3.20 \pm 0.11$, and $1.37 \pm 0.23$, respectively. If we fit this profile with the GNFW-um model, we would get
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Contours for fits of projections for tSZ (left) and kSZ (right) for the GNFW-m model of the ms-red samples. The contours for fits with all GNFW parameters free are shown in red, and the contours for fits fixing an additional parameter are shown in blue. The vertical lines show the input parameter, the best 3D fit from the corresponding sample. The small panels show the 3D profiles computed using the corresponding values for each of the different contours, using a GNFW model, and the same profiles normalized by the TNG profile to more clearly see the differences. This figure shows that the fits are not sensitive to information of the inner profiles, resulting in degeneracies, more significant for pressure. Fixing an additional parameter allows for better constraints on the remaining free parameters.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Comparison of parameter spaces for fits of the 2D projections of different samples, all for the GNFW-m model. An additional parameter is fixed, $x_{c,t}$ for pressure (left) and $\alpha_k$ for density (right), as discussed in Section 4.3. This figure shows the differences among the fits of different samples (ms vs. mh, and red vs. tot), highlighting the importance of these modeling choices.</figcaption></figure>

values for these parameters of 3.29 ± 0.04, 2.95 ± 0.10, and 1.33 ± 0.23, which would be differences of 0.2, 1.7, and 0.1σ. If we fit the profile with the GNFW1h-m model, we would get values for $\log_{10}\rho_0$ and $\beta_k$ of 3.12 ± 0.02 and 2.71 ± 0.06, which are differences of 3.6σ and 3.9σ, with no information on the two-halo term. Lastly, if we fit the profile with the GNFW1h-um model, we would get values for $\log_{10}\rho_0$ and $\beta_k$ of 3.14 ± 0.02 and 2.52 ± 0.05, differences of 3.1σ and 5.6σ, respectively, also with no information on the two-halo term. For this sample, it is clear that the inclusion of a two-halo term to the fitting model is the most important modeling systematic, as the differences between the input and selected values for the GNFW-um model are smaller than the differences of the GNFW1h models.

## 5. Conclusions

In summary, we have developed methods to extract halo information from the IllustrisTNG simulations. We derive halo samples of different properties of interest, such as within a specific mass range or selected by color, and stack the thermodynamic properties of density and pressure to create 3D radial profiles using *Illstack*. We fit these profiles to a GNFW profile, testing the importance of the halo sample modeling uncertainties of mass selection, color selection, mass-distribution matching, and addition of a two-halo term to the fitting model. We project the profiles into a 2D observing space using *Mop-c-GT*, resulting in profiles with units of CMB temperature shifts due to the kSZ and tSZ effects, and perform
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Comparison of parameter spaces for fits of the 2D projections ms-red sample with the GNFW-m model as the fiducial model. The solid lines indicate the values of the fits of the 3D profiles for the same sample (ms-red, GNFW-m). This figure shows that for density the most important modeling choice of the fitting model is the inclusion of a two-halo term, and for pressure both inclusion of distribution weighting and two-halo term are important.</figcaption></figure>

<table><caption>**Table 4**
Parameters from 2D Fits for Each Sample and Model</caption><thead><tr><th colspan="5">Model GNFW-matched<br>(Input Parameters for Fits Below)</th><th colspan="5">Model GNFW-unmatched</th></tr><tr><th></th><th colspan="4">Density</th><th colspan="4">Pressure</th></tr><tr><th></th><th>ms-tot</th><th>ms-red</th><th>mh-tot</th><th>mh-red</th><th></th><th>ms-tot</th><th>ms-red</th><th>mh-tot</th><th>mh-red</th></tr></thead><tbody><tr><td>log<sub>10</sub> ρ<sub>0</sub></td><td>3.67<sup>+0.05</sup><sub>-0.05</sub></td><td>3.28<sup>+0.04</sup><sub>-0.04</sub></td><td>4.34<sup>+0.10</sup><sub>-0.10</sub></td><td>3.39<sup>+0.07</sup><sub>-0.06</sub></td><td>P<sub>0</sub></td><td>4.20<sup>+0.25</sup><sub>-0.23</sub></td><td>4.05<sup>+0.24</sup><sub>-0.22</sub></td><td>2.78<sup>+0.28</sup><sub>-0.25</sub></td><td>2.89<sup>+0.25</sup><sub>-0.22</sub></td></tr><tr><td>β<sub>k</sub></td><td>3.37<sup>+0.11</sup><sub>-0.11</sub></td><td>3.20<sup>+0.11</sup><sub>-0.11</sub></td><td>3.45<sup>+0.16</sup><sub>-0.16</sub></td><td>3.30<sup>+0.16</sup><sub>-0.16</sub></td><td>β<sub>t</sub></td><td>6.10<sup>+0.16</sup><sub>-0.15</sub></td><td>6.16<sup>+0.16</sup><sub>-0.15</sub></td><td>6.36<sup>+0.27</sup><sub>-0.25</sub></td><td>5.83<sup>+0.21</sup><sub>-0.20</sub></td></tr><tr><td>A<sub>k2h</sub></td><td>1.41<sup>+0.23</sup><sub>-0.23</sub></td><td>1.37<sup>+0.23</sup><sub>-0.23</sub></td><td>1.32<sup>+0.23</sup><sub>-0.23</sub></td><td>1.35<sup>+0.23</sup><sub>-0.23</sub></td><td>A<sub>t2h</sub></td><td>0.57<sup>+0.06</sup><sub>-0.06</sub></td><td>0.53<sup>+0.06</sup><sub>-0.06</sub></td><td>0.52<sup>+0.06</sup><sub>-0.06</sub></td><td>0.51<sup>+0.06</sup><sub>-0.06</sub></td></tr><tr><th colspan="5">Model GNFW-unmatched</th><th colspan="4">Model GNFW-matched</th></tr><tr><th></th><th colspan="4">Density</th><th colspan="4">Pressure</th></tr><tr><th></th><th>ms-tot</th><th>ms-red</th><th>mh-tot</th><th>mh-red</th><th></th><th>ms-tot</th><th>ms-red</th><th>mh-tot</th><th>mh-red</th></tr><tr><td>log<sub>10</sub> ρ<sub>0</sub></td><td>3.69<sup>+0.05</sup><sub>-0.05</sub></td><td>3.29<sup>+0.04</sup><sub>-0.04</sub></td><td>4.42<sup>+0.10</sup><sub>-0.10</sub></td><td>3.41<sup>+0.07</sup><sub>-0.07</sub></td><td>P<sub>0</sub></td><td>6.80<sup>+0.40</sup><sub>-0.37</sub></td><td>5.9<sub>σ</sub></td><td>5.98<sup>+0.34</sup><sub>-0.32</sub></td><td>4.9<sub>σ</sub></td><td>7.52<sup>+0.76</sup><sub>-0.66</sub></td><td>6.7<sub>σ</sub></td><td>4.86<sup>+0.41</sup><sub>-0.37</sub></td><td>4.5<sub>σ</sub></td></tr><tr><td>β<sub>k</sub></td><td>3.12<sup>+0.10</sup><sub>-0.10</sub></td><td>2.95<sup>+0.10</sup><sub>-0.10</sub></td><td>3.13<sup>+0.14</sup><sub>-0.14</sub></td><td>3.03<sup>+0.15</sup><sub>-0.14</sub></td><td>β<sub>t</sub></td><td>4.88<sup>+0.10</sup><sub>-0.11</sub></td><td>6.0<sub>σ</sub></td><td>5.00<sup>+0.12</sup><sub>-0.12</sub></td><td>5.7<sub>σ</sub></td><td>4.60<sup>+0.18</sup><sub>-0.17</sub></td><td>5.3<sub>σ</sub></td><td>4.66<sup>+0.16</sup><sub>-0.15</sub></td><td>4.4<sub>σ</sub></td></tr><tr><td>A<sub>k2h</sub></td><td>1.38<sup>+0.23</sup><sub>-0.23</sub></td><td>0.1<sub>σ</sub></td><td>1.34<sup>+0.23</sup><sub>-0.23</sub></td><td>1.29<sup>+0.23</sup><sub>-0.24</sub></td><td>A<sub>t2h</sub></td><td>0.55<sup>+0.06</sup><sub>-0.06</sub></td><td>0.2<sub>σ</sub></td><td>0.51<sup>+0.06</sup><sub>-0.06</sub></td><td>0.2<sub>σ</sub></td><td>0.49<sup>+0.06</sup><sub>-0.06</sub></td><td>0.2<sub>σ</sub></td><td>0.49<sup>+0.06</sup><sub>-0.06</sub></td><td>0.2<sub>σ</sub></td></tr><tr><th colspan="5">Model GNFW-matched</th><th colspan="4">Model GNFW-unmatched</th></tr><tr><th></th><th colspan="4">Density</th><th colspan="4">Pressure</th></tr><tr><th></th><th>ms-tot</th><th>ms-red</th><th>mh-tot</th><th>mh-red</th><th></th><th>ms-tot</th><th>ms-red</th><th>mh-tot</th><th>mh-red</th></tr><tr><td>log<sub>10</sub> ρ<sub>0</sub></td><td>3.47<sup>+0.03</sup><sub>-0.03</sub></td><td>3.12<sup>+0.02</sup><sub>-0.02</sub></td><td>3.94<sup>+0.05</sup><sub>-0.05</sub></td><td>3.13<sup>+0.04</sup><sub>-0.04</sub></td><td>P<sub>0</sub></td><td>3.11<sup>+0.13</sup><sub>-0.13</sub></td><td>3.8<sub>σ</sub></td><td>3.08<sup>+0.13</sup><sub>-0.13</sub></td><td>3.5<sub>σ</sub></td><td>1.79<sup>+0.12</sup><sub>-0.11</sub></td><td>3.1<sub>σ</sub></td><td>1.98<sup>+0.12</sup><sub>-0.11</sub></td><td>3.2<sub>σ</sub></td></tr><tr><td>β<sub>k</sub></td><td>2.86<sup>+0.06</sup><sub>-0.06</sub></td><td>2.71<sup>+0.06</sup><sub>-0.06</sub></td><td>2.79<sup>+0.08</sup><sub>-0.08</sub></td><td>2.62<sup>+0.08</sup><sub>-0.08</sub></td><td>β<sub>t</sub></td><td>5.24<sup>+0.10</sup><sub>-0.10</sub></td><td>4.5<sub>σ</sub></td><td>5.35<sup>+0.10</sup><sub>-0.10</sub></td><td>4.2<sub>σ</sub></td><td>5.07<sup>+0.15</sup><sub>-0.15</sub></td><td>4.0<sub>σ</sub></td><td>4.79<sup>+0.12</sup><sub>-0.12</sub></td><td>4.2<sub>σ</sub></td></tr><tr><th colspan="5">Model GNFW-matched</th><th colspan="4">Model GNFW-unmatched</th></tr><tr><th></th><th colspan="4">Density</th><th colspan="4">Pressure</th></tr><tr><th></th><th>ms-tot</th><th>ms-red</th><th>mh-tot</th><th>mh-red</th><th></th><th>ms-tot</th><th>ms-red</th><th>mh-tot</th><th>mh-red</th></tr><tr><td>log<sub>10</sub> ρ<sub>0</sub></td><td>3.49<sup>+0.03</sup><sub>-0.03</sub></td><td>3.14<sup>+0.02</sup><sub>-0.02</sub></td><td>4.03<sup>+0.05</sup><sub>-0.05</sub></td><td>3.16<sup>+0.04</sup><sub>-0.04</sub></td><td>P<sub>0</sub></td><td>5.12<sup>+0.21</sup><sub>-0.20</sub></td><td>2.9<sub>σ</sub></td><td>4.62<sup>+0.19</sup><sub>-0.18</sub></td><td>2.0<sub>σ</sub></td><td>4.92<sup>+0.33</sup><sub>-0.30</sub></td><td>5.3<sub>σ</sub></td><td>3.39<sup>+0.19</sup><sub>-0.18</sub></td><td>1.7<sub>σ</sub></td></tr><tr><td>β<sub>k</sub></td><td>2.67<sup>+0.05</sup><sub>-0.05</sub></td><td>2.52<sup>+0.05</sup><sub>-0.05</sub></td><td>2.56<sup>+0.07</sup><sub>-0.07</sub></td><td>2.43<sup>+0.07</sup><sub>-0.07</sub></td><td>β<sub>t</sub></td><td>4.26<sup>+0.07</sup><sub>-0.07</sub></td><td>10.5<sub>σ</sub></td><td>4.41<sup>+0.07</sup><sub>-0.07</sub></td><td>9.9<sub>σ</sub></td><td>3.77<sup>+0.10</sup><sub>-0.10</sub></td><td>8.9<sub>σ</sub></td><td>3.90<sup>+0.09</sup><sub>-0.09</sub></td><td>8.4<sub>σ</sub></td></tr></tbody></table>

Note. We use the GNFW-m model as the fiducial model as the data, $\vec{d}$, in the likelihood. We list the marginalized parameter for each fit as the median $\pm 1\sigma$, followed by the number of $\sigma$ away from the true value (of the GNFW-m model). This table shows fits with an additional parameter fixed, $x_{c,t}$ for pressure and $\alpha_k$ for density, described further in Section 4.2.
the same fitting routine. We test the importance of the different modeling uncertainties in 2D by modeling how SO would detect an observed sample (subsample of CMASS, analyzed in Amodeo et al. 2021) of differently-derived halo populations. Below we reiterate the important takeaways of the study.

We show in several ways that mass selection, color selection, mass-distribution matching, and fitting to a model that includes the addition of a two-halo term having varying levels of importance when modeling a sample of halos. These results are seen by the differences in the theoretical 3D profiles (see Figure 4) and the observational 2D profiles (see Figure 6). The main trends of these figures are: 1. the stellar mass-selected samples tend to have higher density and pressure values than the halo mass-selected samples due to the shape of TNG's SHMR, and 2. weighting by the observed mass distribution has a larger effect on the pressure profiles due to its nonlinear dependence on mass, but is still relevant for the density profiles.

To further demonstrate the importance of matching the observed mass distribution, we show that some GNFW parameters fit to the 3D theoretical profiles show dependence on mass, shown in Figure 5. Therefore, modeling these mass trends by weighting to the observed mass distribution is important when modeling observed samples for cross correlations.

For the redshift and mass ranges of the observed CMASS sample, we show that the fits are not sensitive to the properties of the inner/intermediate profiles for both pressure and density. With the error bars that are forecasted for SO, the inclusion of such inner parameters, especially in pressure, results in strong degeneracies with other parameters, shown in Figure 7. After fixing an additional parameter in the GNFW fits, we show that we are able to better constrain the remaining free parameters and better understand the fits in relation to the modeling uncertainties of interest. We show again that different selections in mass and color return either clear separations between populations or a spread of estimates for each parameter.

Lastly, in Figure 9 we show in more detail that fitting to a model including weighting by the mass distribution and addition of a two-halo term are extremely important. Given a fiducial profile, not including either of these modeling choices returns incorrect values for each parameter. Shown by Table 4, for this sample the inclusion of a two-halo term tends to be more important than the inclusion of weighting by the mass distribution for the density profiles, seen by the lower σ values from the fiducial profile's parameters. For the pressure profiles the σ values for the models lacking matching by the mass distribution and a two-halo term are comparable, indicating that both of these systematics need to be addressed.

The uncertainties explored in this study are by no means an exhaustive list to consider when modeling a sample of halos. For example, as briefly discussed previously, measurements of masses on both sides of the SHMR are quite uncertain but in particular the measurements of stellar masses introduce significantly more systematic uncertainties. One could perform the same kind of study for this uncertainty by designing samples corresponding to different methods of estimating an observed sample's stellar mass, e.g., converting luminosity to stellar mass, and quantifying how the 3D and 2D profiles differ accordingly. We did explore the SHMR used in CMASS with respect to TNG by selecting and weighting halos according to

either their halo masses and stellar masses. In Figure 8 we can see the differences in the inferred profile parameters for both the density and pressure. Thus, going forward detailed modeling of the SHMR will be critical. Along similar lines the halo occupation distribution is another modeling concern for future SZ observations (Pandey et al. 2020).

In the future, we plan to expand these studies to explore the prescriptions of different parameters used in simulations other than TNG, and study any potential trends with redshift. To expand the variety of simulations, we plan to use the Cosmology and Astrophysics with MachinE Learning Simulations (CAMELS) suite of simulations (Villaescusa-Navarro et al. 2021). CAMELS is a suite of thousands of simulations varying astrophysical and cosmological parameters, and will be very useful in further studying modeling choices such as the ones discussed in this paper.

With a greater understanding of how to correctly model an observed sample, including the choices discussed in this study and more, we will be more prepared to interpret observations with ever-increasing signal-to-noise and resolution. Better interpretations of observations will allow for a better understanding of the current unsolved theoretical questions of the CGM and galaxy evolution as a whole.

E.M. and N.B. acknowledge support from NSF grant AST-1910021. N.B. acknowledges support from NASA ATP grant 80NSSC18K0695 and the Research and Technology Development fund at the Jet Propulsion Laboratory through the project entitled "Mapping the Baryonic Majority." E.S. is supported by the Chamberlain fellowship Q4 at Lawrence Berkeley National Laboratory. S.F. is supported by the Physics Division of Lawrence Berkeley National Laboratory.

## ORCID iDs

Emily Moser https://orcid.org/0000-0003-1593-1505
Stefania Amodeo https://orcid.org/0000-0002-4200-9965
Nicholas Battaglia https://orcid.org/0000-0001-5846-0411
Marcelo A. Alvarez https://orcid.org/0000-0002-2796-9650
Simone Ferraro https://orcid.org/0000-0003-4992-7854
Emmanuel Schaan https://orcid.org/0000-0002-4619-8927

## References

Abazajian, K., Addison, G., Adshead, P., et al. 2019, arXiv:1907.04473
Ade, P., Aguirre, J., Ahmed, Z., et al. 2019, JCAP, 2019, 056
Agertz, O., & Kravtsov, A. V. 2016, ApJ, 824, 79
Ahn, C. P., Alexandroff, R., Allende Prieto, C., et al. 2014, ApJS, 211, 17
Amodeo, S., Battaglia, N., Schaan, E., et al. 2021, PhRvD, 103, 063514
Baldry, I. K., Glazebrook, K., Brinkmann, J., et al. 2004, ApJ, 600, 681
Battaglia, N. 2016, JCAP, 2016, 058
Battaglia, N., Bond, J. R., Pfrommer, C., & Sievers, J. L. 2012a, ApJ, 758, 75
Battaglia, N., Bond, J. R., Pfrommer, C., & Sievers, J. L. 2012b, ApJ, 758, 74
Battaglia, N., Ferraro, S., Schaan, E., & Spergel, D. N. 2017, JCAP, 2017, 040
Battaglia, N., Hill, J. C., Amodeo, S., et al. 2019, BAAS, 51, 297
Behroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, ApJ, 717, 379
Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57
Bregman, J. N. 2007, ARA&A, 45, 221
Cen, R., & Ostriker, J. P. 2006, ApJ, 650, 560
Chen, H.-W., Lanzetta, K. M., Webb, J. K., & Barcons, X. 1998, ApJ, 498, 77
Cooray, A., & Sheth, R. 2002, PhR, 372, 1
Davé, R. 2009, in ASP Conf. Ser. 419, Galaxy Evolution: Emerging Insights and Future Challenges, ed. S. Jogee et al. (San Francisco, CA: ASP), 347
Di Matteo, T., Springel, V., & Hernquist, L. 2005, Natur, 433, 604
Federrath, C. 2015, MNRAS, 450, 4035
Ford, A. B., Oppenheimer, B. D., Davé, R., et al. 2013, MNRAS, 432, 89
Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Fukugita, M., Hogan, C. J., & Peebles, P. J. E. 1998, ApJ, 503, 518

Fukugita, M., & Peebles, P. J. E. 2004, ApJ, 616, 643

Geha, M., Blanton, M. R., Yan, R., & Tinker, J. L. 2012, ApJ, 757, 85

Gelman, A., & Rubin, D. B. 1992, StaSc, 7, 457

Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS, 445, 175

Hand, N., Addison, G. E., Aubourg, E., et al. 2012, PhRvL, 109, 041101

Hernquist, L. 1990, ApJ, 356, 359

Hill, J. C., Baxter, E. J., Lidz, A., Greco, J. P., & Jain, B. 2018, PhRvD, 97, 083501

Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19

Hummels, C. B., Smith, B. D., Hopkins, P. F., et al. 2019, ApJ, 882, 156

Hummels, C. B., Smith, B. D., & Silvia, D. W. 2017, ApJ, 847, 59

Kravtsov, A. V., Vikhlinin, A. A., & Meshcheryakov, A. V. 2018, AstL, 44, 8

Lanzetta, K. M., Bowen, D. V., Tytler, D., & Webb, J. K. 1995, ApJ, 442, 538

Leauthaud, A., Tinker, J., Behroozi, P. S., Busha, M. T., & Wechsler, R. H. 2011, ApJ, 738, 45

Leauthaud, A., Tinker, J., Bundy, K., et al. 2012, ApJ, 744, 159

Marasco, A., Debattista, V. P., Fraternali, F., et al. 2015, MNRAS, 451, 4223

Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113

Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1

Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Nelson, D., Pillepich, A., Genel, S., et al. 2015, A&C, 13, 12

Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624

Nelson, D., Springel, V., Pillepich, A., et al. 2019, ComAC, 6, 2

Nielsen, N. M., Churchill, C. W., Kacprzak, G. G., Murphy, M. T., & Evans, J. L. 2015, ApJ, 812, 83

Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJL, 660, L43

Oppenheimer, B. D., & Davé, R. 2008, MNRAS, 387, 577

Oppenheimer, B. D., Schaye, J., Crain, R. A., Werk, J. K., & Richings, A. J. 2018, MNRAS, 481, 835

Padmanabhan, N., Schlegel, D. J., Seljak, U., et al. 2007, MNRAS, 378, 852

Pandey, S., Baxter, E. J., & Hill, J. C. 2020, PhRvD, 101, 043525

Peeples, M. S., Corlies, L., Tumlinson, J., et al. 2019, ApJ, 873, 129

Peng, Y.-j., Lilly, S. J., Kovač, K., et al. 2010, ApJ, 721, 193

Persic, M., & Salucci, P. 1992, MNRAS, 258, 14P

Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016a, A&A, 586, A140

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016b, A&A, 594, A13

Rudie, G. C., Steidel, C. C., Trainor, R. F., et al. 2012, ApJ, 750, 67

Scannapieco, C., Tissera, P. B., White, S. D. M., & Springel, V. 2008, MNRAS, 389, 1137

Schaan, E., Ferraro, S., Amodeo, S., et al. 2021, PhRvD, 103, 063513

Schaan, E., Ferraro, S., Vargas-Magaña, M., et al. 2016, PhRvD, 93, 082002

Sehgal, N., Aiola, S., Akrami, Y., et al. 2019, BAAS, Vol. 51, 6

Shankar, F., Guo, H., Bouillot, V., et al. 2014, ApJL, 797, L27

Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51

Somerville, R. S., Hopkins, P. F., Cox, T. J., Robertson, B. E., & Hernquist, L. 2008, MNRAS, 391, 481

Springel, V., Di Matteo, T., & Hernquist, L. 2005, ApJL, 620, L79

Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676

Steidel, C. C., Erb, D. K., Shapley, A. E., et al. 2010, ApJ, 717, 289

Sunyaev, R. A., & Zeldovich, Y. B. 1970, Ap&SS, 7, 3

Sunyaev, R. A., & Zeldovich, Y. B. 1972, CoASP, 4, 173

Sunyaev, R. A., & Zeldovich, Y. B. 1980, MNRAS, 190, 413

Suresh, J., Rubin, K. H. R., Kannan, R., et al. 2017, MNRAS, 465, 2966

Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017, ARA&A, 55, 389

Tumlinson, J., Thom, C., Werk, J. K., et al. 2013, ApJ, 777, 59

Tumlinson, J., Werk, J. K., Thom, C., et al. 2011, ApJ, 733, 111

Vikram, V., Lidz, A., & Jain, B. 2017, MNRAS, 467, 2315

Villaescusa-Navarro, F., Anglés-Alcázar, D., Genel, S., et al. 2021, ApJ, 915, 71

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518

Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291

Werk, J. K., Prochaska, J. X., Thom, C., et al. 2013, ApJS, 204, 17

Werk, J. K., Prochaska, J. X., Tumlinson, J., et al. 2014, ApJ, 792, 8

Zhao, H. 1996, MNRAS, 278, 488