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

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 holes (SMBHs), galactic winds, and magnetic fields based on ideal magnetohydrodynamics (Weinberger et al. 2017; Marinacci et al. 2018; Pillepich et al. 2018). ⁵

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³ Mpc³ and tracks 2 × 455³ 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³ Mpc³ and 2 × 455³ tracked particles. The TNG simulation with the next lowest level of resolution for the same box size is TNG100-2, and has 2 × 910³ 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.

We have identified simulated halos with certain properties of interest and created thermodynamic radial profiles using the publicly available stacking code repository Illstack. ⁶ An already existing public TNG repository illustris-python ⁷ 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 × 10⁻⁴−1 × 10² Mpc split into 25 bins in logspace, resulting in radial bin size ∼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\leqslant {\mathrm{log}}_{10}\left(\tfrac{M}{{M}_{\odot }}\right)\leqslant 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 ±1σ 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.

Zoom In Zoom Out Reset image size

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.

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

2.3.2. Color-splitting Definitions

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 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 × 10¹³ M_⊙ 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 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).

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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., ∑_i w_i Δm = 1, where w_i are the weights of each mass bin and Δ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 ${\mathrm{val}}_{m}=\tfrac{\sum {val}* w}{\sum w}$. In contrast, the unmatched (unweighted) average would be ${\mathrm{val}}_{\mathrm{um}}=\tfrac{\mathrm{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\leqslant {\mathrm{log}}_{10}\left(\tfrac{M}{{M}_{\odot }}\right)\leqslant 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).

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

$$\begin{eqnarray}&&\mathrm{ln}L(\vec{d})=-\displaystyle \frac{1}{2}{\left[\vec{d}-\vec{\mu }(\vec{\theta })\right]}^{{\rm{T}}}{\vec{{\rm{\Sigma }}}}^{-1}\left[\vec{d}-\vec{\mu }(\vec{\theta })\right],\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&P(\vec{\theta })\propto L(\vec{d})\Pr (\vec{\theta }),\end{eqnarray} \tag{ 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{eqnarray}\begin{array}{rcl}\displaystyle \frac{\rho (x)}{{\rho }_{{cr}}(z)} & = & {f}_{b}{\rho }_{0}{\left(x/{x}_{c,k}\right)}^{{\gamma }_{k}}{\left[1+{\left(x/{x}_{c,k}\right)}^{{\alpha }_{k}}\right]}^{-\tfrac{{\beta }_{k}-{\gamma }_{k}}{{\alpha }_{k}}},\\ {\rho }_{{cr}}(z) & = & \displaystyle \frac{3{H}_{0}^{2}}{8\pi G}\left[{{\rm{\Omega }}}_{m}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}\right],\end{array}\end{eqnarray} \tag{ 5 }$$

where scaled radius x ≡ r/r_200c, f_b is the baryon fraction Ω_b /Ω_m , ρ_cr(z) is the critical density of the universe at redshift z, H₀ is the Hubble constant, and G is the gravitational constant.

The GNFW thermal pressure profile has the functional form

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{P(x)}{{P}_{200c}} & = & {P}_{0}{\left(x/{x}_{c,t}\right)}^{{\gamma }_{t}}{\left[1+{\left(x/{x}_{c,t}\right)}^{{\alpha }_{t}}\right]}^{-{\beta }_{t}},\\ {P}_{200c} & = & \displaystyle \frac{200{{GM}}_{200{\rm{c}}}{\rho }_{\mathrm{cr}}(z){f}_{b}}{2{r}_{200{\rm{c}}}},\end{array}\end{eqnarray} \tag{ 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, α is the slope of the GNFW profile at x ∼ 1, β is the outer slope of the profile at x ≫ 1, and γ is the inner slope of the profile at x ≪ 1 (Zhao 1996).

For the gas density profiles (Equation (5)) we perform fits for parameters ρ₀, α_k , and β_k , keeping parameters x_c,k and γ_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 we use the fitting formulas derived in Battaglia (2016) as the starting points for each of the MCMC chains.

Table 1. Values for the Fixed Parameters and Free Parameters in the GNFW Profile Fits

Pressure	Density
Fixed
αt = 1.0	xc,k = 0.5
γt = −0.3	γk = −0.2
Free
P0, xc,t , βt , At2h	ρ0, αk , βk , Ak2h
P0: (0.1,30)	${\mathrm{log}}_{10}{\rho }_{0}:(1.0,6.0)$
xc, t : (0.01,1.0)	αk : (0.1,6.0)
βt : (1.0,10.0)	βk : (1.0,10.0)
At2h : (0.01,5.0)	Ak2h : (0.01,5.0)

Note. A_2h is the amplitude of the two-halo term. Priors for the fits are given for the free parameters.

Download table as:  ASCIITypeset image

For the thermal pressure profiles (Equation (6)) we perform fits for parameters P₀, x_c,t , and β_t , keeping parameters α_t and γ_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{eqnarray}\begin{array}{rcl}{P}_{\mathrm{GNFW}1{\rm{h}}} & = & {P}_{1h},\\ {P}_{\mathrm{GNFW}} & = & {P}_{1h}+{A}_{t2h}{P}_{2h},\end{array}\end{eqnarray} \tag{ 7 }$$

where P_GNFW1h is the model used when just taking the one-halo term into account, and P_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 ≲1.1. We combine the independent chains and derive marginalized estimates for each parameter using the median ±1σ as the error estimates. We select the model resulting in the minimum χ² as our best fit, shown as the fits in Figure 4.

Zoom In Zoom Out Reset image size

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 ≤ M ≤ 11.72, and the halo mass limits of the sample are 12.12 ≤ M ≤ 13.98, where $M={\mathrm{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.

Table 3. Parameters from 3D Fits for Each Sample and Model

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

Note. The best values (minimum χ²) for each fit are listed first, followed by the marginalized parameter listed as the median ±1σ.

Download table as:  ASCIITypeset image

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, ⁸ (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 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 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).

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{{\rm{\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.

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

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 ρ₀ and β_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σ 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.

Zoom In Zoom Out Reset image size

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 β_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.

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 χ². The error bars are the forecasted errors for the SO experiment.

Zoom In Zoom Out Reset image size

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 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.

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 α_k and for pressure we fix x_c,t , as each of these parameters have degeneracies with ${\mathrm{log}}_{10}{\rho }_{0}$ and β_t , respectively, as seen by the contours in this figure.

Zoom In Zoom Out Reset image size

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.

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 ($\vec{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).

Zoom In Zoom Out Reset image size

For pressure, there is a clear separation among the P₀ values favored by the halo- and stellar-mass-selected samples, with the ms samples returning higher P₀ values than the mh samples. The contours show more of a spread in the remaining parameters, with more significant differences in β_t than A_t2h. For density, the normalization parameter ${\mathrm{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 β_k we see a spread similar to β_t for pressure, and for A_k2h 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 differs from Figure 8 in that the fits all use the same input profile $\vec{d}$ for the 2D observable profiles, calculated by the GNFW-m model, but use different fitting models ($\vec{\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 ${\mathrm{log}}_{10}{\rho }_{0}$ and A_k2h, although slightly lower values for β_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 ${\mathrm{log}}_{10}{\rho }_{0}$ and β_k , with no information on a two-halo term.

Zoom In Zoom Out Reset image size

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₀ is the GNFW1h-um model, for β_t the GNFW1h-m model, and for A_t2h 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 σ values of Table 4). For both of the profiles, incorrectly assuming the model induces systematic effects in the resulting parameters.

Table 4. Parameters from 2D Fits for Each Sample and Model

Model GNFW-matched
(Input Parameters for Fits Below)
Density					Pressure
	ms-tot	ms-red	mh-tot	mh-red		ms-tot	ms-red	mh-tot	mh-red
${\mathrm{log}}_{10}{\rho }_{0}$	${3.67}_{-0.05}^{+0.05}$	${3.28}_{-0.04}^{+0.04}$	${4.34}_{-0.10}^{+0.10}$	${3.39}_{-0.06}^{+0.07}$	P0	${4.20}_{-0.23}^{+0.25}$	${4.05}_{-0.22}^{+0.24}$	${2.78}_{-0.25}^{+0.28}$	${2.89}_{-0.22}^{+0.25}$
βk	${3.37}_{-0.11}^{+0.11}$	${3.20}_{-0.11}^{+0.11}$	${3.45}_{-0.16}^{+0.16}$	${3.30}_{-0.16}^{+0.16}$	βt	${6.10}_{-0.15}^{+0.16}$	${6.16}_{-0.15}^{+0.16}$	${6.36}_{-0.25}^{+0.27}$	${5.83}_{-0.20}^{+0.21}$
Ak2h	${1.41}_{-0.23}^{+0.23}$	${1.37}_{-0.23}^{+0.23}$	${1.32}_{-0.23}^{+0.23}$	${1.35}_{-0.23}^{+0.23}$	At2h	${0.57}_{-0.06}^{+0.06}$	${0.53}_{-0.06}^{+0.06}$	${0.52}_{-0.06}^{+0.06}$	${0.51}_{-0.06}^{+0.06}$
Model GNFW-unmatched
Density					Pressure
	ms-tot	ms-red	mh-tot	mh-red		ms-tot	ms-red	mh-tot	mh-red
${\mathrm{log}}_{10}{\rho }_{0}$	${3.69}_{-0.05}^{+0.05}$, 0.3σ	${3.29}_{-0.04}^{+0.04}$, 0.2σ	${4.42}_{-0.10}^{+0.10}$, 0.6σ	${3.41}_{-0.07}^{+0.07}$, 0.3σ	P0	${6.80}_{-0.37}^{+0.40}$, 5.9σ	${5.98}_{-0.32}^{+0.34}$, 4.9σ	${7.52}_{-0.66}^{+0.76}$, 6.7σ	${4.86}_{-0.37}^{+0.41}$, 4.5σ
βk	${3.12}_{-0.10}^{+0.10}$, 1.7σ	${2.95}_{-0.10}^{+0.10}$, 1.7σ	${3.13}_{-0.14}^{+0.14}$, 1.5σ	${3.03}_{-0.14}^{+0.15}$, 1.2σ	βt	${4.88}_{-0.11}^{+0.12}$, 6.0σ	${5.00}_{-0.12}^{+0.12}$, 5.7σ	${4.60}_{-0.17}^{+0.18}$, 5.3σ	${4.66}_{-0.15}^{+0.16}$, 4.4σ
Ak2h	${1.38}_{-0.23}^{+0.23}$, 0.1σ	${1.34}_{-0.23}^{+0.23}$, 0.1σ	${1.29}_{-0.24}^{+0.23}$, 0.1σ	${1.33}_{-0.24}^{+0.23}$, 0.1σ	At2h	${0.55}_{-0.06}^{+0.06}$, 0.2σ	${0.51}_{-0.06}^{+0.06}$, 0.2σ	${0.49}_{-0.06}^{+0.06}$, 0.2σ	${0.49}_{-0.06}^{+0.06}$, 0.2σ
Model GNFW1h-matched
Density					Pressure
	ms-tot	ms-red	mh-tot	mh-red		ms-tot	ms-red	mh-tot	mh-red
${\mathrm{log}}_{10}{\rho }_{0}$	${3.47}_{-0.03}^{+0.03}$, 3.4σ	${3.12}_{-0.02}^{+0.02}$, 3.6σ	${3.94}_{-0.05}^{+0.05}$, 3.6σ	${3.13}_{-0.04}^{+0.04}$, 3.1σ	P0	${3.11}_{-0.13}^{+0.13}$, 3.8σ	${3.08}_{-0.13}^{+0.13}$, 3.5σ	${1.79}_{-0.11}^{+0.12}$, 3.1σ	${1.98}_{-0.11}^{+0.12}$, 3.2σ
βk	${2.86}_{-0.06}^{+0.06}$, 4.1σ	${2.71}_{-0.06}^{+0.06}$, 3.9σ	${2.79}_{-0.08}^{+0.08}$, 3.6σ	${2.62}_{-0.08}^{+0.08}$, 3.8σ	βt	${5.24}_{-0.10}^{+0.10}$, 4.5σ	${5.35}_{-0.10}^{+0.10}$, 4.2σ	${5.07}_{-0.15}^{+0.15}$, 4.0σ	${4.79}_{-0.12}^{+0.12}$, 4.2σ
Model GNFW1h-unmatched
Density					Pressure
	ms-tot	ms-red	mh-tot	mh-red		ms-tot	ms-red	mh-tot	mh-red
${\mathrm{log}}_{10}{\rho }_{0}$	${3.49}_{-0.03}^{+0.03}$, 3.1σ	${3.14}_{-0.02}^{+0.02}$, 3.1σ	${4.03}_{-0.05}^{+0.05}$, 2.8σ	${3.16}_{-0.04}^{+0.04}$, 2.7σ	P0	${5.12}_{-0.20}^{+0.21}$, 2.9σ	${4.62}_{-0.18}^{+0.19}$, 2.0σ	${4.92}_{-0.30}^{+0.33}$, 5.3σ	${3.39}_{-0.18}^{+0.19}$, 1.7σ
βk	${2.67}_{-0.05}^{+0.05}$, 5.8σ	${2.52}_{-0.05}^{+0.05}$, 5.6σ	${2.56}_{-0.07}^{+0.07}$, 5.0σ	${2.43}_{-0.07}^{+0.07}$, 5.0σ	βt	${4.26}_{-0.07}^{+0.07}$, 10.5σ	${4.41}_{-0.07}^{+0.07}$, 9.9σ	${3.77}_{-0.10}^{+0.10}$, 8.9σ	${3.90}_{-0.09}^{+0.09}$, 8.4σ

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 ±1σ, followed by the number of σ 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 α_k for density, described further in Section 4.2.

Download table as:  ASCIITypeset image

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 σ 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 ${\mathrm{log}}_{10}{\rho }_{0}$, β_k , and A_k2h were 3.28 ± 0.04, 3.20 ± 0.11, and 1.37 ± 0.23, respectively. If we fit this profile with the GNFW-um model, we would get 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 ${\mathrm{log}}_{10}{\rho }_{0}$ and β_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 ${\mathrm{log}}_{10}{\rho }_{0}$ and β_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.