<?xml-model href='http://www.tei-c.org/release/xml/tei/custom/schema/relaxng/tei_all.rng' schematypens='http://relaxng.org/ns/structure/1.0'?><TEI xmlns="http://www.tei-c.org/ns/1.0">
	<teiHeader>
		<fileDesc>
			<titleStmt><title level='a'>Evolution of the Thermodynamic Properties of Clusters of Galaxies out to Redshift of 1.8</title></titleStmt>
			<publicationStmt>
				<publisher></publisher>
				<date>03/01/2021</date>
			</publicationStmt>
			<sourceDesc>
				<bibl> 
					<idno type="par_id">10232678</idno>
					<idno type="doi">10.3847/1538-4357/abc68d</idno>
					<title level='j'>The Astrophysical Journal</title>
<idno>0004-637X</idno>
<biblScope unit="volume">910</biblScope>
<biblScope unit="issue">1</biblScope>					

					<author>Vittorio Ghirardini</author><author>Esra Bulbul</author><author>Ralph Kraft</author><author>Matt Bayliss</author><author>Bradford Benson</author><author>Lindsey Bleem</author><author>Sebastian Bocquet</author><author>Micheal Calzadilla</author><author>Dominique Eckert</author><author>William Forman</author><author>Juan David Da González</author><author>Gourav Khullar</author><author>Guillaume Mahler</author><author>Michael McDonald</author>
				</bibl>
			</sourceDesc>
		</fileDesc>
		<profileDesc>
			<abstract><ab><![CDATA[The thermodynamic properties of the hot plasma in galaxy clusters retain information on the processes leading to the formation and evolution of the gas in their deep, dark matter potential wells. These processes are dictated not only by gravity but also by gas physics, e.g., active galactic nucleus feedback and turbulence. In this work, we study the thermodynamic properties, e.g., density, temperature, pressure, and entropy, of the most massive and the most distant (seven clusters at z>1.2) clusters selected by the South Pole Telescope and compare them with those of the nearby clusters (13 clusters at z<0.1) to constrain their evolution as a function of time and radius. We find that thermodynamic properties in the outskirts of high-redshift clusters are remarkably similar to the low-redshift clusters,and their evolution follows the prediction of the self-similar model. Their intrinsic scatteris larger, indicating that the physical properties that lead to the formation and virialization of cluster outskirts show evolving variance.On the other hand, thermodynamic properties in the cluster cores deviate significantly from selfsimilarity, indicating that the processes that regulate the core are already in place in these very high redshift clusters.This result is supported by the unevolving physical scatter of all thermodynamic quantities in cluster cores.Unified Astronomy Thesaurus concepts: Galaxy clusters (584); Intracluster medium (858); Galactic and extragalactic astronomy (563); High-redshift galaxy clusters (2007)]]></ab></abstract>
		</profileDesc>
	</teiHeader>
	<text><body xmlns="http://www.tei-c.org/ns/1.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:xlink="http://www.w3.org/1999/xlink">
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="1.">Introduction</head><p>Clusters of galaxies are the largest gravitationally bound objects in the universe and are ideal laboratories to study how cosmic structures form and evolve in time. While the majority of their mass is in the form of dark matter, the hot fully ionized plasma, i.e., the intracluster medium (ICM), retains most of the baryonic component, with only a small contribution from stars and cold gas (3%-5%; <ref type="bibr">Gonzalez et al. 2013</ref>). The ICM is observable in the X-ray band mainly through its emission via thermal bremsstrahlung and radiative recombination processes. X-ray observationsof clustersof galaxies provide in-depth information aboutthe ICM's thermodynamic properties. The thermal Sunyaev-Zeldovich (SZ) effect, a spectral distortion of the cosmic microwave background caused by the ICM, provides a complementary tool for finding clusters at all redshifts and examining their properties.</p><p>X-ray studies of clusters of galaxies provided constraints on thermodynamic properties of the ICM in nearby clusters with redshifts of <ref type="bibr">&lt;0.3 (e.g., De Grandi &amp; Molendi 2002;</ref><ref type="bibr">Croston et al. 2006;</ref><ref type="bibr">Vikhlinin et al. 2006;</ref><ref type="bibr">Cavagnolo et al. 2009;</ref><ref type="bibr">Arnaud et al. 2010;</ref><ref type="bibr">Pratt et al. 2010;</ref><ref type="bibr">Bulbul et al. 2012</ref>). X-ray observations have also provided the serendipitous detection of single high-redshift clusters (z&#61600;&gt;&#61600;1; <ref type="bibr">&#61600; Fabian et al. 2003;</ref><ref type="bibr">Tozzi et al. 2015;</ref><ref type="bibr">Brodwin et al. 2016)</ref>; however,these studies are prone to X-ray selection biases (e.g., the cool-core bias; <ref type="bibr">Eckert et al. 2011</ref>). The majority of theoretical studies in the literature also focus on predicting thermodynamic properties of the ICM in nearby clusters <ref type="bibr">(Kravtsov &amp; Borgani 2012)</ref>. In recent years, owing to the wide-area sky surveys performed with the current SZ telescopes, e.g., the South Pole Telescope (SPT; <ref type="bibr">Carlstrom et al. 2011)</ref>, the Atacama Cosmology Telescope <ref type="bibr">(Fowler et al. 2007)</ref>, and the Planck mission <ref type="bibr">(Planck Collaboration etal. 2016)</ref>, it has become possible to detect clusters outto much higher redshifts (z&#61600;&#8764;&#61600;1.8) with a simpler selection function, i.e., the SZ signaltightly correlates with mass (Planck <ref type="bibr">Collaboration et al. 2014;</ref><ref type="bibr">Bocquet et al. 2019)</ref>. Therefore, X-ray followup observations of the SZ-selected clusters provide a unique opportunity to study the evolution of ICM properties in a uniform way.</p><p>Integrated X-ray properties of the SPT-selected clusters spanning a large redshift range have been studied in the literature <ref type="bibr">(McDonald et al. 2014;</ref><ref type="bibr">Sanders et al. 2018;</ref><ref type="bibr">Bulbul et al. 2019)</ref>. <ref type="bibr">Bartalucci et al. (2017a</ref><ref type="bibr">Bartalucci et al. ( , 2017b) )</ref> examined the individual thermodynamic properties of the ICM by combining the Chandra and XMM-Newton followup observationsof a handful of high-redshiftclusters(z&#8764; 1) detected by SPT and ACT. Studies of the evolution of the ICM propertiesin large SZ-selected cluster sampleshave become possiblewith large targeted X-ray follow-up programs, e.g., Chandra Large Program (LP). <ref type="bibr">McDonald et al. (2013</ref><ref type="bibr">McDonald et al. ( , 2014) )</ref> studied the stacked thermodynamic properties of SPT-selected clusters in a large redshift range, from 0.3 to 1.2, and in particular reported thatthe evolution in the electron numberdensity is consistent with the self-similarexpectation, where only gravitational forces dominate the formation and evolution of the ICM in the intermediate regions (0.15R 500 -R 500 ) 13 of the SPT-selected clusters of galaxies in the redshift range of 0.2&#61600;&lt;&#61600;z&#61600;&lt;&#61600;1.2. The authors also found a clear deviation from self-similarity in the evolution of the core density of these clusters. Deeper Chandra observations of eight high-redshift SPT-selectedclusters beyond a redshift of 1.2 confirm earlier results of no evolution in the cluster cores, indicating that active galactic nucleus (AGN) feedback is tightly regulated since this early epoch and self-similar evolution are followed in intermediate regions <ref type="bibr">(McDonald et al. 2017, hereafterMD17)</ref>. Recently, <ref type="bibr">Sanders et al. (2018)</ref> reported a self-similarevolution of the thermodynamic properties at all radii for the same large sample but using a different center and a slightly different analysis scheme out to R 500 .</p><p>In this work, we combine deep Chandra and XMM-Newton observations of a sample of the seven highest-redshiftand most massive SPT-selected galaxy clusters beyond a redshift of 1.2 to study the thermodynamic properties of the ICM and their evolution. We take advantage of the sharp point-spread function (PSF) of Chandra to study the small scales (atthis redshift, beyond 1.2,Chandra resolution of 0 5 corresponds to about 5 kpc), while the large effective area of XMM-Newton provides the required photon statistics to measure densities and temperaturesout to large scales. Thus, the combination of Chandra and XMM-Newton allows us to obtain precise and extended density profilesand sufficient photon statistics to measure temperature profiles required to probe the evolution of the ICM properties, e.g., density, temperature,pressure,and entropy, out to the overdensity radius R 500 . The paper is organized as follows: in Section 2 we present the sample properties and the analysis of the XMM-Newton and Chandra data of the sample, in Section 3 we provide our results, the systematic uncertainties are discussed in Section 4, and we summarize our conclusions in Section 5.</p><p>Throughoutthe paper we assume a flat &#923;CDM cosmology with &#937; m &#61600;=&#61600;0.3, &#937; &#923; &#61600;=&#61600;0.7, and H 0 &#61600;=&#61600;70 km s -1 Mpc -1 . All uncertainties quoted correspond to 68% single-parameter confidence intervals unless otherwise stated.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="2.">Cluster Sample and Data Analysis</head></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="2.1.">Cluster Sample</head><p>Our sample consistsof seven SPT-selected high-redshift (z&#61600;&gt;&#61600;1.2) massive clusters of galaxies with signal-to-noise ratio (S/N) greater than 6 and a total SZ-inferred mass greater than 3&#61600;&#215;&#61600;10 14 M e <ref type="bibr">(Bleem et al. 2015)</ref>. The deep XMM-Newton observations of these clusters have been performed in  (PIs E. <ref type="bibr">Bulbul and A. Mantz)</ref>, and Chandra observations were performed in AO-16 through both the XVP program (PIM. McDonald) and two guest observer(GO) programs(PI G. Garmire, S. Murray). The total Chandra and XMM-Newton clean exposure time used in this work is &#8764;2&#61600;Ms (see Table <ref type="table">1</ref>).</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="2.2.">Imaging Analysis</head></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head>2.2.1.XMM-Newton Imaging Analysis</head><p>We follow the data analysis prescription developed by the XMM-Newton Cluster Outskirts Project collaboration (X-COP; <ref type="bibr">Eckert et al. 2017)</ref> with their new background modeling method <ref type="bibr">(Ghirardiniet al. 2018b</ref>).We differ from the X-COP analysis by the fact that we use the mean surface brightness for these high-redshift clusters because it is not really possible to compute the median surfacebrightnessprofile as done in X-COP, since the cells that will be produced will be very few and highly correlated.See Section 4.2 for how these issues influence our results. Thanks to the reduction of the systematic uncertainty on the background below 5% through this method, we are able to measure thermodynamic properties of highredshift clusters outto R 500 . We provide the summary of the analysis below.We use the XMM-Newton Science Analysis System (SAS) and Extended Source Analysis Software (ESAS; <ref type="bibr">Snowden et al. 2008)</ref>, developed to analyze XMM-Newton EPIC observations. In our analysis,we use XMM-SAS v17.0 and CALDB files as of 2019 January (XMM-CCF-REL-362).</p><p>Filtered event files are generated using the XMM-SAS tasks mos-filter and pn-filter.</p><p>The photon countimages are extracted from the filtered event files from three EPIC detectors,MOS1, MOS2, and pn, on board XMM-Newton, in the soft and narrow energy band 0.7-1.2 keV. The choice of this narrow band is to maximize the source-to-background ratio and minimize the systematic uncertainties in the modeling of the EPIC background <ref type="bibr">(Ettori &amp; Molendi 2011)</ref>. To create the total EPIC images,the countimages from the three detectors are summed.Next, we use eexpmap to compute exposure maps by also taking the vignetting effect into account. The exposure maps are also summed using the scaling factors of 1:1:3.44 for MOS1:MOS2:pn detectors, i.e., the ratio between  The high-energy particle background images are generated by using the background images extracted from the unexposed corners of the detectors and rescaling them to the field of view (FOV). After the light-curve cleaning, residual soft protons still contaminate the FOV <ref type="bibr">(Salvetti et al. 2017)</ref>. We measure the soft-proton contamination in the FOV of each observation by calculating the fraction of countrates in the unexposed and exposed portions of the detector in a hard band (7-11.5 keV; <ref type="bibr">Leccardi &amp; Molendi 2008)</ref>. We then generate the 2D soft-proton image <ref type="bibr">(Ghirardini et al. 2018b</ref>, as described in their Appendix A), to model the remaining soft-proton contamination. We constructthe total non-X-ray background (NXB) by summing the high-energy particle background and the residualsoft-proton images.Thus, we obtain total photon images,exposure maps, and total NXB images for each observation.</p><p>To detect and excise point and extended sources in the FOV, we use the XMM-SAS tool ewavelet with a selection of scales in the range of 1-32 pixels with an S/N threshold of 5. We remove all the point sources found by the ewavelet tool from further analysis. We also run CIAO point-source detection tool wavdetect on Chandra images. The sources detected on XMM-Newton and Chandra images are combined to remove missed point sources by ewavelet. See Section 2.2.3 for details on the Chandra analysis.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head>2.2.2.Point-Spread Function Correction for XMM-Newton</head><p>Due to the relatively large size of the PSF of XMM-Newton, some X-ray photons that originate from one particular region on the sky may be detected elsewhere on the detector. XMM-Newton's 5&#8243;-wide PSF (ataim point) needs to be taken into account to correct for this effect due to the small spatial scales of the clusters in our sample <ref type="bibr">(Read et al. 2011)</ref>. To estimate the impact of the PSF on surface brightness profiles, we first create a matrix, PSF i,j , whose value is the fraction of photons originating in the ith annulus in the sky but detected in the jth annulus on the detector.</p><p>In practice,to modelthe PSF and build the PSF matrix, we, following <ref type="bibr">Eckert et al. (2016</ref><ref type="bibr">Eckert et al. ( , 2020))</ref>, build an image of an annulus with a constant value inside the annulus itself and zero outside, with the constant chosen in such a way that the sum of all pixels is 1; this represents the probability density function (pdf) for the true photons generated in the annulus that represents their origin on the plane of the sky. The XMM-Newton mirrors smear this annuluslimited pdf onto a larger fraction of the exposed CCDs. We then use a functional form (e.g., a King profile plus a Gaussian as in <ref type="bibr">Read et al. 2011)</ref> to model the instrumental PSF function in each location of the detector. The observed photons are the result of the convolution of the original sky photons by the PSF function, S b,obs &#61600;=&#61600;PSF&#61600;#&#61600;S b,true . The pdf is no longer limited to the annulus but has spread to the surroundings. The fraction of the pdf, originating from annulus "i," that now is present in annulus "j" is the value that we put in the corresponding line and row of the PSF matrix. An example image of the PSF matrix is given in Figure <ref type="figure">1</ref>. While the majority of the photons that originate from a given annulus are detected in the same region (the largest values are on the diagonal), some fraction of them are detected in a different annulus.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head>2.2.3.Chandra Imaging Analysis</head><p>We process the Chandra observations of the sample using the CIAO&#61600;4.11(Chandra Interactive Analysis of Observations; <ref type="bibr">Fruscione etal. 2006</ref>) and calibration files in CALDB 4.8.2. We filter the data for good time intervals, including the corrections for charge transfer inefficiency <ref type="bibr">(Grant et al. 2005)</ref>. We remove the photons detected in bad CCD columns and hot pixels,compute the calibrated photon energies by applying the ACIS gain maps, and correct for their time dependence. We also remove the time intervals that are affected by the background flares by examining the light curves. We ran wavdetect, the standard CIAO tool to find point sources in Chandra observations, with scales in the range of 1-32 pixels and a threshold for identifying a pixel as belonging to a source of 10 -6 . We merge point sources detected on Chandra images with those detected on XMM-Newton images as described in Section 2.2.1. All point sources detected in this process are excluded from further analysis.</p><p>We extract photon count images in the soft energy band 0.5-2.0 keV, as is routinely done when analyzing Chandra data.</p><p>For the instrumentalbackground we use blank-sky background spectra that are rescaled based on the flux in the hard band 9.5-12 keV to account for variations in the particle background.Exposure maps are generated to correct for the vignetting effect. The particle-background-subtracted, vignetting-corrected images are shown in Figure <ref type="figure">A1</ref>. Due to the small size of Chandra's PSF, 80% of the total encircled counts are detected within 07 from its source. We therefore do not apply any PSF correction to Chandra data.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head>2.2.4.Joint Chandra and XMM-Newton Surface Brightness Analysis</head><p>To compute the surface brightness profile, we first measure the number of photon counts (N c,i ) in concentric annuli around Figure <ref type="figure">1</ref>. Example PSF matrix image used in our analysis. The matrix shows the contribution to the jth annulus from the ith annulus at each position (i, j). The nondiagonaland asymmetric nature ofthe distribution shows that the contribution of the emission from the cluster centerto the outskirts is not negligible and should be corrected for. the cluster center. We find the cluster center by measuring the centroid in a 250-500 kpc aperture on Chandra images following the approach introduced by <ref type="bibr">McDonald et al. (2013)</ref>.This method allows us to find the center of the largescale distribution of the intracluster plasma independent of the core morphology. The widths of the annuliare required to be larger than 2&#8243;, increasing logarithmically,and with at least 30&#61600;counts contained within each annulus. For XMM-Newton the width of these annuli is determined in such a way that each has at least a total of 100 counts and the minimum width is larger than 5&#8243;. We then compute the mean exposure time t i exp, from the exposure map and background counts using the total background map N NXB,i for the two X-ray telescopes. The surface brightnessin each annulus is calculated using the following relation:</p><p>where A reg,i is the area, in arcmin 2 , of each annulus "i."</p><p>From a theoretical point of view, the surface brightness profile is related to the number density through</p><p>where n p and n e are number densities of protons and electrons, respectively, and dl is the integral along the line of sight. We fit the <ref type="bibr">Vikhlinin et al. (2006)</ref> density model to the observed Chandra and XMM-Newton surface brightness data jointly:</p><p>The parameters of the ICM model are constrained by fitting the observed countsN c,i in each annulus against the predicted counts &#956; i (see Equation ( <ref type="formula">5</ref>)) using the following Poisson likelihood:</p><p>( )</p><p>The net number of counts &#956; i inferred by the ICM model in the ith annulus is calculated using the predicted surface brightness, Equation (2),convolved with the PSF matrix, considering the exposed area and time for each annulus, as well as both sky and particle background.</p><p>where t i exp, and A reg,i are the exposure time and area of the annulus "i," respectively, B sky is the cosmic X-ray background (CXB), and N NXB,i are the detector background counts.</p><p>The sum of XMM-Newton and Chandra likelihoods is used as the total likelihood for the fit. We first minimize the c = -&#61682; 2 log 2 using the Nelder-Mead method <ref type="bibr">(Gao &amp; Han 2012)</ref>.Then, we fit using the Bayesian nested sampling algorithm MultiNest <ref type="bibr">(Feroz et al. 2009</ref>) using shallow Gaussian priors centered around the Nelder-Mead method best-fit results and with a standard deviation of 1 (or 2.3 dex)in order to ensure that the fit is not stuck in a local minimum.</p><p>The surface brightness profiles and best-fit models are shown in Figure <ref type="figure">A2</ref>, while the best-fitparameters of the ICM model are given in Table <ref type="table">2</ref>. We note that the emissivity measurements of Chandra and XMM-Newton observatoriesare consistent with each other within 3%; therefore, calibration differences are irrelevant in the measurements of emissivity and number density (as also shown in <ref type="bibr">Bartalucci et al. 2017b</ref>).</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="2.3.">XMM-Newton Spectral Analysis</head><p>We extract spectra using the XMM-ESAS tools mosspectra and pn-spectra <ref type="bibr">(Snowden et al. 2008)</ref>. Redistribution matrices (RMFs) and ancillary response files (ARFs) are created with rmfgen and arfgen,respectively.The point sources (see Section 2.2.1 for details) are excluded from the spectralanalysis.The spectralfitting packageXSPEC v12.10 (Arnaud 1996) with ATOMDB v3.0.9 is used in the analysis <ref type="bibr">(Foster et al. 2012)</ref>. The Galactic column density is allowed to vary within 15% of the measured Leiden/Argentine/Bonn (LAB) Galactic HI survey value in our fits <ref type="bibr">(Kalberla et al. 2005)</ref>. The extended C-statistics are used as an estimator of the goodness of fit <ref type="bibr">(Cash 1979)</ref>. The abundances are normalized to the <ref type="bibr">Asplund et al.(2009)</ref> solar abundance measurements with the mean molecular weight &#956;&#61600;=&#61600;0.5994 and the mean molecular mass per electron &#956; e &#61600;=&#61600;1.1548, and the ratio of the number density of protons to electrons is equal to n p /n e &#61600;=&#61600;0.8527. The MOS spectra are fitted in the energy band of 0.5-12 keV, while we use the 0.5-14 keV energy band for pn. We ignore the energy ranges between 1.2 and 1.9 keV for MOS, and 1.2-1.7 keV and 7.0-9.2 keV for pn due to the presence ofbright and time-variable fluorescence lines. The energy band below 0.5 keV, where the EPIC calibration is uncertain, is eliminated from spectral fits. The source spectrum is modeled with an absorbed single-temperature thermal model apec with varying temperature, metallicity, and normalization. For the clusters with multiple observations the model parameters are tied between multiple spectra and fitted jointly.</p><p>The particle background is determined using the rescaled filter-wheel-closed spectra, which allows us to measure the intensity and the spectral shape. On top of this, we include an additional model component for the residual soft protons <ref type="bibr">(Salvetti et al. 2017)</ref>, modeled as a broken powerlaw with shape fixed (slopes 0.4 and 0.8 and break energy 5 keV; Leccardi&amp; Molendi 2008) and normalization free. Regarding the sky background, we model it as the sum of three components: (i) the CXB with an absorbed power law with photon index fixed to 1.46, (ii) the galactic halo (GH) with an absorbed APEC model with temperature free to vary in the range of 0.1-0.6 keV, and (iii) the Local Bubble (LB) with an APEC model with temperature fixed to 0.11 keV <ref type="bibr">(Leccardi &amp; Molendi 2008;</ref><ref type="bibr">Snowden et al. 2008)</ref>. The normalizations of the CXB, LB, and GH background components are set free. To find the sky parameters, we fit the background region,by extracting a spectrum 5&#8242;(&#8764;5R 500 ) away from the core. We impose Gaussian priors on these parameters with width equal to the parameter uncertainty found in the fitting of the background region.</p><p>We first extract the XMM-Newton spectra within R 500 to measure the redshifts of the clusters from the X-ray data. We fit the spectra within R 500 , so that the statistics are of high quality to determine an accurate X-ray redshift. We fit these spectra using an absorbed single-temperature thermal model with free temperature, metallicity, redshift, and normalization.Taking into account the gain calibration uncertainty of XMM-Newton pn at 3 keV (the redshifted position of the Fe-K line) of 12&#61600;eV (private communication with the XMM-Newton calibration team), we find that the redshifts are consistent with the previously reported photometric <ref type="bibr">(Bleem et al. 2015)</ref> and spectroscopic redshifts <ref type="bibr">(Stalder et al. 2013;</ref><ref type="bibr">Bayliss et al. 2014;</ref><ref type="bibr">Khullar et al. 2019)</ref> within the 2&#963; confidence level for these clusters.A comparison of redshifts based on X-ray data with photometric and spectroscopic redshifts is shown in Figure <ref type="figure">2</ref>. We point out that for SPT-CLJ0459-4947 thepreviously reported redshift (Bocquetet al. 2019) is measured using the position of the Fe-K line from XMM-Newton data from LP by A. Mantz.</p><p>To examine the radial profiles of thermodynamic properties, we next extractthe spectra from concentric annuli with sizes increasing logarithmically around the cluster centroid. The minimum width of annuli is set to be &#8764;15&#8243; to minimize the effect of XMM-Newton's PSF, but still having a large enough statistic to determine the projected temperature. We group the output spectra to ensure having a minimum of 5 counts per bin. The XMM-Newton PSF is taken into account using the crosstalk ARFs generated by the SAS task arfgen. This method allows all the spectra to be cofitted by taking into account the cross-talk contribution to an annulus from another region <ref type="bibr">(Snowden et al. 2008;</ref><ref type="bibr">Ettori et al. 2010)</ref>. The use of flat constant priors on the temperature and metallicity and the use of the "jeffreys" prior on the normalizations (i.e., ( ) = - K K Prior apec apec 1 ) allow us to account for the uncertainty on the sky background, as well as the uncertainty in their free parameters. The spectra are fit using the Markov Chain Monte Carlo (MCMC) implementation in Xspec of the Goodman-Weare algorithm <ref type="bibr">(Goodman &amp; Weare 2010)</ref>, with 50,000 steps and 1000 burn-in period to ensure that we investigate the parameter space and derive the uncertainties on free parameters (temperature,metallicity, and normalization) in our fitting software.At the end of this process,we obtain the best-fit projected temperatures and their covariance matrix, which are easily computed using the MCMC chain.</p><p>To obtain the 3D deprojected temperature profile of each cluster, we project the ICM temperature model on the plane of the sky by taking into account emission weighting to determine spectroscopic-like temperature <ref type="bibr">(Mazzotta et al. 2004)</ref>,</p><p>where &#945;&#61600;=&#61600;3/2, n e is the electron number density, T sl is the predicted 2D spectral temperature, and the temperature model T 3D is a widely used phenomenological model to describe the temperature profiles <ref type="bibr">(Vikhlinin et al. 2006)</ref>.</p><p>The 3D ICM model we used in this work is</p><p>We first minimize the c = -&#61682; 2 log 2 using the Nelder-Mead method <ref type="bibr">(Gao &amp; Han 2012)</ref>. Then, the fit is performed with the MCMC method using the code emcee (Foreman-Mackey et al. 2013) using Gaussian priors centered around the Nelder-Mead method results and with a sigma of 0.5 (or 1.15 dex). We add an additional prior on the temperature fit,by imposing that the pressure derivative decreases monotonically with radius to maintain convective stability.We use 10,000 steps with burn-in length of 5000 steps to have resulting chains independentof the starting position and thinning of 10 to reduce the correlation between consecutive steps. The likelihood adopted in the fit is</p><p>where T obs and T sl are the arrays of the measured spectral temperatures and of the spectroscopic-like projected temperatures as in Equation ( <ref type="formula">6</ref>),respectively,and &#931; i,j is the spectral log-temperature covariance matrix (see Section 2.3). Thus, we use a &#967; 2 -like log-likelihood, where the temperature distribution in each annulus is assumed to be a lognormal (Andreon 2012) and the full covariance between the annuli is considered. The best-fit parametersfor the temperature profile are given in Table <ref type="table">3</ref>. We show an example of the temperature reconstruction process in Figure <ref type="figure">3</ref>.</p><p>Figure <ref type="figure">2</ref>. Comparisonsof X-ray redshifts (in blue) with the photometric redshifts in red <ref type="bibr">(Bleem et al. 2015)</ref> and spectroscopic redshiftsin green <ref type="bibr">(Bayliss et al. 2014;</ref><ref type="bibr">Stalder etal. 2013;</ref><ref type="bibr">Khullar et al. 2019)</ref>. The error bars indicate the sum of statistical and systematic uncertainties at the 1&#963; level.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="3.">Results</head><p>In this section we explore thermodynamic properties (e.g., density, temperature, pressure, and entropy) of the high-redshift SPT clusters in our sample taking advantage of the SPT-SZ survey's clean selection function and its high sensitivity. We further compare the thermodynamic properties of the ICM of the clusters in our sample with the X-COP sample to investigate theirevolution with redshift. The X-COP sample is selected based on the Planck S/N including only lowredshift clusters with z&#61600;&lt;&#61600;0.1 <ref type="bibr">(Ghirardini et al. 2018a, hereafter G18)</ref>. In G18, the authors were able to recover ICM properties out to the virial radius using the jointX-ray and SZ analysis, adding on to the previous studies that probe the region within R 500 by joining X-ray and SZ observations (e.g., <ref type="bibr">Ameglio et al. 2007;</ref><ref type="bibr">Bonamente et al. 2012;</ref><ref type="bibr">Hasler et al. 2012;</ref><ref type="bibr">Eckert et al. 2013a</ref><ref type="bibr">Eckert et al. ,2013b;;</ref><ref type="bibr">Shitanishiet al. 2018)</ref>.We further remark that the analysis done for the high-redshift clusters is almost identical to the analysis applied in G18 for the X-COP cluster sample,allowing for controlled measurement of the evolution in the thermodynamic quantities. We remark that even though in X-COP three clusters have been excluded from the sample because ofdisturbed morphology, the X-COP sample is not biased toward relaxed and cool-core clusters: in fact, only 4 of the 12 X-COP clusters can be considered as relaxed; thus, the fraction of cool cores is very similar to whatis found in SZselected cluster samples (e.g., <ref type="bibr">Rossetti et al. 2017)</ref>.</p><p>The self-similar model <ref type="bibr">(Kaiser 1986</ref>), which assumes purely gravitational collapse, predicts a particular evolution with redshift of the cluster properties once they are scaled based on their common quantities,e.g., mass within an overdensity radius <ref type="bibr">(Voit et al. 2005)</ref>. We therefore measure the mass of our clusters and rescale ourthermodynamic quantities with this mass within R 500 . In the next section we describe our method for the mass reconstruction under the assumption of hydrostatic equilibrium (HE) and then show the thermodynamic profiles and describe their properties.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="3.1.">Total Cluster Mass Reconstruction</head><p>A common way to measure the totalmass M 500 is to use mass proxies calibrated with an X-ray or SZ observable, e.g., L -M or &#958;-M scaling relations (e.g., <ref type="bibr">Pratt et al. 2009;</ref><ref type="bibr">Bocquet et al. 2019;</ref><ref type="bibr">Bulbul et al. 2019)</ref>. However, to avoid introducing a bias into our results by using the evolution in a specific scaling relation, we directly measure the cluster total mass using X-ray observations. The direct measurements based on X-ray data can be obtained from the thermodynamic properties using the assumption of HE and spherical symmetry, i.e.,</p><p>where G is the gravitationalconstant,m p is the mass of the proton, and &#961; g is the gas density. There are several methods that are used in the literature to solve the previous equation (see <ref type="bibr">Ettori et al. 2013</ref>, for a review). Throughoutthis work, we adopt a "forward" modeling approach to obtain a measurement of M 500 , the total cluster mass within R 500 . We make use of the best-fitting density and temperature profiles as recovered in Sections 2.2 and 2.3, respectively,propagating them through the HE equation to recover the mass profile. We point out that we forced the pressure profile to be decreasing at all radii. This method has the advantage of starting from smooth thermodynamic profiles, where the large number of parameters in these functional forms allows us to reproduce the density and temperature profiles over a large radial range.We direct the reader to the Appendix for comparison with literature results and with other mass reconstruction techniques we have employed to solve Equation (9).</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="3.2.">Density,Temperature, Pressure, and Entropy Profiles</head><p>The deprojected electron density profile n e (r) (see Section 2.2.4) obtained from surface brightness analysis is first converted into gas density &#961;(r)&#61600;=&#61600;&#956; e m p n e (r) and then rescaled by the critical density of the universe ( )</p><p>, where H(z)&#61600;=&#61600;H 0 E(z) and ( ) <ref type="figure">4</ref>(a)shows the gas density profiles of the sample. We notice that, in the outskirts, the profiles of the SPT-selected high-z and Planck-selected nearby X-COP clusters are fully consistent with each other, while in the core the SPT-selected high-z profiles are a factor of a few smaller. In the core, the observed scatter (measured as in Equation ( <ref type="formula">6</ref>) in G18) is an orderof magnitude in both the SPT-selected high-z and the Planck-selected nearby X-COP clusters, due to the cool-core/noncool-core states in both samples, i.e., the effect of this dichotomy mostly dominates the scatter near the core. The scatter becomes minimal around 0.4R 500 , at the samelocation where X-COP clusters reach their minima in the scatter.Toward R 500 in the outskirts the scatter increases again in both samples. The increase the high-redshift sample is faster, reaching the value of about 0.35 at R 500 , while the scatter in the X-COP sample remains at 0.2 at the same radius. A comparison of the scatter is seen in Figure <ref type="figure">4(c)</ref>.</p><p>To be able to measure the slope of the density profiles, we perform a piecewise power-law fitting technique as described in detail in G18. Comparing our sample with the nearby X-COP clusters, we find that in the core the slope in our sample is flatter compared to the X-COP clusters, while in the outskirts (&gt;0.3R 500 ) the mean slopes are consistent with each other (see Figure <ref type="figure">4(b)</ref>).</p><p>Next, we study the temperature profiles of the SPT clusters and compare them with the nearby X-COP clusters. For this comparison, the spectroscopic temperature profiles (see Section 2.3 for details) are scaled by the self-similar T 500 , also used in Equation (10) of G18 for the X-COP clusters <ref type="bibr">(Voit et al. 2005)</ref>: where the total mass M 500 is measured in Section 3.1 and used self-consistently in calculations of R 500 and T 500 . In Figure <ref type="figure">5</ref>(a),we compare the rescaled temperature profiles of the SPT high-z clusters with the nearby X-COP clusters. We find that the scaled temperature profiles in the two samples are consistent with each other in the entire radial range out to R 500 . The size of the PSF of XMM-Newton is comparable to the size of the core of these high-redshift clusters; therefore, we cannot resolve well temperatureswithin &lt;150&#61600;kpc,or 0.1R 500 . Performing a piecewise power-law analysis in two radial bins, we obtain similar slopes and the intrinsic scatter in the temperature profiles when comparing them with the X-COP cluster results.</p><p>The pressure profiles are obtained by combining the deprojected density and temperatureprofiles as P&#61600;=&#61600;n e T e . Pressure profiles can be constrained from both X-ray and SZ observations and used for constraining astrophysical properties and the total mass of clusters out to their virial radius <ref type="bibr">(Bonamente et al. 2012;</ref><ref type="bibr">Ghirardini et al. 2018b</ref>).We remark that the pressureand entropy profiles within 0.1R 500 are obtained by combining the resolved density profile with the unresolved temperature profile; hence,results on evolution of pressureprofile within this radius depend heavily on the temperaturemodel adopted. The <ref type="bibr">Vikhlinin et al. (2006)</ref> temperature model is able to reproduce a variety of cluster temperature profiles in the core, and the large uncertainty in the inner partof the profile reflects the large width of the central temperature bin. Therefore, in the relevant figures we warn the readerabout the possible model-dependent sensitivity of our results using gray shadow areas.</p><p>We rescaled the pressure using the self-similar pressure P 500 as described in <ref type="bibr">Nagai et al. (2007)</ref>:</p><p>( ) Figure <ref type="figure">6</ref>(a) shows a comparison of the rescaled pressure profiles of our sample of high-z clusters with the X-COP sample.We find that in the core of SPT high-z clusters the rescaled pressure is on average lower and flatter compared with what is measured in nearby clusters. In the outskirts,pressure becomes consistent with the finding of low-redshift X-COP clusters.The scatter is also fully consistent between high-and low-redshift clusters in all our radial points except the outermost at R 500 , when at high redshift it is 20% higher. Another thermodynamic property thatcould be extracted from X-ray observations is the entropy. Entropy is often used to constrain the clumpiness and self-similarity in cluster outskirts <ref type="bibr">(Urban et al. 2011;</ref><ref type="bibr">Walker et al. 2012;</ref><ref type="bibr">Bulbul et al. 2016</ref>). The entropy profiles are obtained using the relation = - K Tn e 2 3 . Similarly, the entropy is rescaled with the self-similarvalue K 500 for comparison (see <ref type="bibr">Voit et al. 2005)</ref>: In Figure <ref type="figure">7</ref> we show the entropy profiles of the sample, the slope of the entropy, and the intrinsic scatter.An excess is observed in the entropy compared to self-similarity within 0.3R 500 near the core. We attribute this excess to nongravitational processes (e.g., AGN feedback,infalling substructures, merging activities) in the cores. A similar entropy excess in the core was reported in nearby low-redshift clusters <ref type="bibr">(Urban et al. 2014;</ref><ref type="bibr">Bulbul et al. 2016;</ref><ref type="bibr">Ghirardini et al. 2018a;</ref><ref type="bibr">Walker et al. 2019)</ref>, but smaller than the entropy excess observed in these high-z clusters. The high-z entropy excess may be due to the increased incidence ofnongravitationaleffects, e.g., galaxy and clusterformation, and minor mergers athigher redshifts that trigger AGN activity <ref type="bibr">(Hlavacek-Larrondo et al. 2012;</ref><ref type="bibr">McDonald et al.2016;</ref><ref type="bibr">B&#238;rzan et al. 2017)</ref>.</p><p>The entropy profiles are flatin the cores and steepen and become consistent with the self-similar model beyond &#8764;0.2R 500 , similarly to and fully consistent with the entropy profiles in the outskirts of nearby clusters (fora review see <ref type="bibr">Walker et al. 2019, and references therein)</ref>. The intrinsic scatter is comparable for both samples.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="3.3.">Evolution of Thermodynamic Properties with Redshift</head><p>In this section we investigate the redshift evolution of thermodynamic properties of the ICM and measure the deviation from self-similarity of our sample. Following a similar approach described in MD17, we determine the evolution of the density in different radial bins. We characterize the evolution of the thermodynamicquantities using the functions given below:</p><p>where C &#961;,T,P,K representthe deviations with respectto selfsimilar values for the evolution <ref type="bibr">(Kaiser 1986</ref>) of density, temperature, pressure, and entropy.Starting from the density, temperature,pressure,and entropy profiles of the nearby X-COP sample, we infer the expected profiles at the redshifts Figure <ref type="figure">6</ref>. Same as Figure <ref type="figure">4</ref>, but for the pressure profiles. The gray shaded area represents the location within which the temperature profiles are unresolved, wher presented pressure profiles depend on the temperature model adopted.</p><p>of the SPT high-z sample assuming a simple deviation from the self-similar evolution,as indicated in Equation ( <ref type="formula">13</ref>). We then compare these profiles with the thermodynamic profiles of the SPT clusters using a log-likelihood c = -&#61682; log 2 2 to fit and to determine the best-fit evolution parametersC &#961;,T,P,K . The best-fit parametersof these fits are given in Table <ref type="table">4</ref>. The uncertainties of the X-COP profiles, as well as their measured scatter, and the uncertainties on R 500 and Q 500 (see Equations ( <ref type="formula">10</ref>)-( <ref type="formula">12</ref>)) are propagated through the fit. We also include the systematic uncertainties related to our observations in our measurements (see Section 4 for details). The systematic and statistical uncertaintiesare summed in quadratureto estimate the total uncertainty.</p><p>We note that the cluster centers are determined from the Chandra data and initial results are obtained using the centroid of the large-scale ICM emission in this analysis. The choice of cluster center plays an important role especially when measuring the evolution of the central cluster properties <ref type="bibr">(Sanders etal. 2018</ref>). To investigate the effectof the center location,we determine the evolution in density using both the centroid and the X-ray peak. The evolution in density, temperature, pressure, and entropy profiles obtained using both the centroids (red) and X-ray peaks (green) is shown in Figure <ref type="figure">8</ref>. We find no evolution in the density at small radii (&#8764;0.3R 500 ) using large-scale centroids. The self-similar evolution in cluster cores is excluded significantly by &#8764;11&#963;. Using the X-ray peaks instead of the centroids,the evolution values move slightly toward self-similarity in the core. However, the departure from self-similarity is still significant at a &#8764;9&#963; confidence level. We also note thatthe intrinsic scatterin density of high-redshift clusters, shown in Figure <ref type="figure">4</ref>(c), at small radii is similar to that of the low X-COP redshift clusters. Nongravitational phenomena (e.g., AGN feedback, sloshing) dominate the physical processes in cluster cores and can affect the evolution in the core of the clusters. Thus, our finding may suggest that nongravitational physical processesthat regulate cluster cores were already in place since a redshift of 1.8 (with a look-back time of &#8764;10 Gyr). Our results in cluster cores are consistent with the results in MD17 at the 1&#963; confidence level. However, the uncertainties in the measurements are reduced at least by a factor of two. <ref type="bibr">Sanders etal. (2018)</ref> suggestthat use of the X-ray peak instead of centroids could mimic a potential evolution in cluster cores and bias the results in evolution studies. Changing the cluster center does not significantly affect our results.</p><p>At large radii, the evolution in density becomes consistent with the self-similar expectation around 0.1R 500 and remains fully consistent out to R 500 . MD17 also reported the best-fit Figure <ref type="figure">7</ref>. Same as Figure <ref type="figure">4</ref>, but for the entropy profiles. The gray shaded area represents the location within which the temperature profiles are unresolved, where presented entropy profiles depend on the temperature model adopted.</p><p>evolution consistent with the self-similarity; however,due to the limited statistics, the authors could not rule out no evolution scenario. We tightly constrain self-similarity in cluster outskirts and confirm it with a higher significance level. We also observe an increase of the scatter on cluster density profiles (see Figure <ref type="figure">4</ref>) in cluster outskirts. This may imply that although the cluster-to-cluster variance in the outskirts increases because of larger mass accretion rates and merger activity at higher redshifts <ref type="bibr">(Wechsler et al. 2002;</ref><ref type="bibr">Fakhouri &amp; Ma 2009;</ref><ref type="bibr">Tillson et al. 2011;</ref><ref type="bibr">Avestruz et al. 2016)</ref>, the average evolution in density,however,remains consistent with this self-similarity.</p><p>In the case of temperature profiles, we do not measure any significant deviation from self-similarity from the cluster cores out to R 500 . The intrinsic scatter is also consistent with that of the low-redshift clusters within uncertainties (see Figure <ref type="figure">5</ref>). Therefore, the cluster temperature evolution and the cluster-tocluster variance do not seem to change from low to high redshifts.The change of the cluster center makes a very small difference and does not change the results. This is not surprising considering the large uncertainties on temperature measurements.</p><p>We observe a mild evolution in pressure profiles in cluster cores.Similarly, the evolution becomes consistent with selfsimilarity at &#8764;0.1R 500 and larger scales. At small scales, pressure profiles deviate significantly from self-similar evolution at a 6&#963; level.Using the X-ray peak as the cluster center does not change the results significantly.</p><p>Interestingly, in the core, a mildly significant (&#8764;3&#963; confidence)evolution is observed forthe entropy,if we use the centroid as the cluster center. Changing the cluster center to the X-ray peak reduces significantly the observed evolution. In the outskirts the evolution becomes fully consistent with selfsimilarity, regardless of the center used.</p><p>It is important to remind the reader that the evolution measured in clustercores for pressure and entropy is quite dependent on the adopted cluster temperature model, because the first temperature bin is very large, encapsulating the entire cluster core, &#61600;&lt;0.1R 500 . Note. In each single table the first two columns represent the inner and outer radial ranges in which we have looked for the evolution. The third column represents t measured evolution with redshifts, along with its statistical and systematic uncertainty. The fourth column represents the significance measured in number of sigma the difference between the measured evolution and the evolution predicted by the self-similar expectation.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="3.4.">Polytropic Index</head><p>The global structure of the ICM can be effectively described by a polytropic equation ofstate P e &#61600;=&#61600;K&#961; &#915; , where the polytropic index is indicative of stratification of the ICM <ref type="bibr">(Shaw et al. 2010)</ref>. Both simulations <ref type="bibr">(Komatsu &amp; Seljak 2001;</ref><ref type="bibr">Ostriker et al. 2005;</ref><ref type="bibr">Ascasibaret al. 2006;</ref><ref type="bibr">Capelo et al. 2012)</ref> and observations <ref type="bibr">(Markevitch et al. 1998;</ref><ref type="bibr">Sanderson et al. 2003;</ref><ref type="bibr">Bulbul et al. 2010;</ref><ref type="bibr">Eckert et al. 2015;</ref><ref type="bibr">Ghirardini et al. 2019)</ref> find that the stratification of the ICM, especially in the outer part, is well represented by a polytropic equation of state with &#915; in the range of 1.1-1.3. In particular,the X-COP collaboration reports that the value of &#915; in cluster outskirts,where &#961;/&#961; c &#61600;&#61600;&#61600;400, is &#915;&#61600;=&#61600;1.17&#61600;&#177;&#61600;0.01 at redshifts below 0.1. However,the polytropic index in the highredshift universe, or its evolution, has never been investigated. We find that the polytropic index (see also Figure <ref type="figure">9</ref>), is 1.19&#61600;&#177;&#61600;0.05 in low-density regions, i.e., in the cluster outskirts. This value is fully consistent with the value measured at low redshifts in the X-COP clusters,indicating that there is no significant evolution with redshift,i.e., the ICM stratification is the same at low and high redshift. In high-density regions, i.e., in the core, we are not able to resolve the index owing to the large size of the XMM-Newton PSF.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="4.">Systematics</head><p>In this section we examine several systematic uncertainties that affect our results on the evolution of the thermodynamic properties of clusters, evaluating their magnitudes. The variation of the thermodynamic property Q can be converted into the systematic uncertainty on the evolution following the formula below:</p><p>where &#7825; is the average redshiftof our sample and k is the systematic uncertainty on the evolution of each thermodynamic property Q. Solving this equation for the systematic uncertainty k gives the following equation:</p><p>(</p><p>Evolution in the density profiles as a function of redshift obtained using the centroid (in red) and X-ray peak (in blue). The red shaded region around ou data points represents the sum in quadrature of the statistical and systematic uncertainties (see Section 4 for details). The yellow shaded area represents the same as found by MD17. Zero values of 2&#61600;+&#61600;C &#961; indicate no evolution with redshift. The self-similar evolution of C &#961; &#61600;=&#61600;0 (corresponding to ( ) r &#181; E z 2 ) is represented by a horizontal dashed line. The other panels are the same but for (b) temperature with self-similar predicted evolution corresponding to ( ) &#181; T E z 2 3 , (c) pressure with selfsimilar predicted evolution corresponding to ( ) &#181; P E z 2 3 , and (d) entropy with self-similar predicted evolution corresponding to ( ) &#181;</p><p>-K E z 2 3 . Moreover, for pressure and entropy, below 0.1R 500 the values of the evolution are extrapolated because temperature measurements are not resolved on smaller scales. The vertical dashed represents the location of R 500 in all panels. Moreover, in the panels where entropy and pressure are presented, we mark with gray shaded areas the core region, whe the temperature profiles are unresolved.</p><p>We consider the systematic uncertainties related to hydrostatic mass bias, clumping factor, cluster progenitors, and calibration differencesbetween XMM-Newton and Chandra below. In Table <ref type="table">5</ref> we show the amplitude of the mass bias on each thermodynamic quantity in the core and at R 500 .</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="4.1.">Mass Bias</head><p>The thermodynamic profiles and their evolution depend on the mass thatis used to rescale the profiles. However,given that the low-redshift X-COP sample and high-redshiftSPT sample have very similar selection criteria,i.e., a selection based on SZ S/N,and the masses are obtained in both cases assuming HE, the mass rescaling is expected to affect the two samples by the same amount; hence, the evolution should not be affected. In this section we search for possible systematics in the hydrostatic mass measurement that affect differently the low-and high-redshift clusters.</p><p>An estimate of this mass systematic bias can be obtained by measuring the average ratio between several massmeasurements. In Section 3.1 we have described our reference method of solving the HE equation to measure M 500 , and in the Appendix we compare this measure with other techniques and other masses in the literature obtained from scaling relations <ref type="bibr">(Bleem et al. 2015;</ref><ref type="bibr">Bocquet et al. 2019)</ref>. Figure <ref type="figure">10</ref> shows the cluster masses obtained through thesemethods.To estimate the bias, we calculate the average ratio between the different mass available and the mass obtained using the reference method, described in Section 3.1. Since the error bars are not homogeneous, we apply a bootstrap method, i.e., we measure the mass bias 10 6 times, where each time a new distribution of masses is drawn from the masses shown in Figure <ref type="figure">10</ref>. This method yields a mass bias of 1&#61600;+&#61600;b&#61600;=&#61600;1.12&#61600;&#177;&#61600;0.01. The result implies that high-redshift clusters have potentially 12% higher hydrostatic masses compared to the nearby clusters. Given that clusters athigh redshifts are still forming and not yet thermalized, an increase in the nonthermalpressure support due to gas motions in their outskirts and elevated AGN activity in their cores, resulting in an increase in mass bias with redshift, is expected.</p><p>If the hydrostatic masses we use in this work are biased (with respect to the low-redshift masses) by a factor of (1&#61600;+&#61600;b), this bia translates to a bias in the fiducial radius that can be written as</p><p>And it translates into an uncertainty on a rescaled thermodynamic property Q as</p><p>where Q&#61600;=&#61600;T,&#61600;P,&#61600;K.</p><p>Figure <ref type="figure">9</ref>. Rescaled temperature against rescaled density in high-redshift cluster sample (in red) and in low-redshift clusters (in black; <ref type="bibr">Ghirardini et al. 2019</ref>). The lines represent the best-fit broken power law to the data. In particular, we find that the slope in the relation is consistent in low-and high-redshift clusters in the low-density regime, i.e., in cluster outskirts,supporting again the selfsimilar model of cluster evolution. Note. The thermodynamic biases are, from left to right, (1) hydrostatic bias caused by how the profiles are rescaled, (2) clumping bias caused by the presence of unresolved clumps, (3) bias caused by the fact that SPT high-z clusters are not exactly the progenitors of the redshift 0 clusters we are comparing them with, and ( calibration bias caused by difference between Chandra and XMM-Newton temperatures.</p><p>Figure <ref type="figure">10</ref>. Mass comparison for the object in our sample. In red are the masses from the SPT catalog <ref type="bibr">(Bleem et al. 2015)</ref> and the massesfrom the SPT cosmological analysis <ref type="bibr">(Bocquet et al. 2019)</ref>.In black are the masses recovered by MD17 using the M gas -M tot scaling relation. In green are the NFW best-fit masses in the two cases described in the text. In blue are the forward best-fit M 500 computed using a functional form to fit the temperature and density profiles.</p><p>Using the mass bias obtained above, we then estimate the corresponding systematicbias in the evolution. This bias affects both x-and y-axes, except for density, where the rescaling on the y-axis is independent of mass.The bias is translated into</p><p>on the x-axis and</p><p>on the y-axis; then, by summing up in quadrature these two values and applying Equations ( <ref type="formula">14</ref>) and ( <ref type="formula">15</ref>), we measure the systematic uncertainty on the evolution of the thermodynamic quantities caused by the mass bias.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="4.2.">Clumping Factor</head><p>Unresolved clumps in ourobservations can lead to higher local densities measured and can bias the observed thermodynamic quantities. In G18 the authors correct the density for the presence of these clumps by both removing the extended sources contaminating the FOV and computing the median of surface brightness distribution in each annulus, which has been shown to be unbiased by the presenceof high-density unresolved substructures <ref type="bibr">(Roncarelli et al. 2013;</ref><ref type="bibr">Zhuravleva et al. 2013)</ref>. In particular, to compute the median, a Voronoi tessellation algorithm needs to be performed <ref type="bibr">(Diehl &amp; Statler 2006)</ref> to produce cells containing surface brightness elements.In this work, we eliminate the detected pointand extended sources from our analysis. Due to low counts observed and the small extension of the clusters on the sky, the cells produced via the Voronoi tessellation algorithm would be very few and highly correlated with each other. Therefore, it is not possible to compute the median of the surface brightness distribution in the same way as applied to the X-COP sample. Instead,we estimate this bias by adopting the upper limit of 10% within R 500 measured in a sample of ROSAT clusters in <ref type="bibr">Eckert et al. (2015)</ref>. We find that the density profiles are biased by a systematic uncertainty of&#916;&#961;/&#961;&#61600;=&#61600;0.10. This translates into a systematic on the density measurements of &#8764;0.06 (see Equation ( <ref type="formula">15</ref>)).</p><p>For the other thermodynamic properties, we combine the effect aforementioned with the bias of 5% in the pressure arising by the presence of clumps (as measured in simulations by <ref type="bibr">Khedekar et al. 2013</ref>, where the 5% refers to the upper limit within R 500 ). This translatesinto a bias of 5% on the temperature, consistentwith the predicted theoretical bias by <ref type="bibr">Avestruz et al.(2016)</ref>,and a -2% bias on the entropy.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="4.3.">Progenitors</head><p>It is possible thatthese SPT-selected high-redshift clusters are notthe progenitors of the low-redshift clusters in X-COP. In fact, the predicted mass of the SPT clusters is expected to be greaterthan 10 15 M e at redshift zero when the mass growth curve is taken into account <ref type="bibr">(Fakhouriet al. 2010)</ref>. Therefore, the SPT-selected clusters are more massive than the X-COP clusters <ref type="bibr">(Ettori et al. 2018)</ref>, where the reported masses are less than 10 15 M e . We treat the effect due to the fact that the X-COP clusters could be evolved from a different population of clusters than the SPT clusters as a systematic uncertainty.</p><p>To estimate this bias, we assume that the gas density follows the dark matter density as a first approximation. We then use a concentration-mass-redshift relation in <ref type="bibr">Amodeo etal. (2016)</ref> to calculate the relative change in the density from a cluster with a mass of 15&#61600;&#215;&#61600;10 14 M e , i.e., the expected mass of SPT clusters ata redshiftof 0 <ref type="bibr">(Fakhouri et al. 2010)</ref>,to a mass of 7&#61600;&#215;&#61600;10 14 M e , i.e., the average mass of X-COP clusters. Assuming thatpressure follows the universal pressure profile, we estimate the thermodynamic quantities. These values are then propagated as systematic errors as shown in Equation ( <ref type="formula">15</ref>). The results in the core and in the outskirts are given in Table <ref type="table">5</ref>. We note that the self-similar model predicts an evolution that is independentof mass. Therefore, once the thermodynamic quantities are rescaled with their self-similar value, the fact that they are too massive to be the progenitors of the X-COP clusters is of minor importance, especially at large radii.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="4.4.">Calibration Difference between Chandra and XMM-Newton</head><p>Calibration differences between Chandra and XMM-Newton are described in the literature. Temperature measurements can be biased up to 40% depending on the energy band used and cluster temperature(e.g., <ref type="bibr">Schellenberger et al. 2015, and references therein)</ref>.On the other hand, density measurements by Chandra and XMM-Newton are fully consistent within the uncertainties (see <ref type="bibr">Bartalucci et al. 2017a</ref>; see also Section 2.2.4).</p><p>To quantify the bias due to calibration differences, we extract both XMM-Newton and Chandra spectra of the region within R 500 and fit the spectra using a single-temperature thermal apec model. We note that in the case of the SPT high-redshift clusters it is not possible to measure temperature profiles using Chandra observations in several radial bins owing to the limited statistics. A comparison of measured single temperatures is shown in Figure <ref type="figure">11</ref>. We find that the temperaturemeasurements are consistent with each other within statistical uncertainties. However, we note that the uncertainties on the Chandra measurements are large because of limited statistics.</p><p>To estimate the systematic uncertainty on each thermodynamic quantity Q caused by this discrepancy in the temperature is not trivial. The increase of the temperature would lead to an increase in the total massby the same amount, if the slope of the temperature profile does not change. <ref type="bibr">Schellenberger et al. (2015)</ref> report that temperature measurements based on Chandra data are, on average, 17% higher than those derived from XMM-Newton for the average mass of the clusters in our SPT sample. Thus,a systematic of 17% on the temperature becomes a 17% systematic on the mass, and thus a 5.7% bias on R 500 (one-third considering the propagation of uncertainty) and 11.3% bias on Q 500 (twothirds considering that all self-similar quantities depend on mass with power of 2/3). Thus, the variation on each thermodynamic quantity is We point out that, for the last two terms,the variation on the rescaled thermodynamic quantity from the radial and the Q 500 rescaling is in the opposite direction with respect to the systematic bias on the quantity Q.Thus,by computing the slope at each radius,we get the relative rescaled thermodynamic variation at each radii, and finally, using Equation ( <ref type="formula">15</ref>), we obtain the systematic bias affecting the evolution of each thermodynamic quantity as given in Table <ref type="table">5</ref>.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head n="5.">Conclusions</head><p>In this paper we have studied the thermodynamic profiles for the seven mostmassive clusters at redshift above 1.2 in the SPT-SZ survey. These clusters were observed by both Chandra and XMM-Newton for a totalclean exposure time of about 2 Ms. We combine the data from these two telescopes to recover density, temperature,pressure, and entropy profiles and examine their evolution with redshift from cluster cores to outskirts. Furthermore, we measure the temperature profiles of a complete setof SPT-selected high-redshift clusters for the first time, allowing us to reconstructthe total cluster masses under the assumption of HE. Our results include the systematic uncertainties that are extensively studied in Section 4.</p><p>Deep XMM-Newton observations of the SPT-selected clusters have sufficient statistics to determine the redshifts from the X-ray data alone. The Fe-K line at 6.7 keV (rest frame) is clearly detected in the spectrum of each cluster in the sample. The centroids of these emission lines are used to measure the redshifts. We show that the redshifts obtained from the X-ray data of the SPT high-z clusters are consistent with the previously reported redshifts obtained through optical photometry and spectroscopy <ref type="bibr">(Bayliss et al. 2014;</ref><ref type="bibr">Bleem et al. 2015;</ref><ref type="bibr">Khullar et al. 2019)</ref>.</p><p>Combination of Chandra's high spatial resolution and XMM-Newton's large FOV and effective area is the most powerful way to measure thermodynamic profiles of clusters at high redshifts, z&#61600;&gt;&#61600;1.2 from their cores (0.01R 500 ) to the outskirts (R 500 ). Accurate measurements of temperature profiles enable a few key measurements for these clusters, e.g., total mass, pressure, and entropy. We are able to measure their total mass through the HE assumption with relatively small uncertainties(10%-20%) at these redshifts.The hydrostatic masses are generally in good agreement with reported masses in the literature obtained through SZ S/N and scaling relations <ref type="bibr">(Bleem et al.2015)</ref>.</p><p>We further measure the density, temperature, pressure, and entropy profiles of the high-z SPT cluster sample and compare their distributions with the previously reported thermodynamic properties of the nearby X-COP clusters. The scatters of all the thermodynamic quantities are similar in low-and high-redshift clusters in small spatial scales near the cores. At large radii, the scatterincreases more steeply in the sample of high-redshift clusters. This may be due to an increased frequency of merger events and higher mass accretion rate at high redshifts <ref type="bibr">(Wechsleret al. 2002;</ref><ref type="bibr">Fakhouri &amp; Ma 2009;</ref><ref type="bibr">Tillson et al. 2011)</ref>.</p><p>The average profiles of density, temperature, pressure, and entropy of high-z clusters are self-similar and consistent with those of the X-COP clusters at large spatialscales near R 500 . Temperature profiles of high-redshift clusters are self-similar at all radii. We also report that the polytropic index (1.19 &#177; 0.05) is fully consistentwith that measured at low-redshiftclusters, indicating that there is no significantevolution with redshift. The high observed scatter in density, pressure, and entropy in cluster cores is due to the cool-core/non-cool-core dichotomy in these cluster samples.The scatter in the thermodynamic properties becomes minimal at 0.4R 500 and increases toward R 500 in the SPT-selected high-z clusters. The increase in the mergerfrequency and mass accretion rate in high-z clusters may contribute to high scatterin cluster outskirts <ref type="bibr">(Wechsler et al. 2002;</ref><ref type="bibr">Fakhouri &amp; Ma 2009;</ref><ref type="bibr">Tillson et al. 2011)</ref>.</p><p>We are also able to constrain the evolution in density and temperatureprofiles of the cluster. Measurementsof the evolution in entropy and pressure profiles with redshift also become available owing to precise temperature constraints for the first time. We find that the evolution in thermodynamic profiles deviates significantly from the self-similar evolution in cluster cores, while in the outskirts the profiles are on average in agreementwith the prediction from the self-similar model. We find no evidence for evolution in the density in cluster cores,confirming the results in MD17.We point out that the analysis performed in this paperand the one in MD17 are different in how self-similarity has been probed. We have considered two high-S/N cluster samplesat low and high redshift, while in MD17 the authors have considered &#8764;100 low-S/N clusters. Therefore,it is striking that two different analyses on two different samples yield the same results on the evolution of cluster density profiles.We observe only mild evolution in pressure and entropy profiles in cluster cores. On the other hand, the evolution of temperature profilesis in agreement with self-similarity. Utilization of the X-ray peak instead of the centroid of the large-scale emission does not significantly affect our results (it changes the measured evolution in the core toward self-similarity but does not change significantly the significance).</p><p>Planned and future X-ray telescopes with sufficiently small spatial resolution and large effective area (e.g., Athena, Lynx) will provide sufficientstatistics to precisely measure temperature and density profiles down to kiloparsec scales in the cores of a large sample of clusters <ref type="bibr">(Nandra et al. 2013;</ref><ref type="bibr">Gaskin et al. 2019)</ref>. These measurements will allow us to probe in detail the role of AGN feedback in the first clusters formed and to shed light on the accretion processes in cluster outskirts and the structure formation in the universe. Since the large bin size of the annuli is caused by the large XMM-Newton PSF of about 15&#8243;, which corresponds to a physical size of 150 kpc, the constraints on the concentration parameter are very weak, meaning that the concentration is almost unconstrained. The fit is done using the code emcee (Foreman-Mackey et al. 2013), starting from a standard maximum likelihood fit, &#967; 2 minimization using the Nelder-Mead method (Gao &amp; Han 2012), using 10,000 steps with burning length of5000 steps to have resulting chains independent from the starting position,and thinning of 10 in order to reduce the correlation between consecutive steps.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head>&#180;+</head></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head>A.3. Reconstructed Mass</head><p>Our reconstructed M 500 , using the method described above and in Section 3.1,are shown in Figure <ref type="figure">10</ref>,and displayed in Table <ref type="table">6</ref>. We compare ourmass reconstruction among themselves,and with the SPT masses as calculated in the catalog <ref type="bibr">(Bleem et al. 2015)</ref> using the M-&#950; fixed scaling relation, with the masses calculated from the scaling relations obtained for the SPT cosmological results <ref type="bibr">(Bocquet et al. 2019)</ref>, and with the masses used in MD17, which come from the M gas -M tot scaling relations <ref type="bibr">(Vikhlinin et al. 2009)</ref>.Overall the masses we measureare consistentwith all the other masseswe are comparing with, with two peculiar cases:(1) SPT-CLJ0459-4947, for which the masses coming from the forward reconstruction agree with the othermassesin the literature, i.e., the two SPT masses and the masses in MD17, but the NFW reconstructionprefers a higher mass. This can potentially indicate that the NFW mass model could not be the best model to describe the dark matter potentialfor this object. (2) SPT-CLJ2341-5724,which has all the massescoming from our analysis consistent within 1&#963;; however,when comparing with the literature masses, we find that these are much higher than what we measure, indicating the possibility that this cluster does not fall on the scaling relations used to determine the literature masses. The recovered mass of SPT-CLJ0205-5829 has very large uncertainties. This is because the XMM-Newton 55 ks observation 0803050201 is highly flared, with only about 10 ks remaining after flare removal, and on top of that this cluster has a point source very close to the cluster center, thus decreasing the photon statistics, with the resulting effect being larger error bars for the temperature, translating into large error bars on the mass since M&#61600;&#8764;&#61600;T.   <ref type="formula">2006</ref>) functional (solid black line) form plus a constant sky background (horizontal dotted line); it is convolved with the instrumental PSF and is shown with a blue line. In the case of Chandra the PSF is simply a diagonal matrix with ones on the diagonal, while for XMM-Newton it is calculated as in Section 2.2.2. In the bottom panels we show the residuals (</p><p>). The dashed vertical line represents the location of R 500 , as measured by solving the HE equation (Equation ( <ref type="formula">9</ref>)) using the "forward T" method (see Section 3.1).  </p></div><note xmlns="http://www.tei-c.org/ns/1.0" place="foot" xml:id="foot_0"><p>The Astrophysical Journal, 910:14 (24pp),2021 March 20 Ghirardini et al.</p></note>
			<note xmlns="http://www.tei-c.org/ns/1.0" place="foot" xml:id="foot_1"><p>The Astrophysical Journal, 910:14 (24pp),2021 March 20Ghirardini et al.   Figure A1. (Continued.)   </p></note>
			<note xmlns="http://www.tei-c.org/ns/1.0" place="foot" xml:id="foot_2"><p>The Astrophysical Journal, 910:14 (24pp),2021 March 20Ghirardini et al.   Figure A2. (Continued.)   </p></note>
		</body>
		</text>
</TEI>
