# Probing galaxy cluster and intra-cluster gas with luminous red galaxies

Yan Gong, Yin-Zhe Ma, Hideki Tanimura

► To cite this version:

Yan Gong, Yin-Zhe Ma, Hideki Tanimura. Probing galaxy cluster and intra-cluster gas with luminous red galaxies. Monthly Notices of the Royal Astronomical Society, 2019, 486 (4), pp.4904-4916. 10.1093/mnras/stz1177. hal-02144076

HAL Id: hal-02144076

https://hal.science/hal-02144076v1

Submitted on 6 Sep 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Probing galaxy cluster and intra-cluster gas with luminous red galaxies

Yan Gong<sup>*</sup>, 1,2★ Yin-Zhe Ma<sup>2,3★</sup> and Hideki Tanimura<sup>4★</sup>

Accepted 2019 April 26. Received 2019 April 26; in original form 2018 October 2

## ABSTRACT

We use the cross-correlation between the thermal Sunyaev–Zeldovich (tSZ) signal measured by the Planck satellite and the luminous red galaxy (LRG) samples provided by the SDSS DR7 to study the properties of galaxy cluster and intra-cluster gas. We separate the samples into three redshift bins $z_1 = (0.16, 0.26)$, $z_2 = (0.26, 0.36)$, $z_3 = (0.36, 0.47)$, and stack the Planck y-map against LRGs to derive the averaged y-profile for each redshift bin. We then fit the stacked profile with the theoretical prediction from the universal pressure profile (UPP) by using the Markov Chain Monte Carlo method. We find that the best-fitting values of the UPP parameters for the three bins are generally consistent with the previous studies, except for the noticeable evolution of the parameters in the three redshift bins. We simultaneously fit the data in the three redshift bins together, and find that the original UPP model cannot fit the data at small angular scales very well in the first and third redshift bins. The joint fits can be improved by including an additional parameter $\eta$ to change the redshift dependence of the model (i.e. $E(z)^{8/3} \rightarrow E(z)^{8/3+\eta}$) with best-fitting value as $\eta = -3.11^{+1.09}_{-1.13}$. This suggests that the original UPP model with less redshift dependence may provide a better fit to the stacked tSZ profile.

**Key words:** galaxies: clusters: general – cosmic background radiation – large-scale structure of Universe.

## 1 INTRODUCTION

The clusters of galaxies are the important objects in learning the structure formation and cosmology in general. Around hot clusters, the gas is generally ionized with temperature above $10^6$ K (McCarthy et al. 2007) and provides high density of hot electrons. The intra-cluster medium (ICM) can be studied by using direct X-ray imaging or the Sunyaev–Zeldovich effect (SZ effect; Sunyaev & Zeldovich 1972, 1980). Intensity of the former is proportional to the square of the electron density ~$n_e^2$ which could be sensitive to the centre baryons. In contrast, observation of the SZ effect depends on the integration of pressure profile, which could be used to trace down the gas distribution in the lower density region. These regions, such as filament, sheets and voids, are the locations believed to host most of the baryons (Haider et al. 2016; Tanimura et al. 2019a). For this reason, there has been a growing interest in predicting and measuring the SZ effect by using various techniques in microwave

and radio bands (Birkinshaw & Gull 1978; Birkinshaw 1999; Carlstrom, Holder & Reese 2002; Ma et al. 2015).

The thermal SZ (tSZ) effect is a secondary anisotropy in the cosmic microwave background (CMB) radiation, which is caused by the hot free electrons around ICM scattering off the CMB photons, boosting the low-frequency CMB into higher frequency, i.e.

$$ \frac{\Delta T}{T_{\text{CMB}}} = \left[ \eta \frac{e^{\eta} + 1}{e^{\eta} - 1} - 4 \right] y \equiv g_{\nu} y, \quad (1) $$

where $g_{\nu} \equiv [\eta(e^{\eta} + 1)/(e^{\eta} - 1)] - 4$ captures the frequency dependence, and

$$ \eta = \frac{h\nu}{k_B T_{\text{CMB}}} = 1.76 \left( \frac{\nu}{100 \text{ GHz}} \right). \quad (2) $$

The dimensionless Comptonization parameter $y$ is

$$ y = \int n_e(\mathbf{r}) \sigma_T \frac{k_B T_e(\mathbf{r})}{m_e c^2} dl, \quad (3) $$

where $\sigma_T$ is the Thomson cross-section, $k_B$ is the Boltzmann constant, $m_e$ is the electron rest mass and the integral is taken along the line of sight. The sign of $g_{\nu}$ determines whether $\Delta T$ is positive or negative (an increment or decrement) for the CMB temperature. With $T_{\text{CMB}} = 2.725$ K, we have $g_{\nu} \ge 0$ for $\nu \ge 217$ GHz,

<footnote n="1"><sup>1</sup>Key Laboratory for Computational Astrophysics, National Astronomical Observatories, Chinese Academy of Sciences, 20A Datun Road, Beijing 100012, China</footnote>

<footnote n="2"><sup>2</sup>NAOC-UKZN Computational Astrophysics Centre (NUCAC), University of KwaZulu-Natal, Durban 4000, South Africa</footnote>

<footnote n="3"><sup>3</sup>School of Chemistry and Physics, University of KwaZulu-Natal, Westville Campus, Private Bag X54001, Durban 4000, South Africa</footnote>

<footnote n="4"><sup>4</sup>Institut d'Astrophysique Spatiale, CNRS (UMR 8617), Université Paris-Sud, Bâtiment 121, F-91405 Orsay, France</footnote>

<footnote>* E-mail: gongyan@bao.ac.cn (YG); ma@ukzn.ac.za (Y-ZM); hideki.tanimura@ias.u-psud.fr (HT)</footnote>
and vice versa. We can define the electron pressure profile as $P_e(\mathbf{r}) = n_e(\mathbf{r}) k_B T_e(\mathbf{r})$ in equation (3), and thus the pressure profile $P_e$ is an essential factor in determining the Compton y-parameter. Therefore, modelling and testing the pressure profile is very crucial to obtain an accurate profile of observed y-parameter.

One the other hand, by accurately measuring the Compton y-parameter, a more precise modelling of the pressure profile can be obtained. In 2010, Arnaud et al. (2010) investigated the cluster pressure profile with 33 local ($z < 0.2$) samples of clusters drawn from REXCESS catalogue observed by XMM-Newton. These samples span the mass range of $10^{14}M_{\odot} < M_{500} < 10^{15}M_{\odot}$ where the $M_{500}$ is the mass of the centre of cluster within the radius of density equal to $500\rho_{crit}$. Arnaud et al. (2010) proposed the universal pressure profile (UPP) for the clusters and determined the parameter values by using the 33 massive clusters. Following Arnaud et al. (2010), Planck Collaboration V (2013) used the Planck's 14-month nominal survey and 62 massive nearby clusters, and re-investigated the cluster UPP, and refined the parameter values. Besides these observational probes, Le Brun, McCarthy & Melin (2015) used the cosmo-OWLS simulation and tested the UPP against different AGN models, and refined the UPP model parameters.

In these studies, the tSZ effect is a very useful tool to explore the high-redshift galaxy clusters. Because the signal y is expected to be nearly independent of the redshift, the effect does not diminish with the increasing redshift. Therefore, the effect is very suitable for finding high-redshift clusters (Planck Collaboration VIII 2011) and cross-correlation with other large-scale structure tracers [e.g. cross-correlation with weak-lensing measurement (Ma et al. 2015; Hojjati et al. 2017)]. In addition, luminous red galaxies (LRGs) are early-type massive galaxies consisting mainly of old stars with little ongoing star formation. Thus LRGs are good tracers of galaxy clusters and underlying dark-matter distribution, and should be correlated with thermal SZ signal (Hoshino et al. 2015).

In this work, we will study this cross-correlation and explore the UPP of clusters. Since the Sloan Digital Sky Survey (SDSS) released its data, there have been increasingly large LRG catalogues that span more orders of magnitude in mass and occupy different distances. The samples we use here are DR7 LRG samples, which have mass range to $3 \times 10^{12}$ to $3 \times 10^{14} M_{\odot}$ and redshift range of $0.16 < z < 0.47$, and significantly overlap with the Planck all-sky Compton y-map. Therefore, providing the stacking results of the LRG samples at cosmological distances and fitting the UPP model is the main aim of this paper. We will provide robust tests of the UPP model parameters by this large LRG sample from SDSS DR7.

This paper is organized as follows. In Section 2, we present the Planck y-maps from the CMB and LRG catalogues from SDSS-DR7 used in this study, and the stacking results of LRG y-profile. In Section 3, we present the modelling of the stacked profiles by calculating the one-halo and two-halo terms of the stacked profile. In Section 4, we discuss the data analysis procedure and $\chi^2$ study in our work. In Section 5, we present the results of our numerical fitting and compare our results with other works. The summary and discussion are presented in the last section.

Throughout this paper, we adopt a spatially flat fiducial $\Lambda$CDM cosmology with Planck cosmological parameter values (Planck Collaboration XVI 2014), which are the fractional baryon density $\Omega_b = 0.0490$, fractional matter density $\Omega_m = 0.3175$, spectral index $n_s = 0.9624$, rms matter fluctuation amplitude $\sigma_8 = 0.834$, and $h = 0.6711$ defined as $H_0 = 100 h \text{ km s}^{-1} \text{ Mpc}^{-1}$, where $H_0$ is the Hubble constant. We note that our constraint results are not

very sensitive to the current determined values of cosmological parameters.

## 2 DATA

We adopt the same data sets used in Tanimura et al. (2019b) in our analysis: the LRG catalogue from the Sloan Digital Sky Survey seventh data release Kazin et al. (2010) and the Planck Compton y-map (Planck Collaboration XXII 2016) from the 2015 data release (Planck Collaboration I 2016). Each data is described briefly in this section.

### 2.1 LRG catalogue

LRGs are early-type, massive galaxies with little ongoing star formation, which are selected based on magnitude and colour cut. They are typically located in the centres of galaxy groups and clusters and have been used to detect and characterize the remnants of baryon acoustic oscillations at low to intermediate redshift (Eisenstein et al. 2005; Kazin et al. 2010; Anderson et al. 2014). The SDSS DR7 LRGs are selected in Kazin et al. (2010) and cover $\sim 20$ per cent of the sky with almost a flat distribution. The LRG catalogue provides 105 831 LRGs with galaxy positions, magnitudes, and spectroscopic redshifts. Stellar masses of the LRGs are estimated in the New York University Value-Added catalogue (NYU-VAGC)<fnref n="1" /> using the K-correct software<fnref n="2" /> of Blanton & Roweis (2007) based on a stellar initial mass function of Chabrier (2003) and stellar evolution synthesis models of Bruzual & Charlot (2005). They are given in unit of $h^{-2}M_{\odot}$ and we adopt $h = 0.6711$ (Planck Collaboration XVI 2014).

Since not all LRGs are located in the centres of galaxy clusters, some LRGs may reside in them as ‘satellite’ galaxies. To minimize the satellite LRGs in our sample, we select locally (geometrically) most-massive LRGs (based on stellar mass) using a criterion that is analogous to that used in Planck Collaboration XI (2013): we remove a galaxy if a more massive galaxy resides within a tangential distance of 1.0 Mpc and within a radial velocity difference of $|c\Delta z| < 1000 \text{ km s}^{-1}$. After this selection, 101 407 locally most-massive LRGs are left that are likely to be ‘central’ LRGs. The detail of the sample is described in Tanimura et al. (2019b).

### 2.2 Planck y-map

The Planck tSZ map from the Planck 2015 data release is provided in HEALPix format (Gorski et al. 2005)<fnref n="3" /> with a pixel resolution of $N_{side} = 2048$. Two types of algorithms, MILCA (Hurier, Macias-Perez & Hildebrandt 2013) and NILC (Remazeilles, Aghanim & Douspis 2013), are applied for the Planck band maps to extract the tSZ signal. They are based on ILC (internal linear combination) techniques aiming to preserve an astrophysical component for which the electromagnetic spectrum is known (Bennett et al. 2003). Our analysis is based on the MILCA y-map, but we find consistent results with the NILC y-map.

The 2015 data release also provides sky masks suitable for analysing the y-maps. The mask covers a point-source mask and galactic masks that exclude 40 per cent, 50 per cent, 60 per cent, and 70 per cent of the sky. We combine the point source mask

<footnote n="1">¹http://sdss.physics.nyu.edu/vagc/</footnote>

<footnote n="2">²http://howdy.physics.nyu.edu/index.php/Kcorrect</footnote>

<footnote n="3">³http://healpix.sourceforge.net/</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** The halo mass distributions of the LRG samples in different redshift ranges. The blue, green, and red dashed curves are for $z_1$, $z_2$, and $z_3$ bins, respectively. The mass distribution of the whole redshift range is denoted in solid purple curve. The grey hatched region shows the selected sub-samples, which are used in our analysis, with the same mass distribution for the three z-bins.</figcaption></figure>

with the 40 per cent galactic mask, which excludes ~50 per cent of the sky. The mask is used in our stacking process such that for a given LRG, masked pixels in the y-map near that LRG are not accumulated in the stacked image. Since the mask may bias the y-profile, we accept 74 681 LRGs, for which ≥80 per cent of the region within a 40 arcmin circle around each LRG is available.

## 2.3 Stacking y-map centred on LRGs

We describe our procedure for stacking the Planck y-map at the positions of LRGs. We set a two-dimensional (2D) angular coordinate system of -40 arcmin < $\Delta l$ < 40 arcmin and -40 arcmin < $\Delta b$ < 40 arcmin, divided in 80 × 80 bins. For each LRG in the catalogue, we place it at the centre of the coordinate and subtract its local background signal of the mean tSZ signal in the annular region between 30 and 40 arcmin. This procedure nulls a large scale mode of fluctuation in the y-map, not due to the target object. We repeat this procedure for our sample of LRGs, stack them and finally divide the stack by the number of our sample.

In order to investigate a redshift evolution of tSZ signal (evolution of pressure profile), we divide our LRG sample into three redshift bins of sub-samples, $z_1$ bin (0.16 < $z$ < 0.26), $z_2$ bin (0.26 < $z$ < 0.36), and $z_3$ bin (0.36 < $z$ < 0.47). These redshift bins are selected to have an equal interval in redshift to probe the redshift evolution of pressure profiles around $z \sim 0.2, 0.3$, and 0.4. The number of LRGs in each redshift bin is 18 083 in $z_1$ bin, 29 586 in $z_2$ bin, and 27 012 LRGs in $z_3$ bin.

To study the redshift evolution of pressure profile, the effect due to differences of the mass distributions in different redshift bins (see Fig. 1) must be removed, and it needs to be probed by comparing the y-profiles with the ‘same’ mass distribution. The original mass distributions of LRG haloes at the three redshift bins are shown in Fig. 1 in halo mass ($M_{500}$), for which halo masses are estimated by the stellar-to-halo mass relation in Wang et al. (2016) obtained from

gravitational lensing measurements. Tanimura et al. (2019b) shows the y-profile of the LRGs can be described using this stellar-to-halo relation. We find that our constraint results are not quite sensitive to the stellar-to-halo relation. To obtain the same mass distributions in the three redshift bins, we find and select the minimum number of LRGs in a mass bin in the three z-bins (or the overlapped region of the three mass distributions, see Fig. 1), and adopt these LRGs as our sub-sample used in the following discussion. The derived mass distribution for the three z-bins is shown as hatched area in Fig. 1 and it provides 11 660 LRGs for each redshift bin.

Fig. 2 shows the average y-map stacked against 11 660 LRGs at $z_1$, $z_2$, and $z_3$ bin, respectively, with the same mass distribution. The corresponding y-profiles as a function of angular scale $\theta$ with their $1\sigma$ statistical uncertainties are shown in Fig. 3. The uncertainties are estimated by a bootstrap resampling by drawing a random sample of LRGs in a redshift bin. For example, in $z_1$ bin, 11 660 LRGs are resampled from the total of 18 083 LRGs with a replacement allowed. For the mass distribution to be unchanged, we replace an LRG with another LRG of the same mass. We repeat this process 1000 times and the bootstrapped data produce 1000 average y-profiles. The rms fluctuation of the profiles is shown as an uncertainty in Fig. 3. We find that the average y-profiles have central peaks of $y = (3.3 \pm 0.2) \times 10^{-7}$ in $z_1$ bin, $y = (1.7 \pm 0.2) \times 10^{-7}$ in $z_2$ bin, $y = (1.3 \pm 0.2) \times 10^{-7}$ in $z_3$ bin, and $y = (2.1 \pm 0.1) \times 10^{-7}$ for all data summed over the three redshift bins. As can be seen, we detect the tSZ signal out to ~30 arcmin that is well beyond the Planck beam of 10 arcmin. Note that there are correlations between different angular scales when deriving y-profiles. The correlation coefficient matrices in different redshift ranges are shown in Fig. 4. Therefore, instead of using the independent errors shown in Fig. 3, we adopt the full covariance matrix to estimate the $\chi^2$ in the fitting process.

We also notice that, given the same angular scale $\theta$, $y(\theta)$ actually cover different physical scales of clusters at different redshifts. To investigate the stacked tSZ signal at the same cluster physical scales, and as a comparison, we also derive the stacked y-profiles as a function of $r/R_{500}$, where $r$ is the physical scale of radius. The corresponding $y(r/R_{500})$ results are shown in the Appendix. In principle, the results of constraints on the UPP parameters do not depend on whether physical or angular scale is chosen, as long as the same theoretical halo model is used to compare with the same data sets. Since it has great convenience for the $y(\theta)$ when convolving with the beam function (as shown in Sec. 3.2), and considering accuracy and time consumption in the computation process, we choose to use $y(\theta)$ in our following discussion.

# 3 STACKED TSZ PROFILE MODELLING

## 3.1 Compton y-parameter

As shown in equation (3), the Compton y-parameter is an integral of the cluster's pressure profile along the line of sight. Therefore, for projected angle $\theta$ from the centre of the profile, the y-parameter can be calculated as (Komatsu et al. 2011)

$$y(\theta) = \frac{2\sigma_T}{m_e c^2} \int_0^{\sqrt{r_{\text{out}}^2 - \theta^2 D_A^2}} P_e \left( \sqrt{l^2 + \theta^2 D_A^2} \right) dl, \quad (4)$$

where $D_A(z)$ is the angular diameter distance to the LRG sample. The reason that there is a factor of 2 in front of the integral and the lower limit is zero is because of the spherical symmetry. The upper limits of the integral is $\sqrt{r_{\text{out}}^2 - \theta^2 D_A^2}$, where $r_{\text{out}}$ is a truncated
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The stacked y-intensity maps for $z_1$ bin (upper left), $z_2$ bin (upper right), $z_3$ bin (bottom left), and the whole redshift range with all LRG data (bottom right). The number of selected LRGs stacked in each map is 11 660 with the same mass distribution as shown in Fig. 1. The mean tSZ intensity in the annular region between 30 and 40 arcmin has been subtracted as the local background.</figcaption></figure>

scale and we always have $r_{out} \gg \theta D_A$ (Komatsu et al. 2011). The electron UPP $P_e$ can be written as

$$P_e(r) = P'_{500} \mathbb{P}(x), \qquad (5)$$

where $x \equiv r/R_{500}$ and $R_{500}$ is the radius within which the average density is 500 times higher than the critical density of the Universe. The $P'_{500} = P_{500}F_{500}$, where $P_{500}$ is a characteristic pressure, based on the standard self-similar model for the variation of cluster mass (Arnaud et al. 2010)

$$P_{500} = 1.65 \times 10^{-3} E(z)^{8/3} \times \left[ \frac{M_{500}}{3 \times 10^{14} h_{70}^{-1} M_{\odot}} \right]^{2/3} h_{70}^2 \text{ keV cm}^{-3}, \qquad (6)$$

where $h_{70} = h/0.7$ and $E(z) = H(z)/H_0$. The $F_{500}$ is a correction factor of $P_{500}$, reflecting the deviation of the standard self-similar scaling, given by Arnaud et al. (2010) and Planck Collaboration V (2013)

$$F_{500} = \left[ \frac{M_{500}}{3 \times 10^{14} h_{70}^{-1} M_{\odot}} \right]^{0.12}. \qquad (7)$$

The $\mathbb{P}(x)$ in equation (5) is the scaled pressure profile, normalized by $P_{500}$. By adopting the generalized NFW profile (Navarro, Frenk & White 1997; Nagai, Kravtsov & Vikhlinin 2007), it can be characterized as

$$\mathbb{P}(x) = \frac{P_0}{(c_{500}x)^\gamma [1 + (c_{500}x)^\alpha]^{(\beta-\gamma)/\alpha}}, \qquad (8)$$

where $P_0$ is a normalization factor, $c_{500}$ is the concentration factor at $R_{500}$, and $\gamma$, $\alpha$, and $\beta$ are the slopes for the central, intermediate, and outer regions of cluster, respectively. Following Planck Collaboration V (2013), we fix $\gamma = 0.31$, and fit the rest four parameters in the fitting process.

## 3.2 Theoretical stacked tSZ profile

The observed angular cross-correlation function of tSZ signal and galaxy cluster distribution, i.e. stacked tSZ profile, $y_{cross}(\theta)$ can be predicted by expanding it into Legendre polynomials considering beam function of an experiment (see e.g. Komatsu & Kitayama 1999)

$$y_{\text{th}}^{\text{cross}}(\theta) = \frac{1}{4\pi} \sum_{\ell} (2\ell + 1) C_{\ell}^{\text{yc}} P_{\ell}(\cos \theta) B_{\ell}, \qquad (9)$$
<figure><img src="image_4.png" /><figcaption>**Figure 3.** The angular cross-correlations of tSZ signal and galaxy cluster distribution in different redshift ranges. The data points are derived from the stacked intensity maps shown in Fig. 2. The error bars are simply derived from the diagonal elements of the covariance matrix in each case. The solid, dotted, and dashed curves are the best fits of total, one-halo, and two-halo correlation functions, respectively. All of data points and curves are rescaled based on the average values of corresponding intensities between 30 and 40 arcmin.</figcaption></figure>

where $P_ℓ(x)$ are the Legendre polynomials, $C_ℓ^{yc}$ is the angular cross power spectrum of tSZ and cluster distribution, and the beam function $B_ℓ = \exp(-ℓ^2\sigma_b^2/2)$. Here, $\sigma_b = \theta_{FWHM}/\sqrt{8\ln2}$, where we take $\theta_{FWHM} = 10$ arcmin as the Planck beam size. By adopting the flat-sky and Limber approximations, the $C_ℓ^{yc}$ can be expressed as (Cole & Kaiser 1988; Komatsu & Kitayama 1999; Fang, Kadota & Takada 2012)

$$C_{\ell}^{yc} = C_{\ell}^{yc,1h} + C_{\ell}^{yc,2h}, \qquad (10)$$

where $C_{\ell}^{yc,1h}$ and $C_{\ell}^{yc,2h}$ are the one-halo and two-halo terms, respectively.

Theoretically, following Fang et al. (2012), in given redshift and cluster mass ranges, the one- and two-halo terms take the forms as

$$C_{\ell}^{yc,1h} = \frac{1}{\bar{n}_{2D}} \int_{z_1}^{z_u} dz \frac{c \chi^2}{H(z)} \int_{M_1}^{M_u} dM \frac{dn}{dM} y_{\ell}(M, z), \qquad (11)$$

$$C_{\ell}^{yc,2h} = \frac{1}{\bar{n}_{2D}} \int_{z_1}^{z_u} dz \frac{c \chi^2}{H(z)} P_m(k, z) W^c(z) W_{\ell}^y(z). \qquad (12)$$

Here, $z_1$ and $z_u$, $M_1$ and $M_u$ are the lower and upper bounds of redshift and mass ranges, respectively. $\chi$ is the comoving distance, $b(M, z)$ and $dn/dM$ are the halo bias and mass function (Sheth & Tormen 1999), and $P_m(k, z)$ is the linear matter power spectrum, where $k = (l + 1/2)/\chi$. The $\bar{n}_{2D}$ is the 2D angular number density
<figure><img src="image_5.png" /><figcaption>**Figure 4.** The correlation coefficient matrices between angular scales for the y-profiles $y(\theta)$ in different redshift ranges. The angular scale $\theta$ has been divided into 15 bins with $\Delta\theta = 2$ arcmin.</figcaption></figure>

of galaxy cluster in given redshift and cluster mass ranges that can be estimated by (Fang et al. 2012)

$$ \bar{n}_{2D} = \int_{z_1}^{z_u} dz \frac{c\chi^2}{H(z)} \int_{M_1}^{M_u} dM \frac{dn}{dM}(M, z). \quad (13) $$

The functions $W^c(z)$ and $W^y_\ell(z)$ are defined as

$$ \begin{aligned} W^c(z) &= \int_{M_1}^{M_u} dM \frac{dn}{dM}(M, z) b(M, z), \\ W^y_\ell(z) &= \int dM \frac{dn}{dM}(M, z) b(M, z) y_\ell(M, z), \end{aligned} \quad (14) $$

where the second integral should cover all possible halo masses as it is the influence from all other correlated haloes. The $y_\ell(M, z)$ is the Fourier transform of the Compton y-parameter, written as

$$ y_{\ell}(M, z) = \frac{a}{\chi^2} \left( \frac{\sigma_{\text{T}}}{m_e c^2} \right) u_{\text{p}} \left( k = \frac{\ell + 1/2}{\chi} \mid M, z \right), \quad (15) $$

and

$$ u_{\text{p}} = \int dr \, 4\pi r^2 \frac{\sin(kr)}{kr} P_{\ell}(r|M, z). \quad (16) $$

Note that the $P_\ell(r|M, z)$ here is based on virial mass $M_{\text{vir}}$ instead of $M_{500}$ used in equation (5), so we need to calculate this $P_\ell$ by converting $M_{\text{vir}}$ to $M_{500}$ (e.g. see appendix B in Planck Collaboration LIII 2017). Then, we can estimate the angular cross-power spectrum $C_\ell^{\text{yc}}$ theoretically using equations (11)–(16) with the help of halo model (Cooray & Sheth 2002).

## 3.3 Predicted stacked tSZ profile

### 3.3.1 Separating into mass and redshift bins

Since we have large number of observed LRGs in SDSS DR7 (in total 74 681, and 18 083, 29 586, and 27 012 samples in $z_1, z_2$, and $z_3$ redshift bins, respectively), it is quite time-consuming to calculate $\bar{y}(\theta)$ in a redshift range. In order to speed up the calculation, we further divide each redshift bin ($z_1, z_2$, and $z_3$) into four sub-bins and divide the total mass range of clusters into 10 different mass bins.

We compute the average redshift $\bar{z}_i$ for a redshift bin $i$, and average cluster mass $\bar{M}_{500}^j$ and radius $\bar{R}_{500}^j$ for a mass bin $j$, and estimate the $P_{e,i}^{(j)}(r|\bar{M}_{500}^j, \bar{R}_{500}^j)$ and $y_i^{(j)}(\theta|\bar{M}_{500}^j, \bar{R}_{500}^j)$ using equations (5) and (3) for redshift bin $i$ and mass bin $j$ (below, we shorten it as 'the $ij$th bin'). We also count the number of samples $n_i^{(j)}$ within the $ij$th bin. Then the mean Compton y-parameter can be calculated as

$$ \bar{y}(\theta) = \frac{1}{N_c} \sum_{i}^{N_z} \sum_{j}^{N_M} f_M^{(j)} n_i^{(j)} y_i^{(j)}(\theta). \quad (17) $$

Here, $N_z = 4$ is the number of redshift sub-bins for either $z_1, z_2$, or $z_3$ bin ($N_z = 12$ for the whole sample). $N_M = 10$ is the number of cluster mass bins. $f_M^{(j)}$ is the number fraction of the selected sample in a mass bin. $N_c$ is the total number of clusters, given by

$$ N_c = \sum_{i}^{N_z} \sum_{j}^{N_M} n_i^{(j)}. \quad (18) $$

Using such a binning strategy, the computation speed becomes faster since we do not need to calculate the $O(10^4)$ number of $y(\theta)$ profile individually. To verify the accuracy of such approach, we compute the $\bar{y}(\theta)$ for such binning strategy and the individual sum up, and the difference is only at $\sim 1$ per cent at most in different redshift ranges. Hence, we will use this scheme to compute the one-halo term as follows (see equation 19).

### 3.3.2 One-halo term

Since we have the information of LRG redshift $z$, $M_{500}$, and $R_{500}$, we can directly predict the one-halo term for each cluster, instead of using the halo mass function $dn/dM$ in equation (11).

We first count the number of LRGs for each sub-bin of mass and redshift, i.e. counting $n_i^{(j)}$, and then calculate the $y_i^{(j)}(\theta)$ profile for each sub-bin. Then we use equation (17) to calculate the average of the profile within each redshift bin ($z_1, z_2, z_3$ or the whole redshift range). This profile is the unconvolved, averaged profile for a certain redshift bin. Because for real data, it is the original Compton y-profile convolved with Planck beam, so we need to calculate the convolved averaged profile with Planck beam function. For this reason, we first Fourier transform the averaged profile into the $\ell$-space of one-halo term $C_\ell^{\text{yc,1h}}$, i.e.

$$ C_{\ell}^{\text{yc,1h}} = 2\pi \int_{-1}^{1} P_{\ell}(\cos\theta) \bar{y}(\theta) \, d\cos\theta, \quad (19) $$

where $P_\ell$ are the Legendre polynomials. Then we calculate the following two-halo term by using halo model.

### 3.3.3 Two-halo term

Similarly, making use of the same binning method, we can also calculate the two-halo term $C_\ell^{\text{yc,2h}}$, and we have

$$ C_{\ell}^{\text{yc,2h}} = \frac{1}{N_c} \sum_{i}^{N_z} \sum_{j}^{N_M} f_M^{(j)} n_i^{(j)} C_{\ell}^{\text{yc,2h}}(\bar{M}_j, \bar{z}_i), \quad (20) $$

where $C_\ell^{\text{yc,2h}}(\bar{M}_j, \bar{z}_i)$ is the two-halo term in the mass bin $j$ and redshift bin $i$. Note that $\bar{M}_j$ is an average virial mass, and can be obtained by $\bar{M}_{500}^j$. The $C_\ell^{\text{yc,2h}}$ can be computed using equation (12) considering binning. In each redshift and mass bin, we notice that $\int_{z_1}^{z_u} dz \to \Delta z_i$ and $\int_{M_1}^{M_u} dM \to \Delta M_j$, where $\Delta z_i$ and $\Delta M_j$ are
<table><caption>**Table 1.** The best fits and $1\sigma$ errors of the free parameters in the UPP model. The parameter $\gamma$ is fixed to be 0.31 in the fitting process (Planck Collaboration V 2013). The results from Arnaud et al. (2010, A10) and Planck Collaboration V (2013, Planck13) are also shown at the bottom as comparison.</caption><thead><tr><th>Case</th><th>P<sub>0</sub></th><th>c<sub>500</sub></th><th>α</th><th>β</th><th>γ</th><th>η</th><th>χ<sup>2</sup><sub>red</sub></th></tr></thead><tbody><tr><td>z<sub>1</sub> bin</td><td>3.35<sup>+8.18</sup><sub>-1.35</sub></td><td>1.45<sup>+0.56</sup><sub>-0.41</sub></td><td>4.31<sup>+0.66</sup><sub>-2.88</sub></td><td>5.14<sup>+3.00</sup><sub>-0.86</sub></td><td>0.31</td><td>—</td><td>1.65</td></tr><tr><td>z<sub>2</sub> bin</td><td>10.46<sup>+1.54</sup><sub>-6.85</sub></td><td>1.88<sup>+0.53</sup><sub>-0.55</sub></td><td>4.03<sup>+0.96</sup><sub>-2.79</sub></td><td>4.39<sup>+3.45</sup><sub>-0.79</sub></td><td>0.31</td><td>—</td><td>0.45</td></tr><tr><td>z<sub>3</sub> bin</td><td>6.15<sup>+5.25</sup><sub>-3.67</sub></td><td>2.23<sup>+2.42</sup><sub>-0.80</sub></td><td>1.24<sup>+3.68</sup><sub>-0.50</sub></td><td>3.00<sup>+0.43</sup><sub>-0.30</sub></td><td>0.31</td><td>—</td><td>0.52</td></tr><tr><td>All data</td><td>2.18<sup>+9.02</sup><sub>-1.98</sub></td><td>1.05<sup>+1.27</sup><sub>-0.47</sub></td><td>1.52<sup>+1.47</sup><sub>-0.58</sub></td><td>3.91<sup>+0.87</sup><sub>-0.44</sub></td><td>0.31</td><td>—</td><td>1.88</td></tr><tr><td>3 bins</td><td>3.04<sup>+15.53</sup><sub>-1.93</sub></td><td>2.12<sup>+0.82</sup><sub>-1.01</sub></td><td>5.80<sup>+3.89</sup><sub>-4.46</sub></td><td>4.83<sup>+2.20</sup><sub>-0.78</sub></td><td>0.31</td><td>—</td><td>1.61</td></tr><tr><td>3 bins with η</td><td>2.99<sup>+3.44</sup><sub>-1.57</sub></td><td>1.16<sup>+0.79</sup><sub>-0.29</sub></td><td>2.66<sup>+1.67</sup><sub>-0.97</sub></td><td>5.48<sup>+2.39</sup><sub>-1.38</sub></td><td>0.31</td><td>-3.11<sup>+1.09</sup><sub>-1.13</sub></td><td>1.39</td></tr><tr><td>A10</td><td>8.403 h<sub>70</sub><sup>-3/2</sup></td><td>1.177</td><td>1.0510</td><td>5.4905</td><td>0.3081</td><td>—</td><td>—</td></tr><tr><td>Planck13</td><td>6.41</td><td>1.81</td><td>1.33</td><td>4.13</td><td>0.31</td><td>—</td><td>—</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** The 1D PDFs of $P_0$, $c_{500}$, $\alpha$, and $\beta$ in different redshift ranges. The blue dash-dotted, green dashed, and red dotted curves are the results in $z_1$, $z_2$, and $z_3$ bins, respectively. The solid purple curve is for the whole redshift range with all stacked data.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Left: The fitting result of $y^{cross}(\theta)$ by simultaneously fitting the three redshift bins without parameter $\eta$. Right: The same as left-hand panel but $\eta$ included. As can be seen, the data cannot be well fitted within $\theta < 5$ arcmin in the $z_1$ and $z_3$ bins in the left-hand panel, and it can be improved by including $\eta$ as shown in the right-hand panel.</figcaption></figure>

the bin widths of the redshift bin $i$ and mass bin $j$, respectively. Therefore, we can simplify the equation as

$$C_{\ell}^{\text{yc,2h}}(\bar{M}_j, \bar{z}_i) = b(\bar{M}_j, \bar{z}_i)P_m\left(k = \frac{\ell + 1/2}{\chi(\bar{z}_i)}, \bar{z}_i\right) W_{\ell}^{y}(\bar{z}_i). \quad (21)$$

Therefore, the predicated stacked tSZ profile can be calculated by adding together equation (20) with equation (19) (i.e. equation 10), and then Fourier transform it back to real space with multiplication of Planck beam function through equation (9). Besides, in Section 2.3, we showed that in order to subtract background we subtract the mean value of the measured profile in the range of 30–40 arcmin. Therefore, in order to make accurate comparison, we also need to subtract the mean value for theoretical profile, i.e.

$$\tilde{y}_{\text{th}}^{\text{cross}}(\theta) = y_{\text{th}}^{\text{cross}}(\theta) - \tilde{y}_{\text{th}}^{\text{cross}}|_{\text{BG}}, \quad (22)$$

where $\tilde{y}_{\text{th}}^{\text{cross}}|_{\text{BG}}$ is the mean value of $y_{\text{th}}^{\text{cross}}(\theta)$ in the range of 30–40 arcmin as the background value. We calculate $\tilde{y}_{\text{th}}^{\text{cross}}(\theta)$ in the whole redshift range and three main redshift bins, respectively, for sampling of the UPP parameters ($P_0, c_{500}, \alpha, \beta$), and fit the corresponding data sets to derive the probability distribution functions (PDFs) of the these parameters.

## 4 FITTING METHOD

In order to fit the measured stacked tSZ profile $y_{\text{obs}}^{\text{cross}}(\theta)$ and extract the parameter values of the UPP $P_e(r)$, we adopt the $\chi^2$ statistic method, given by

$$\chi^2 = \sum_{m,n}^{N_\theta} \left[ y_{\text{obs}}^{\text{cross}}(\theta_m) - \tilde{y}_{\text{th}}^{\text{cross}}(\theta_m) \right] (C_{\theta\theta'})^{-1}_{mn} \times \left[ y_{\text{obs}}^{\text{cross}}(\theta_n) - \tilde{y}_{\text{th}}^{\text{cross}}(\theta_n) \right], \quad (23)$$

where $N_\theta$ is the number of data at different angular scales, and $C_{\theta\theta'}$ is the covariance matrix. Then we can calculate the likelihood function as $\mathcal{L} \sim \exp(-\chi^2/2)$.

We employ the Markov Chain Monte Carlo (MCMC) method to perform constraints on the free parameters in the UPP. The

Metropolis-Hastings algorithm is adopted to determine the accepting probability of the new chain points (Metropolis et al. 1953; Hastings 1970). The proposal density matrix is obtained by a Gaussian sampler with adaptive step size (Doran & Muller 2004). We assume uniform flat priors for all the free parameters, and their ranges are set to be $P_0 \in (0, 20)$, $c_{500} \in (0, 10)$, $\alpha \in (0, 10)$, and $\beta \in (0, 10)$. We also add a free parameter $\eta$ on the power-law index of $E(z)$ in equation (6), when simultaneously fitting the three redshift bins together. This parameter can adjust the redshift-dependence of the electron pressure profile, and we set $\eta \in (-10, 10)$. We run fifteen parallel chains for each case we explore, and obtain $10^5$ chain points for one chain after it reaches the convergence criterion (Gelman & Rubin 1992). After performing the burn-in and thinning the chains, we merge all chains together and obtain about 10 000 chain points to illustrate one-dimensional (1D) and 2D PDFs of the free parameters.

## 5 RESULTS OF CONSTRAINTS

In Fig. 3, we show the fitting results for different redshift ranges. The solid, dotted, and dashed curves are the best fits of total, one-halo, and two-halo correlation functions, respectively. Note that we use the covariance matrix between different angular scales instead of the error bars shown in Fig. 3 in the fitting process. These error bars are simply derived from the diagonal elements of the corresponding covariance matrix. As shown in Table 1, We find that the minimum reduced chi-square, which is defined by $\chi^2_{\text{red}} = \chi^2_{\text{min}}/(N - M)$, where N and M are the number of data and free parameters in the model for the three main redshift bins and the whole range with all LRG data are 1.65, 0.45, 0.52, and 1.88, respectively. This indicates that the data can be well fitted in each case, especially for the $z_2$ and $z_3$ bins.

The corresponding best fits and $1\sigma$ errors of $P_0, c_{500}, \alpha$, and $\beta$ are shown in Table 1, and the 1D PDFs are given in Fig. 5. We find that the $P_0$ is not well constrained by the data, and its 1D PDFs extend in large parameter space with a wide peak between 2 and 11 in different redshift ranges. The constraint result of $c_{500}$ in each case is in a good agreement that the probability peaks
<figure><img src="image_8.png" /><figcaption>**Figure 7.** The contour maps and 1D PDF of free parameters when fitting the data of all three bins with and without $\eta$. The 68.3 per cent, 95.5 per cent, and 99.7 per cent confidence levels are shown. We can find that the best fit of $\eta = -3.11^{+1.09}_{-1.13}$, which significantly changes the power-law index 8/3 of $E(z)$ in equation (6). This indicates that our result prefers weaker redshift dependence of the electron pressure profile than Arnaud et al. (2010).</figcaption></figure>

are around 1.5, although the width of the PDF of the $z_3$ bin is wider than others. We do not obtain stringent constraints on $\alpha$ in the three z-bins separately. As can be seen, the PDFs have flat tops extending from 1.5 to 4.5 for the $z_1$, $z_2$, and $z_3$ bins. The constraint is significantly improved for the all-data case with a peak at 1.5. For $\beta$, we find that the results of the $z_2$ and $z_3$ bins are well consistent, while the best-fitting value is apparently smaller in the $z_3$ bin.

We also find that, generally, the fitting results of the four cases are consistent with that given by Arnaud et al. (2010) and Planck Collaboration V (2013), especially for the all-data case (e.g. the constraints on $c_{500}$, $\alpha$, and $\beta$). This implies that our method is feasible and effective for the studies of the cluster electron pressure profile. Beside, as indicated in Table 1, there is noticeable evolution of the parameters, i.e. $c_{500}$, $\alpha$, and $\beta$, in the three redshift bins. We can see that $\alpha$, and $\beta$ become smaller and smaller as the redshift increases, while $c_{500}$ tends to be larger at high redshift.

In order to suppress the fitting uncertainties of the free parameters and check the consistency of the UPP model in different redshift bins, we also try to simultaneously fit the data in the three redshift bins using $\chi^2_{3\text{bins}} = \chi^2_{z1} + \chi^2_{z2} + \chi^2_{z3}$. As shown in the left-hand panel of Fig. 6, we find that the data can be fitted with $\chi^2_{\text{red}} = 1.61$ (see Table 1), adopting the usual model of electron UPP with four free parameters, i.e. $P_0$, $c_{500}$, $\alpha$, and $\beta$.<fnref n="4" /> However, we can see that, except for the $z_2$ bin, the theoretical curves cannot fit the data very well within angular scales less than 5 arcmin in the $z_1$ and $z_3$ bins.

Since we can find that, by comparing to the data, the predicated curve is lower in the $z_1$ bin, and higher in the $z_3$ bin, we can adjust the redshift dependence of the model by adding a free parameter $\eta$ on the power-law index of $E(z)$ in equation (6) (i.e. $E(z)^{8/3} \rightarrow E(z)^{8/3 + \eta}$). In the right-hand panel of Fig. 6, we show the fitting result when including $\eta$. We can find that the fitting results are

<footnote n="4">⁴We also check the fitting result by setting $\gamma$ as a free parameter, and the result is almost the same (but with wider PDFs).</footnote>
significantly improved with $\chi^2_{\text{red}} = 1.39$, that is ~5 smaller for the $\chi^2_{\text{min}}$ than the case without $\eta$ (see Table 1). We find that the best fit of $\eta$ is $-3.11_{-1.13}^{+1.05}$, which significantly changes the previous power-law index (i.e. 8/3) of $E(z)$. This indicates that our result prefers weaker redshift dependence for the cluster gas pressure model.

The 1D PDFs and 2D contour maps of the free parameters for the three bins with and without $\eta$ are shown in Fig. 7. The $1\sigma$ (68.3 per cent), $2\sigma$ (95.5 per cent), and $3\sigma$ (99.7 per cent) confidence levels are shown here. We can find that the fitting results of the three bins with and without $\eta$ are basically consistent with each other. Besides, by comparing to fig. 5 in Planck Collaboration V (2013), we find that our constraint results (contours and 1D PDFs) are in a good agreement with their results, and our method can even offer more stringent constraints on $c_{500}$ and $\beta$. The constraint can be further improved in the future with more LRG sample included.

# 6 SUMMARY AND DISCUSSION

In this work, we make use of the cross-correlation between the tSZ signal measured by Planck satellite and LRGs from SDSS DR7 to study the UPP of galaxy clusters. We first stack the Planck y-map against the LRGs and derive the mean y-profile. The LRG sample is given by SDSS DR7 with $0.16 < z < 0.47$ and $3 \times 10^{12} \lesssim M_{500} \lesssim 3 \times 10^{14} M_{\odot}$. In order to study the redshift evolution of the properties of intra-cluster gas, we divide the LRG sample into sub-samples in three redshift bins. We find that the peaks of the mass distributions in the three redshift bins move towards higher halo mass as the redshift increases. In order to remove the effect of non-matching mass distribution in our analysis, we select sub-sample from each redshift bin with the same mass distribution.

Then we derive theoretical stacked tSZ y-profile with the help of the halo model. The Planck beam function is also considered in the estimates. To obtain more realistic and accurate predictions, we take into account of the information from measurements, such as the redshift, $M_{500}$, and $R_{500}$, and simplify the calculation by dividing the LRG sample into sub-redshift and mass bins. We adopt the MCMC technique to illustrate the probability distribution of the parameters in the gas pressure profile, and set wide parameter ranges as prior distributions.

We separately fit the y-profile data obtained from the stacked y-map in the $z_1$, $z_2$, and $z_3$ redshift bins, and the whole redshift range with all stacking data. We find that the fitting results of the four parameters, i.e. $P_0$, $c_{500}$, $\alpha$, and $\beta$, are mainly consistent with one another. They are also in a good agreement with the results from Arnaud et al. (2010) and Planck Collaboration V (2013), especially for the all-data case. We find that there is evolution for $c_{500}$, $\alpha$, and $\beta$ from the low to high redshift that the best fits of $\alpha$, and $\beta$ become smaller, while $c_{500}$ becomes greater as the redshift increases.

In order to further investigate the redshift evolution and check the consistence of the UPP model at different redshift bins, we fit the data in the three redshift bins simultaneously by summing up the $\chi^2$ of the three bins together. Interestingly, we find that the UPP model cannot provide good fits on the data at $\theta < 5$ arcmin in the $z_1$ and $z_3$ bins. After checking the results, we propose to add a parameter $\eta$ on the power-law index of $E(z)$ to change the redshift dependence of the model. The best-fitting value of $\eta$ is $-3.11_{-1.13}^{+1.09}$, which suggests that the power index of redshift evolution of the UPP profile should be equal to $\eta + 8/3 = -0.44_{-1.13}^{+1.09}$. This result implies that the UPP profile may be less redshift dependent than its original form. Physically, this result indicates that the cluster pressure profile is less evolved than it was thought to be, and Compton y-profile

could be nearly redshift independent. By including this factor, we find the $\chi^2_{\text{red}}$ decreases from 1.61 to 1.39, which is ~5 smaller in $\chi^2_{\text{min}}$ than that given by the usual UPP model. By comparing the 1D and 2D PDFs with Planck Collaboration V (2013), we find our results from the three bins with and without $\eta$ cases can match theirs very well, and can even offer more stringent constraints on $c_{500}$ and $\beta$. This indicates that our method can provide reliable results that prefer less redshift dependence of the UPP mode. We will include more samples and further confirm our result in the future work, e.g. analysing the SDSS DR12 data, and provide more accurate constraints on the cluster pressure profile.

Besides, this study also can be an important step towards fully quantify the distribution of missing baryons. This is because, a significant amount of baryons are associated with filaments, voids and sheets which have much weaker signals of SZ effect than haloes (Haider et al. 2016). As one can see in Tanimura et al. (2019a), the stacked SZ signal of the filaments is usually entangled with the halo contribution. Therefore, improvement on the halo model's pressure profile will lead to a more precise subtraction of the halo contribution, and will result in better measurement of the signals from filaments and sheets. Although this is out of the scope of this paper, our study can eventually contribute to the more precise determination of the baryons within filaments and sheets.

# ACKNOWLEDGEMENTS

YG acknowledges the support of NSFC-11822305, NSFC-11773031, NSFC-11633004, the Chinese Academy of Sciences (CAS) Strategic Priority Research Program XDA15020200, the NSFC-ISF joint research programme no. 11761141012, and CAS Interdisciplinary Innovation Team. YZM is supported by the National Research Foundation of South Africa with grant no. 105925, no. 110984, and NSFC-11828301. This research has been also supported by funds from the European Research Council (ERC) under the Horizon 2020 research and innovation programme grant agreement of the European Union: ERC-2015-AdG 695561 (ByoPiC, https://byopic.eu).

# REFERENCES

Anderson L. et al., 2014, MNRAS, 441, 24

Arnaud M., Pratt G. W., Piffaretti R., Böhringer H., Croston J. H., Pointecouteau E., 2010, A&A, 517, A92

Bennett C. L. et al., 2003, ApJS, 143, 1

Birkinshaw M., 1999, Phys. Rep., 310, 97

Birkinshaw M., Gull S. F., 1978, Nature, 274, 111

Blanton M. R., Roweis S., 2007, AJ, 133, 734

Bruzual G., Charlot S., 2003, MNRAS, 344, 1000

Carlstrom J. E., Holder G. P., Reese E. D., 2002, Ann. Rev. Astron. Astrophys., 40, 643

Chabrier G., 2003, PASP, 115, 763

Cole S., Kaiser N., 1988, MNRAS, 233, 637

Cooray A., Sheth R., 2002, Phys. Rep., 372, 1

Doran M., Muller C. M., 2004, J. Cosmol. Astropart. Phys., 09, 003

Eisenstein D. J. et al., 2005, ApJ, 633, 560

Fang W., Kadota K., Takada M., 2012, PRD, 85, 023007

Gelman A., Rubin D., 1992, Stat. Sci., 7, 457

Gorski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759

Haider M., Steinhauser D., Vogelsberger M., Genel S., Springel V., Torrey P., Hernquist L., 2016, MNRAS, 457, 3024

Hastings W. K., 1970, Biometrika, 57, 97

Hojjati A. et al., 2017, MNRAS, 471, 1565

Hoshino H. et al., 2015, MNRAS, 452, 998
Hurier G., Macias-Perez J. F., Hildebrandt S. R., 2013, A&A, 558, 118

Kazin E. A. et al., 2010, ApJ, 710, 1444

Komatsu E. et al., 2011, ApJS, 192, 18

Komatsu E., Kitayama T., 1999, ApJ, 526, L1–L4

Le Brun A. M. C., McCarthy I. G., Melin J.-B., 2015, MNRAS, 451, 3868

Ma Y.-Z., Van Waerbeke L., Hinshaw G., Hojjati A., Scott D., Zuntz J., 2015, J. Cosmol. Astropart. Phys., 9, 046

McCarthy I. G. et al., 2007, MNRAS, 376, 497

Metropolis N., Rosenbluth A. W., Rosenbluth M. N., Teller A. H., Teller E., 1953, J. Chem. Phys., 21, 1087

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

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

Planck Collaboration VIII, 2011, A&A, 536, A8

Planck Collaboration V, 2013, A&A, 550, A131

Planck Collaboration XI, 2013, A&A, 557, A52

Planck Collaboration XVI, 2014, A&A, 571, A16

Planck Collaboration I, 2016, A&A, 594, A1

Planck Collaboration XXII, 2016, A&A, 594, A22

Planck Collaboration LIII, 2017, A&A, 617, A48

Remazeilles M., Aghanim N., Douspis M., 2013, MNRAS, 430, 370

Sheth R. K., Tormen G., 1999, MNRAS, 308, 119

Sunyaev R. A., Zeldovich Y. B., 1972, Comments Astrophys. Space Phys., 4, 173

Sunyaev R. A., Zeldovich I. B., 1980, MNRAS, 190, 413

Tanimura H. et al., 2019a, MNRAS, 483, 223

Tanimura H. et al., 2019b, preprint (arXiv:1903.06654)

Wang W., White S. D. M., Mandelbaum R., Henriques B., Anderson M. E., Han J., 2016, MNRAS, 456, 2301

## APPENDIX: STACKED TSZ PROFILE IN PHYSICAL SCALES

In order to stack the tSZ signal at the same physical scales from each LRG, we stack the *Planck* y-map at the positions of LRGs in physical coordinates, instead of angular coordinate. We reset a 2D physical coordinate system of $-20 < r/R_{500} < 20$ and $-20 < r/R_{500} < 20$, divided in $80 \times 80$ bins. The local background region is also redefined in physical scale, which is an annular region between $|r/R_{500}| = 15$ and $20$ for each LRG. Then, we follow the same procedure as described in Section 2.3 for stacking. After removing the LRGs with $\ge 20$ per cent masked region within a $r/R_{500} = 20$ circle, we obtain $18\,117, 29\,074$, and $26\,401$ LRGs, respectively, for the $z_1, z_2, z_3$ bins, and total $73\,592$ in the whole redshift range. We finally select $11\,926$ LRGs in each redshift bin with the same mass distribution. This mass distribution is quite similar with that in the $y(\theta)$ case.

The y-maps, cross-correlations of tSZ signal and galaxy cluster distribution, and correlation coefficient matrices between physical scales in the $y(r/R_{500})$ case are shown in Figs A1–A3, respectively. By comparing to the $y(\theta)$ case, as shown in Figs 3 and A2, we find that there is a bump feature around $r/R_{500} \sim 10$ in the $y(r/R_{500})$ profiles, especially for the $z_1$ bin, which is mainly due to the two-halo term. This feature is smoothed out in the $y(\theta)$, since it is obtained by stacking different cluster physical scales at the same angular scale $\theta$.

<figure><img src="image_9.png" /><figcaption>**Figure A1.** The stacked y-intensity maps for $z_1$ bin (upper left), $z_2$ bin (upper right), $z_3$ bin (bottom left), and the whole redshift range with all LRG data (bottom right). The number of selected LRGs stacked in each map is 11 926 with the same mass distribution. The mean tSZ intensity in the annular region between $r/R_{500} = 15$ and $20$ has been subtracted as the local background.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure A2.** The cross-correlations of tSZ signal and galaxy cluster distribution as a function of $r/R_{500}$ in different redshift ranges. The error bars are simply derived from the diagonal elements of the covariance matrix in each case. All of data points are rescaled based on the average values of corresponding intensities between $r/R_{500}$ = 15 and 20.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure A3.** The correlation coefficient matrices between physical scales for the y-profiles y(r/R<sub>500</sub>) in different redshift ranges. The scale r/R<sub>500</sub> has been divided into 15 bins.</figcaption></figure>

This paper has been typeset from a TeX/LA TeX file prepared by the author.