,

Context. The thermodynamical properties of the intracluster medium (ICM) are driven by scale-free gravitational collapse, but they also reflect the rich astrophysical processes at play in galaxy clusters. At low masses ( ∼ 10 14 M ⊙ ) and high redshift ( z ≳ 1), these properties remain poorly constrained, observationally speaking, due to the di ffi culty in obtaining resolved and sensitive data. Aims. We aim to investigate the inner structure of the ICM as seen through the Sunyaev-Zel’dovich (SZ) e ff ect in this regime of mass and redshift. We focused on the thermal pressure profile and the scaling relation between SZ flux and mass, namely the Y SZ − M scaling relation. Methods. The three galaxy clusters XLSSC 072 ( z = 1 . 002), XLSSC 100 ( z = 0 . 915), and XLSSC 102 ( z = 0 . 969), with M 500 ∼ 2 × 10 14 M ⊙ , were selected from the XXL X-ray survey and observed with the NIKA2 millimeter camera to image their SZ signal. XMM-Newton X-ray data were used as a complement to the NIKA2 data to derive masses based on the Y X − M relation and the hydrostatic equilibrium. Results. The SZ images of the three clusters, along with the X-ray and optical data, indicate dynamical activity related to merging events. The pressure profile is consistent with that expected for morphologically disturbed systems, with a relatively flat core and a shallow outer slope. Despite significant disturbances in the ICM, the three high-redshift low-mass clusters follow the Y SZ − M relation expected from standard evolution remarkably well. Conclusions. These results indicate that the dominant physics that drives cluster evolution is already in place by z ∼ 1, at least for systems with masses above M 500 ∼ 10 14 M ⊙ .


Introduction
The presence of a diffuse hot gas component that permeates galaxy clusters, the intracluster medium (ICM), was revealed by X-ray observations in the 1970s (see Sarazin 1986 andBiviano 2000 for historical reviews).Thanks to subsequent observational and theoretical achievements, it is now established that clusters are dominated by dark matter (∼ 80%), that the ICM accounts for most of their baryonic matter (∼ 12%), and that galaxies only provide the remaining few percent of the total cluster mass.Clusters form at the intersection of cosmic filaments and trace the growth of cosmic structures, as peaks in the matter density field.As such, they are recognized as unique astrophysical laboratories and as important cosmological probes (e.g., Vikhlinin Send offprint requests to: Rémi Adam (remi.adam@oca.eu) ⋆ The FITS file of the published NIKA2 maps are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr(XXXXX) or via XXXXX ⋆⋆ Deceased et al. 2009b;Planck Collaboration et al. 2014;Bocquet et al. 2019).We invite the reader to consult Allen et al. (2011) for a review.
The formation of galaxy clusters and their ICM is mainly driven by the gravitational collapse of large-scale structures (Kravtsov & Borgani 2012).This implies that, to first order, clusters are self-similar objects (Kaiser 1986) whose observational properties follow well-predicted behaviors once rescaled in mass and redshift.However, clusters are also affected by complex astrophysical processes related to gas dynamics and the formation of galaxies.The ICM is believed to be established in the early phase of cluster formation history and to be continuously fed by the merging of smaller groups and the accretion of the surrounding material; the infalling gas kinetic energy is mostly converted into heat via shocks and turbulent cascades.In parallel, a fraction of the baryons condensates into stars, eventually inducing significant supernovae or active galactic nucleus (AGN) feedback onto the ICM (Fabian 2012;Hlavacek-Larrondo et al. 2022).These processes should leave an imprint both in the inner structure of the ICM and in the scaling relation between global integrated properties of the clusters (Lovisari & Maughan 2022).Characterizing the ICM thermodynamics is therefore an excellent way to address the nature of the astrophysical processes associated with cluster formation; it is a necessary step when using clusters to probe the growth of large-scale structures (Voit 2005).
The thermal Sunyaev-Zel'dovich (SZ) effect (Sunyaev & Zeldovich 1970, 1972) provides an independent and complementary probe of the ICM to X-ray observations (Birkinshaw 1999;Mroczkowski et al. 2019, for reviews).It is due to the inverse Compton scattering of cosmic microwave background (CMB) photons on ICM electrons.Its surface brightness is independent of redshift, which makes it particularly competitive for distant objects provided that sufficiently high angular resolution and sensitivity are available (e.g., Korngut et al. 2011, Kitayama et al. 2016, and Adam et al. 2016 for the use of relevant facilities).Unlike X-ray observations, which rely on the combination of gas density and temperature (measured using spectroscopy) to infer the pressure, the SZ effect directly measures the thermal ICM electron pressure.The pressure profile is an excellent tracer of the matter distribution because it reflects how the gas is compressed in the potential well.Similarly, the SZ flux, Y SZ , tracks the total mass with a low intrinsic scatter since it measures the thermal energy directly, itself related to the depth of the potential well (Nagai et al. 2007;Pratt et al. 2019).Therefore, the SZ effect is also an excellent way to detect clusters, with a clean selection function, and it has been given much attention not only for studying the SZ signal directly, but also in follow-up observations of SZ-selected samples (e.g., Sanders et al. 2018; Bartalucci et al. 2019).For all these reasons, the Y SZ − M relation and the pressure profile are key tools that require detailed characterization if one wants to fully benefit from the statistical power of SZ surveys for cluster cosmology and astrophysics.
The thermal pressure profile and the Y SZ − M scaling relation have been deeply investigated and calibrated up to intermediate redshifts (e.g., Arnaud et al. 2010;Planck Collaboration et al. 2013;Ghirardini et al. 2019;Bonamente et al. 2008;Planck Collaboration et al. 2014;Medezinski et al. 2018).However, nontrivial redshift evolution is expected (Le Brun et al. 2017) because of the changes in the cluster mass-accretion rate and the merger activity (Fakhouri et al. 2010;Fakhouri & Ma 2010) or the enhanced star formation and AGN activity (Alberts et al. 2016).Yet, current attempts, either in SZ or in follow-up X-ray observations of SZ-selected samples, did not report significant nonstandard evolution of the bulk ICM properties for massive systems out to z ∼ 2 (see McDonald et al. 2017, Bartalucci et al. 2017, Mantz et al. 2018 (hereafter XXL Paper XVII), and Ghirardini et al. 2021).At lower masses, however, clusters are expected to be more affected by the gas dynamics and the interaction with galaxies due to their shallower potential well, so supplementary deviations in the ICM properties are expected in this regime (M ≲ 10 14 M ⊙ , e.g., Pop et al. 2022).Although low-mass clusters at high redshifts will represent a large fraction of the detections of ongoing and future cluster surveys (e.g., eROSITA, Euclid, LSST, CMB-Stage4, see Pillepich et al. 2012, Bulbul et al. 2022, Euclid Collaboration et al. 2019, LSST Science Collaboration et al. 2009, and Abazajian et al. 2016), their detailed SZ observational properties remain nearly unexplored to date due to the difficulty in obtaining sufficiently high-quality data (see Dicker et al. 2020, for MUSTANG2 observations).Given their expected high sensitivity to cluster astrophysics, dedicated SZ follow-ups with resolved observations are thus becoming essential in efforts to further advance cluster cosmology and astrophysics.
In this paper, we report on SZ observations of three low-mass (M 500 ∼ 2 × 10 14 M ⊙1 ) high-redshift (z ∼ 1) galaxy clusters selected from the XXL X-ray survey (Pierre et al. 2016, hereafter XXL Paper I): XLSSC 072 (z = 1.002),XLSSC 100 (z = 0.915), and XLSSC 102 (z = 0.969).They were imaged with the New IRAM KIDs Array 2 (NIKA2) millimeter camera (Adam et al. 2018a) from the Institut de Radio Astronomie Millimétrique (IRAM) 30-meter telescope at 150 GHz and 260 GHz.Given the nature of these objects in terms of mass and redshift, we aim to test the standard evolution of the ICM as calibrated on nearby massive systems in a regime that has not been explored with resolved SZ data yet and where astrophysical processes should be more effective.The data, complemented with X-ray and optical observations are used to investigate the dynamical state of the clusters.The SZ images are used to derive the thermal pressure profiles and extract the SZ fluxes, which are compared with standard evolution expectations once the masses are extracted under different assumptions using SZ and/or X-ray data.This work extends over the earlier multiwavelength analysis of XLSSC 102 (Ricci et al. 2020), hereafter XXL Paper XLIV2 , by adding new NIKA2 observations for two clusters, the use of new X-ray data for XLSSC 102, and by extending the analysis methodology to recover the cluster pressure profile and their masses.
The paper is organized as follows.In Section 2, we present the different datasets used in this work.In Section 3, we describe the multiwavelength morphological comparison and discuss the dynamical state of the clusters.We present the modeling and the data analysis methodology in Section 4. The results are reviewed in Section 5, and the main conclusions are summarized in Section 6.Throughout the paper, we assume a flat ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 and Ω m = 0.3.The Hubble parameter at redshift z normalized to the present-day value is defined as E(z) = H(z)/H 0 .

Data
In this section, we discuss the selection of the clusters and present the main data that were used.

Cluster sample
The target clusters were selected from the XXL survey (see also Ricci 2018 for details), performed with XMM-Newton in the Xrays (Paper I).Thanks to its selection function, the XXL survey allowed us to identify clusters at low masses and high redshift (Pacaud et al. 2016, hereafter XXL Paper II).We only considered the most securely detected XXL clusters (classified as C1) from the northern part of the XXL survey (XXL-N) that are observable from the IRAM 30m telescope.We requested detections in the optical using the galaxy overdensity to make sure that the clusters were also confirmed independently from the ICM content (see Section 2.4 for details).Only detections with robust    spectroscopic redshift estimates were accounted for (Paper XX).
According to these criteria, we selected the three clusters at redshift z ∼ 1 with X-ray data of sufficient quality to allow a reliable combination with NIKA2, both in terms of images and radial profiles.Although the number of objects is limited, we expect them to be representative of other systems with similar parameters.All targets are part of the 100 brightest XXL clusters (Paper II): XLSSC 072, XLSSC 100, and XLSSC 102.Different mass estimates were obtained from the XXL survey.In Paper XX, count rates together with scaling relations were used iteratively to infer the mass and the temperature without relying on X-ray spectroscopy.Spectroscopic temperatures were extracted within 300 kpc of the cluster center from XMM-Newton (Giles et al. 2016;Adami et al. 2018, hereafter XXL Paper III and XXL Paper XX).The masses were then estimated according to the mass-temperature relation calibrated using weak lensing data (Lieu et al. 2016, XXL Paper IV).A similar approach was also performed by Umetsu et al. (2020), using Paper XX temperatures, where more reliable weak lensing masses were measured thanks to Hyper Suprime-Cam data (Aihara et al. 2022) and to the use of a less restrictive prior on the concentration-mass relation.
The three selected clusters are located in the footprint surveyed by the Atacama Cosmology Telescope (ACT, Hilton et al. 2018Hilton et al. , 2021)).The ACT cluster catalog reports detections for which the signal-to-noise ratio (S/N) is larger than four.The masses were obtained by matching the normalization of the universal pressure profile (UPP) as calibrated by Arnaud et al. (2010) to the ACT data.Masses rescaled using a richnessbased weak-lensing mass calibration factors are also provided.XLSSC 102 is reported in both the Hilton et al. (2018) and Hilton et al. (2021) catalogs with a S/N of 8.3 and 12.1, respectively.XLSSC 100 is not reported in either catalog.XLSSC 072 is only reported in Hilton et al. (2018), with a S/N of 6.5.
In the case of XLSSC 072, a dedicated analysis was performed by Duffy et al. (2022), hereafter XXL Paper XLVIII, using a deep XMM-Newton follow-up of the target, and thus clearly of higher quality than the survey data used in the previous XXL analysis.They obtained T 300 kpc = 4.9 +0.8  −0.6 keV, T R 500 = 4.5 ± 0.6 keV and an hydrostatic mass M 500 = 2.4 +1.7 −0.8 × 10 14 M ⊙ thanks to the NFW mass model fit to the data and a temperature profile resolved in four bins.We note that the temperature estimate of Paper XX is surprisingly much lower than that of Paper XLVIII, the later being more reliable given the better data quality.The mass estimate of Umetsu et al. (2020), which used Paper XX temperatures, is consequently lower than what would have been obtained with alternative temperature measurements.
Our targets have masses M 500 ∼ 2 × 10 14 M ⊙ according to XXL estimates.This is in rough agreement with the mass of XLSSC 102 as already measured in Paper XLIV, M 500 ∼ (1 − 2)×10 14 M ⊙ , depending on the analysis choices (see Table 4 of Paper XLIV).We note that the ACT masses are significantly larger than those from XXL, depending on the assumptions, especially for XLSSC 102.The main properties of the cluster sample are summarized in Table 1.

NIKA2
The clusters XLSSC 072, XLSSC 100, and XLSSC 102 were observed from January 2018 to February 2020 with the NIKA2 camera (Adam et al. 2018a) at the IRAM 30m telescope, under projects 179-17, 094-18, 208-18, 093-19, 218-19, and 076-20.The observation scheme and the data reduction are similar for the three targets.We refer the reader to Paper XLIV for further details, where the XLSSC 102 data at 150 GHz were already presented.
In brief, the beam was monitored using Uranus observations.Pointing corrections were checked using nearby quasars about every hour.The observing conditions were overall stable with average zenith opacity for the period (for the methodology, see Catalano et al. 2014).The absolute calibration uncertainty was estimated using the dispersion of the flux measured from the observations of Uranus that bracket the clusters observations.Table 2 summarizes the observational details for all three clusters, which are in line with the characteristics of the instrument given in Perotto et al. (2020).
In the case of XLSSC 072, part of the data (October 2018; 25% of the total observing time) could not be used with the standard calibration procedure due to a failure in the software control.For these data, the 150 GHz channel calibration was performed using the sufficiently bright radio source FIRST J021511.4-034309(or XXL-GMRT J021511.4-034309,∼ 30 mJy at 150 GHz), located about 3 arcmin west of the cluster X-ray peak, assuming a constant flux as measured over the rest of the observing time.At 260 GHz, the source was too faint for proper measurement and the data could not be used.The resulting absolute calibration accuracy at 150 GHz was estimated to be 30% for this subset, and the pointing accuracy was verified to be within a few arcsec using the dispersion of fluxes over the reliable scans used for cross-calibration.More information on the data validation can be found in Appendix A.
The data were reduced as described in Adam et al. (2015), by combining the individual detector time lines to remove the contribution from the electronic and atmospheric noise.The individual scan maps were checked and flagged based on the presence of large correlated noise residuals, prior to co-adding them using inverse variance weighting.The astrophysical signal filtering induced by the data reduction was estimated by injecting a fake signal in the data and comparing the input and output as a function of angular scale (the cluster signal is filtered at scales larger than the field of view of about 6.5 arcmin; see Adam et al. 2015 for details).The statistical properties of the noise were derived by computing the power spectrum of half-difference maps obtained by dividing the dataset into equal parts.Monte Carlo Notes.† Peak S/N at an effective resolution of 27 arcsec FWHM at 150 GHz.
(MC) noise realizations were computed using this power spectrum and preserving the noise standard deviation as a function of coordinates, following Adam et al. (2016).
In Figure 1, we present the NIKA2 surface brightness images of the targets at 150 and 260 GHz.All three clusters are well detected and show an extended decrement at 150 GHz, as expected for the SZ effect.At an effective resolution of 27 arcsec, the peak S/N is −9.7, −9.2, and −6.9, for XLSSC 072, XLSSC 100, and XLSSC 102, respectively.Given the 150 GHz data, the SZ signal is expected to peak at about 0.1 mJy/beam at 260 GHz and thus be well below the noise level.Several infrared and radio galaxies are also visible at both frequencies.They were identified and subtracted in the rest of the study, as discussed in Appendix B, and we do not expect their contamination to play a significant role in the presented work.

XMM-Newton
We analyzed the XMM-Newton data around the position of the three clusters of interest using the XMM-Newton Science Analysis Software (XMMSAS) v16.1.We used the pipeline developed for the XMM-Newton Cluster Outskirts Project (X-COP, Eckert et al. 2017) to analyze the data.Namely, we applied the standard event selection criteria by running the XMMSAS tasks emchain and epchain.We then filtered out regions of enhanced soft proton background using the mos-filter and pn-filter executables to create clean event lists.Then we extracted X-ray photon maps for each observation in the [0.5 − 2] keV bands by selecting all the valid events in the energy band of interest.We used the eexpmap task to extract effective exposure maps, which allowed us to take the telescope's vignetting into account.Finally, we used the mos-spectra and pn-spectra executables to create maps of the non-X-ray background by rescaling filter-wheel-closed data to the count rates measured in the unexposed corners of the telescope.Next, we stacked the extracted count maps, exposure maps and background maps for the three detectors (MOS1, MOS2, and PN).For more details on the data reduction technique, we refer the reader to Ghirardini et al. (2019).
We applied the aforementioned processing to all the observations of the survey and then co-added all the observations to create mosaic images of the entire XXL area.From the resulting mosaic images, we extracted cutouts centered on the clusters of interest.

Optical and near infrared
Optical and near-infrared data are used to compare the ICM distribution to that of the galaxies.Here, we followed the procedure presented in Paper XLIV, to which we refer the reader for details.In brief, we used the galaxy photometric catalogs from the Canada-France-Hawai Telescope Legacy Survey (CFHTLS, Gwyn 2012) and selected possible member galaxies based on Notes. (a) Paper XXVIII. (b) Paper XV. (c) Identified in Paper XLIV, near the galaxy number density peak; coordinates extracted using the HSC I filter (this work).
their photometric redshift, magnitude and type (elliptical).We then produced density maps using a Gaussian kernel, leading to a map resolution of 54 arcsec (FWHM of the Gaussian kernel).
We then subtracted the contribution from local background and normalized the maps in level of signal-to-background (S/B).To investigate morphology, we used 1000 MC realizations of the galaxy density maps per cluster field, computed by Poissonian realizations of an idealized model fit to the data.
As a complement, we used public data from the Hyper Suprime-Cam Subaru Strategic Program (Aihara et al. 2018(Aihara et al. , 2022) ) for visual purposes 3 .In particular, filters R, I, and Z were combined to produce color images of the three target cluster regions.
The locations of the brightest cluster galaxies (BCG) in XXL clusters were identified in Lavoie et al. (2016), hereafter XXL Paper XV, and Ricci et al. (2018), hereafter XXL Paper XXVIII, using slightly different criteria.As noted in Ricci (2018) the identified BCGs agree for XLSSC 072 but are different for XLSSC 100, for which we distinguished the two candidate BCGs.The BCG of XLSSC 102 is not reported in Paper XV.However, a second BCG associated with a sub-cluster was discussed in Paper XLIV.We extracted its coordinates using the HSC I filter.The coordinates of the candidate BCGs are listed in Table 3.

Dynamical state estimates
The ICM pressure profile and the Y SZ − M scaling relation are expected to depend on the cluster dynamical state (Arnaud et al. 2010;Yu et al. 2015).It is therefore essential to have information about the dynamical state if one wants to test the standard evolution of these SZ-related observables.In this section, we determined the dynamical state of the clusters according to their morphology and the comparison between the different tracers of the cluster components, including their centers.

Morphology
The morphology of the clusters is estimated qualitatively via individual SZ, X-ray, and optical data, which we used to trace the ICM thermal pressure, the ICM thermal density, and the galaxy population, respectively.The comparison of these different datasets is also informative given their different sensitivities to the cluster components (e.g., Komatsu et al. 2001;Rasia et al. 2013;Nurgaliev et al. 2013;Donahue et al. 2016;Cialone et al. 2018;De Luca et al. 2021).
Radio and submillimeter contaminating point sources have been subtracted from SZ images prior to analysis and it should have a negligible impact on the results (see Appendix B).Xray contaminating point sources (essentially AGNs) have been masked conservatively using a 30 arcsec radius aperture according to the catalog presented in Chiappetti et al. (2018), hereafter XXL Paper XXVII.The X-ray and SZ images have been smoothed to an effective resolution of 27 arcsec (FWHM), while the optical density map kernel is two times larger to ensure a sufficient S/B (54 arcsec FWHM).The R, I, and Z band HSC images were slightly smoothed and combined to produce a color image that visually highlights the cluster member galaxies, given their redshift.
In Figure 2, we present the multiwavelength view of XLSSC 072, XLSSC 100, and XLSSC 102.This follows the analysis and results presented in Paper XLIV, to which we refer the reader for a more detailed comparison that focuses on XLSSC 102 using a similar strategy.The three clusters are welldetected in all bands.The signal compares well both in terms of amplitude and extension for all sources, in agreement with them 3 https://hsc-release.mtk.nao.ac.jp/doc/ having similar masses 4 .Deviation from spherical symmetry is observed in all clusters, which suggests the presence of disturbances in the gas and galaxy distribution.The SZ and X-ray signals overlap well on large scales but may show a significant difference on smaller scales (see Section 3.2 for a quantitative comparison of the centroid and peak coordinates) that could indicate local compressions caused by merging events (as in, e.g., Adam et al. 2014).Both SZ and X-ray signals present relatively flat surface brightness distributions for all clusters, which is consistent with them being dynamically disturbed systems.Accordingly, we do not observe prominent X-ray peaks that would indicate the presence of a dense cool core associated with a relaxed system (Rossetti & Molendi 2010).As already investigated and reported in Paper XLIV, XLSSC 102 presents a bimodal galaxy number density distribution, while its ICM pressure and density are maximized in between the two peaks.It was interpreted as the result of two merging subclusters.Although the agreement between the galaxy and the gas distribution is better in XLSSC 100, a large offset is observed between the two, with the galaxy density extending more toward the southeast, possibly indicating that the gas is stripped in the direction of a passing subcluster.The ICM and the galaxy density match each other well in XLSSC 072 but they are both elongated in the east-west direction.The presence of multiple BCGs in XLSSC 100 (and possibly XLSSC 102) provides another indication of dynamical activity.In both clusters, one of the BCGs agrees well with the ICM location, while the other is largely offset, as expected for merging clusters with asymmetric mass ratios.Moreover, the elongation of the ICM is roughly aligned with the axis defined by the two BCGs, which  agrees with this scenario.In XLSSC 072, the BCG position is well aligned with the ICM and the galaxy distribution center.All three clusters present morphological signatures of large dynamical activity related to merging events.This is the case in terms of their properties seen in all individual bands for the three clusters and also given the difference observed in the ICM and the galaxy tracers for XLSSC 100 and XLSSC 102.The better agreement between these tracers, in the case of XLSSC 072, could be due to line-of-sight projection effects, in which case the merging event would be mostly oriented along the line-of-sight.According to the qualitative morphological study, we classify XLSSC 100 and XLSSC 102 as dynamically disturbed systems, and XLSSC 072 as likely dynamically disturbed.We note that in Paper XLVIII, XLSSC 072 is classified as disturbed according to the centroid shift estimate, but it would be classified as relaxed according to the BCG-X-ray offset, which is in good agreement with our findings.In Appendix D, we confirm these conclusions on the cluster dynamical state using the entropy profile, which is another excellent indicator of the ICM thermal state (e.g., Pratt et al. 2010).

Peak and centroid position
The measurement of the offsets between the peaks and centroid of the SZ, X-ray, galaxy density, and the BCGs positions is also an interesting way to address the cluster dynamical state (e.g., Lin & Mohr 2004;Hudson et al. 2010;Rossetti et al. 2016;Lopes et al. 2018;Zenteno et al. 2020;De Luca et al. 2021).To do so, we reproduced the analysis done for XLSSC 102 in Paper XLIV, but for the full sample.We estimate the peak as the coordinates of the maximum S/N of the signal, taken at the effective resolution of 27 arcsec (FWHM) for the SZ and X-ray data, and 54 arcsec (FWHM) for the galaxy number density.The centroid is obtained by fitting a 2D Gaussian function on the images 5 .Uncertainties are computed by running the same procedure on 1000 MC realizations of each data (except for the BCG coordinates, which have negligible uncertainties).We refer the reader to Paper XLIV for more details on the procedure.
The posterior likelihood distribution in the R.A.-Dec.plane are reported in Figure 3 for the peaks and the centroids, respectively (see also Table 4 for numerical results).As expected, better constraints on the centroid are obtained compared to the peak position.Given the uncertainties, the recovered peak and centroid coordinates of XLSSC 072 are consistent for all the data, with the only exception being the tension between the SZ and the X-ray centers, albeit only at a level of about 2σ.We note that the BCG is located at the intersection of the SZ and X-ray confidence intervals.In the case of XLSSC 100 and XLSSC 102, however, significant differences are observed between the ICM and at least one of the BCGs, both for the peak and the centroid.The other BCG is generally located closer to the X-ray peak and further from the SZ one, in agreement with the scenario in which a merger event induced a local boost of the pressure aside from the remnant denser core of the clusters, next to which the BCG is sitting.In the two clusters, the location of the centroid of the SZ and X-ray agree better than that of the peak, despite the smaller error bars.Again, this supports the fact that a merger event disturbed the cluster cores, while the ICM on large-scale was only weakly affected.In all cases, the uncertainties in the peak and centroid of the galaxy density distributions are too large to draw firm conclusions, although they agree well with the merger sce- 5 We note that the X-ray centroids roughly coincide with that of the XXL reference coordinates given the detection algorithm (see Table 1).
nario by tracking better the BCG that present the largest offset to the ICM.
The offsets between the different cluster components support XLSSC 100 and XLSSC 102 being merging clusters.They also favor -although the evidence is weaker-XLSSC 072 being dynamically perturbed.

Search for discontinuities and substructures
As a complementary investigation of the cluster dynamical state, we searched for substructure and discontinuities in SZ signal using the methodology developed in Adam et al. (2018b) and the Gaussian gradient magnitude and difference of Gaussian filtering.Given the limited S/N of the data and the compactness of the signal at these redshift and mass, we did not find any significant features in the data.This agrees with the results of Adam et al. (2018b), which state that the signature from merger event is difficult to identify at S/N ≲ 10.

Modeling and analysis procedure
This section presents the modeling and analysis methodology developed to extract the pressure profile and the location of our targets on the Y SZ − M relation.After discussing the SZ and Xray observables, we present the methodology used to extract the density profile, extract the pressure profile, and estimate masses using the hydrostatic equilibrium (HSE) assumption and scaling relations.In this work, we essentially considered HSE masses, despite the fact that our targets present indication for dynamical activity.This choice was motivated by the fact that the cluster pressure profile and the Y SZ − M relation, which we aim to test at high redshift and low mass, were calibrated based on HSE mass measurements (or scaling relations themselves calibrated with HSE masses, e.g., Arnaud et al. 2010 andPlanck Collaboration et al. 2014).Our approach is very similar to these works to ease the comparison.Moreover, given the mass and redshift of our targets, no other reliable individual mass estimates are available.
In Appendix D, we also discuss the reliability of the HSE assumption in light of thermodynamical indicators.

Sunyaev-Zel'dovich and X-ray observables
The SZ effect surface brightness is given by (Birkinshaw 1999) where I 0 is the CMB intensity.The characteristic SZ spectrum, f (ν), does only depend on the frequency in the nonrelativistic approximation, which applies well in the case of our sample.The amplitude of the distortion is given by the Compton parameter, which depends on the line-of-sight integration of the thermal electron pressure, P e , as where σ T is the Thomson cross-section and m e c 2 the electron rest mass.Given the NIKA2 beam and bandpass at 150 GHz, a Compton parameter y = 10 −4 corresponds to a surface brightness of ∆I ν = −1.19± 0.09 mJy/beam (Ruppin et al. 2018).The X-ray surface brightness is expressed as (Sarazin 1986) where n e is the thermal gas electron number density.The cooling function, Λ, depends weakly on the temperature, T , and on the metallicity, Z.

Extraction of the thermal electron density profile
The cluster electron density profiles are extracted as in Paper XLIV.In brief, we used the pyproffit package (Eckert et al.  2020) 6 , which is the Python implementation of the proffit soft-ware (Eckert et al. 2011).The X-ray surface brightness was extracted in radial bins of 5 arcsec width by accumulating photon counts within each annulus and correcting the vignetting by dividing by the local exposure map.As in Section 3, point sources from the XXL catalog were masked by excluding circles of 30 arcsec radius, corresponding to an encircled energy fraction of ∼ 90%.The multi-scale decomposition developed in Eckert et al. (2016) (hereafter XXL Paper XIII) was used to deproject the thermal electron number density profile assuming spherical symmetry.A single-temperature APEC model (Smith et al. 2001) absorbed by the Galactic N H was used to convert from X-ray count rate to emission measure, with temperature fixed to the ones from Paper III listed in Table 1.The model was convolved with the XMM-Newton point spread function and fitted to the data using a Poisson likelihood in PyMC with the No U-Turn Sampler (Salvatier et al. 2016).In the end, we obtained the best-fit electron number density profile together with 1000 realizations that we used to compute uncertainties.We note the presence of small differences compared to the profile presented in Paper XLIV.They are due to better modeling of the point spread function, accounting for the exact location of the cluster in the field of view and the combination of the multiple XXL pointings.The profiles are reported in Appendix C.

Extraction of the thermal pressure profile
We modeled the ICM thermal pressure via the MINOT software (Adam et al. 2020) using several different approaches described hereafter.MINOT allows us to easily produce SZ maps, given a pressure profile, projected on the same grid as the data.The maps are convolved with the NIKA2 beam and the transfer function that describes the filtering induced by the data reduction procedure (Adam et al. 2015).The surface brightness profiles are extracted in bins of 5 arcsec in width and up to a distance of 3 arcmin from the cluster center.As a reference, the XXL detection center is used, corresponding roughly to the X-ray centroid.
We discuss this choice in Section 5. Given a model of the pressure profile, the parameters are fitted to the data using a Markov chain Monte Carlo (MCMC) technique.The parameter space is sampled with the emcee package (Foreman-Mackey et al. 2013).
A Gaussian likelihood function was used to compare the model and the data.We account for the full covariance matrix, computed using MC realizations of the noise (see Adam et al. 2016 for the procedure).In addition to the NIKA2 data, we also impose a Gaussian prior on the total SZ flux (see Equation 12, integrated up to 5R 500 ).To do this, we use the Planck measurement obtained by fitting a Gaussian function with a 10 arcmin FWHM (i.e., the Planck beam size) on the Compton parameter map (Planck Collaboration et al. 2016b) at the location of the unresolved targets (see Table 5).As for the density profiles, we propagate the uncertainties on the pressure profile using 1000 pressure profile models randomly taken from the MCMC chains.
We consider the two following ways of modeling the pressure profile (see also Section 4.4 for the direct modeling of the mass profile, which is also an alternative way of modeling the pressure profile indirectly).
The first is the Generalized Navarro-Frenk-White model.As a baseline, the pressure is described according to the generalized Navarro-Frenk-White (gNFW) model (Nagai et al. 2007), expressed as The parameter P 0 is a normalization, r p is a characteristic radius, and c, a, and b describe the slope of the profile from the core to the outskirts.The fit parameters include all the pressure profile parameters from Equation 4plus the map zero level as a nuisance parameter.We use flat priors on the normalization and scale radius (P 0 > 0 and 5R 500 > r p > 0, with R We also consider the modeling of the pressure by fixing the normalization of the profile at given radii (see, e.g., Ruppin et al. 2017, for a similar method).This is the binned model.The full profile is then interpolated in logarithmic space before line-ofsight integration and projection, as implemented in MINOT.We define the values of the radii as five bins logarithmically spaced from 50 kpc to 1 Mpc, r i ≡ [50,106,224,473,1000] kpc.This allows us to sample the profiles where NIKA2 is sufficiently sensitive and obtain reliable constraints in each bin.The five pressure model parameters are given by P i ≡ P(r i ), to which the zero level of the map is added as a nuisance parameter.The pressure parameters are restricted to verify P i > 0 in all bins i.Assuming that the ICM is in hydrostatic equilibrium and spherically symmetric, the total mass enclosed within the radius r is given by where µ gas = 0.61 is the gas mean molecular weight, m p is the proton mass, and G the Newton constant.We extract the hydrostatic equilibrium mass profile using the two following methods.First, the HSE mass profile given in Equation 5is obtained by combining the density and the pressure profiles measured independently from the X-ray and SZ data, as discussed in Sections 4.2 and 4.3, respectively.Given the critical density of the Universe at the cluster's redshift, we derive the overdensity contrast by integrating the mass profile, which we use to obtain R 500 and thus compute M HSE,500 .Uncertainties are propagated from the pressure and the density profiles by combining 1000 model realizations randomly taken from the MCMC chains.
Alternatively, we directly model the total mass density profile as the sum of the gas density, which we know from Xray data (see Equation D.4), and a Navarro-Frenk-White (NFW, Navarro et al. 1996) model to describe the other components (essentially the dark matter).We note that in practice, we have checked that modeling the total mass with a single NFW model does not significantly affect our results since the gas is absorbed in the NFW component.The NFW density model is written as In the case where the NFW model accounts for the total mass (gas included), the characteristic radius can be simply written as r s = R 500 /c 500 and with ρ c being the critical density of the Universe at the cluster's redshift and c 500 the concentration.The enclosed hydrostatic mass is given by Taking advantage of the MINOT code implementation, we combine this model with the density profile and Equation 5 to compute the pressure profile model, P e (r) = P e (r 0 ) + r 0 r Gµ gas m p n e (r ′ )M HSE (r ′ ) with r 0 being a radius taken sufficiently far, at which point the pressure is negligible 7 .The pressure profile model is then compared to the NIKA2 data as in Section 4.3.However, in this case, the density profile is randomly sampled from the 1000 MC realization available at each step of the MCMC to account for the associated uncertainty.The fit parameters that we use are the mass M HSE,500 and the concentration c 500 , which are related to ρ 0 and r s .While the main goal of this method is to directly describe the mass profile with a physically motivated model, this also provides an alternative pressure profile model that complements the methods described in Section 4.3.We refer the reader to Eckert et al. (2022) and Muñoz-Echeverría et al. (2023), for example, for a detailed description of this approach.
4.5.Mass estimates from scaling relations and fitting the universal pressure profile In addition to direct mass measurements based on the HSE assumption, it is useful to compute and compare masses estimated using the global cluster properties that are usually easier to obtain.We thus also use the UPP normalization as a mass proxy.

Mass estimation from the universal pressure profile
The UPP (e.g., Arnaud et al. 2010, which we follow here, or any other calibration of the profile) depends exclusively on the cluster mass and redshift.In this scenario, Equation 4can be expressed as with P 500 ≡ P 500 (M 500 ) being the self-similar normalization (Nagai et al. 2007) and f M = M 500 3×10 14 M ⊙ 0.12 a small mass dependence correction.Following Arnaud et al. (2010), the parameters of the profile are set to (P 0 , c 500 , a, b, c) = (8.403,1.177, 1.0510, 5.4905, 0.3081).We use Equation 10 as in Section 4.3 to fit the NIKA2 data with the mass M 500 and the map zero level as the only free parameters.This is similar to the methodology used by Hilton et al. (2018Hilton et al. ( , 2021) ) to extract ACT masses.Given the fact that the clusters appear as dynamically active, we also reproduce this work by using the mean profile of morphologically disturbed clusters from Arnaud et al. (2010).We note that the UPP was calibrated using nearby clusters and assumes standard evolution, which is what we aim to test in the present work, as we discuss in Section 5.Although they are not directly obtained from the HSE assumption, the masses used in the UPP calibration were obtained from the Y X − M relation, itself calibrated using the direct HSE masses of relaxed clusters (Arnaud et al. 2007).

Mass estimates from the Y SZ − M relation
The SZ flux Y SZ,500 is tightly correlated with the mass.The bestfit scaling relation as calibrated in Planck Collaboration et al. ( 2014) is given by .(11) It was obtained using masses derived using the Y X − M relation, itself calibrated using HSE masses computed from X-ray observations.This relation is used to estimate the mass according to our SZ flux measurement.To do so, we compute the spherically integrated SZ flux given the pressure profile as Equation ( 12) is integrated up to R 500 to obtain Y SZ,500 , and thus depends on M 500 .Therefore, we perform the measurement by iterating about the scaling relation (convergence is obtained within less than 1% after a few iterations).The full probability distribution in the Y SZ − M plane is obtained by repeating the measurement with 1000 pressure profiles randomly taken from the MCMC chains.By default, we use the pressure profile obtained from the gNFW fit to the data.We note that although it is possible to use Equation ( 11) to estimate the cluster's mass, this relation, as calibrated using nearby clusters, is precisely what we aim to test in the present work.We discuss it in Section 5.

Mass estimates from the Y X − M relation
The X-ray analog of the SZ flux, Y X,500 ≡ k B T X M gas,500 (Kravtsov et al. 2006), is an excellent mass proxy (Arnaud et al. 2007;Vikhlinin et al. 2009a).Here, we use the best-fit scaling relation from Arnaud et al. (2010) in order to obtain a high-quality mass estimate based on direct X-ray-only measurement: The masses used to calibrate this relation were obtained applying the HSE on relaxed clusters observed in X-ray.We estimate Y X,500 using the measured X-ray temperatures listed in Table 1 (Paper III).The gas mass profile is computed from the gas density profile as with µ e = 1.16 being the electron mean molecular weight.As for Y SZ,500 , the integration is performed up to R 500 to obtain M gas,500 .
Since the estimate of Y X,500 depends on R 500 , and thus M 500 , we perform the measurement by iterating about the scaling relation.
The full probability distribution, and thus the uncertainty on the mass, is obtained by repeating the measurement with 1000 density profiles taken from the MC realizations and simultaneously sampling the temperature within its uncertainty.Again, as discussed in Section 5, we note that Equation (13) was calibrated using nearby clusters and assumes standard evolution, which is what we aim at testing in the present work.

Results and discussions
In this section, we present the results of the analysis.After discussing the mass measurements, we focus on the pressure profile and the Y SZ − M scaling relation.

Direct HSE mass profiles
Figure 4 presents the HSE mass profiles obtained either from combining the density profiles with the gNFW pressure model, or directly modeling the mass with an NFW profile.The NFW model leads to the smallest uncertainties due to fewer parameters involved, but the two methods show excellent agreement over the full radial range, highlighting the robustness of the measurement.The quality of the recovered profiles is remarkable given the low masses and the high redshifts of these clusters.
The three clusters present comparable mass profiles.They flatten at r ≳ R 500 for XLSSC 072 and XLSSC 102, but it keeps rising for XLSSC 100.This feature is due to XLSSC 100 having a flatter outer pressure profile and a slightly steeper outer density profile than the other two clusters.It implies significantly larger uncertainties on the recovered value of M 500 for XLSSC 100 than XLSSC 072 and XLSSC 102.The numerical results on the mass are reported in Table 6, where we also give the corresponding SZ flux.We refer to Appendix D for further investigation of the reliability of the mass profile using thermodynamics diagnosis.

Comparison between direct measurements, estimates from scaling laws, and the literature
The masses derived with the methods presented in Section 4.4 are listed in Table 6.Only the direct HSE masses are independent of any calibration at low redshift.We also report the SZ flux enclosed within R 500 .In Figure 5, we compare the masses derived in the present work to those obtained in the literature (Table 1).The masses obtained from XXL scaling relations reflect the large uncertainty in the mass proxies on which they rely and the precision of the scaling relation.This is also the case for our Y X − M and Y SZ − M masses, although they rely on more precise mass proxies given the data in hand and on scaling relations that are expected to be very tightly related to the mass.We note that these relations are calibrated using X-ray data (Arnaud et al. 2010), but we do not correct for any mean HSE bias here.The masses derived through the UPP normalization should match the Y SZ − M ones perfectly in the case of a perfect UPP profile.The HSE masses that we derive are the only direct measurements.However, they are affected by systematics in the modeling and by the hydrostatic mass bias.In Appendix D, we discuss possible biases in the recovered masses in light of thermodynamics diagnosis.
In Figure 5, we observe a very good general agreement between the different mass measurements despite the very different methodologies and assumptions involved.Only XLSSC 102 presents a ≳ 2σ tension between the direct HSE masses and the masses obtained from the Y X − M and Y SZ − M scaling relations and the UPP normalization fit.It could be due to the ongoing merger activity that affects the HSE assumption.In Appendix D,  1 and Table 6.The first block (gray points) corresponds to survey measurements (from XXL and ACT), the second block to masses derived using low-redshift, higher mass calibration proxies in the present work (purple points), and the last block to direct HSE measurements from the present work (magenta points).We also report the value from Paper XLIV for XLSSC 102 when assuming a similar method and center.The small difference is mainly due to the updated X-ray density profile used in the present work.using the temperature reported in the detailed analysis of Paper XLVIII instead of the one from Paper III, when computing Y X (see also Appendix D).The UPP-based masses that we derive agree very well with the Y SZ − M masses, which indicates that the shape of the pressure profiles does not strongly deviate from that of the UPP.However, an ∼ 2σ tension between our UPPbased mass and that obtained from ACT data (Hilton et al. 2018) is observed for XLSSC 102 despite the similar methodology employed.When changing the UPP to morphologically disturbed or cool-core models, the changes in the mass are about 1σ.
In the following, we use these masses to compare our NIKA2 measurements with the pressure profile and Y SZ − M expected from standard evolution.Given the precision in the masses and the underlying assumptions that they involve, we use the Y X − M and the direct HSE (NFW-based) masses for reference.

Pressure profile
The SZ surface brightness profiles of XLSSC 072, XLSSC 100, and XLSSC 102 and their corresponding pressure profiles are shown in Figure 6 for the case of the gNFW model.The SZ decrement is detected up to about R 500 for XLSSC 072 and XLSSC 100.The S/N is slightly lower for XLSSC 102.The best-fit models describing the data and their uncertainties are obtained as discussed in Section 4. We report the 68% interval allowed by the data as a gray band.While we use the full covariance matrix in the analysis, the uncertainties only provide the diagonal of the covariance matrix.We refer the reader to Paper XLIV for more details about the computation of the covariance matrix.The residuals between the best-fit model and the data are provided in Appendix E at the map level.
The comparison with the models from Arnaud et al. (2010), namely the UPP, the mean morphologically disturbed profile, and the mean cool-core profile, is performed using Y X − M and direct NFW masses.The surface brightness models have been convolved with the instrument response function.The derived pressure profiles reflect the same behavior but are deconvolved from the instrument response and projection effects.All the measured profiles agree best with the morphologically disturbed model in terms of shape.They are significantly different than the averaged cool-core pressure profile, but they still agree with the UPP within error bars.The choice of the mass is highly relevant to the comparison in terms of amplitude.For instance, the models describing XLSSC 102 reach 2σ lower when using direct HSE mass measurements, while the match is excellent with the Y X − M mass.Similarly, a much better match would be obtained for XLSSC 072 by using the temperature from Paper XLVIII instead of the one from Paper III to compute Y X .As already mentioned regarding the mass profile, XLSSC 100 presents a pressure profile outer slope that is more shallow than expected from the models at the ∼ 2σ level.
Assuming standard evolution, the shape of the pressure profile is in excellent agreement with the dynamical state analysis of Section 3. It suggests that XLSSC 072, XLSSC 100, and XLSSC 102 are increasingly disturbed systems, with even XLSSC 072 showing evidence of disturbance.Alternatively, given the prior knowledge of the cluster dynamical states inferred in Section 3, the data are in good agreement with the pressure profile calibrated on low-redshift clusters (Arnaud et al. 2010) and scaled to low mass and high redshift using standard evolution.The agreement would be even better when accounting for the intrinsic scatter in the expected profile.
In Figure 7, we compare the pressure profiles recovered using the three methods described in Section 4. Despite the very different methodologies, all profiles show excellent agreement within uncertainties at all radii.In the case of these systems, the NIKA2 data are most sensitive to the pressure profile in the range from about 100 kpc to 600 kpc.It is remarkable that reliable constraints on the pressure profile can be obtained nearly up to 2R 500 for such high-redshift and low-mass clusters. .NIKA2 constraints on the thermal pressure profile.Top: gNFW constraints on the SZ surface brightness.Middle: gNFW constraints on the thermal electron pressure profile.The gray band provides the 1σ constraint on the model.The green, blue, and red lines give the expected model according to the UPP, the cool-core (CC) pressure profile, and the morphologically disturbed (MD) pressure profile according to Arnaud et al. (2010) given the Y X − M masses.The corresponding dashed lines give the same model assuming higher or lower masses by 1σ.Bottom: Same as the middle row, but computing the expected models given the masses derived from the direct NFW fit to the pressure profile plus the density profile.20 local clusters.The masses used for the calibration are thus HSE masses, as given in Equation 13.The measured intrinsic scatter (7%) is reported as the gray band on the figure, together with the best-fit relation (Equation 11).For comparison, we also show the location of other NIKA and NIKA2 observed clusters (0.5 < z < 0.9, Adam et al. 2015Adam et al. , 2016;;Ruppin et al. 2017Ruppin et al. , 2018;;Kéruzoré et al. 2020) on the relation, but we stress that the flux and masses were not derived in a homogeneous way for those.The masses (and thus SZ fluxes, see Section 4) used for the XXL sample are either direct HSE masses obtained from the combination of the NIKA2 and XMM-Newton data, or those obtained using the Y X − M relation from Arnaud et al. (2010).Therefore, in the latter case, we implicitly tested the Y SZ − Y X relation at high redshift and low mass.

The Y SZ − M relation
As we can observe, our sample sits at the low-mass end of the Planck Collaboration et al. ( 2014) calibration sample, but our clusters are located at redshift z ∼ 1 instead of z ∼ 0.2.Nonetheless, thanks to the quality of the data, we were able to obtain comparable uncertainties on the flux, but uncertainties on the mass remain larger by a factor of two or more.Despite the different regime that we probed and the fact that these clusters are significantly disturbed (implying a likely higher intrinsic scatter, Yu et al. 2015), the XXL clusters follow the scaling relation remarkably well when using the Y X − M relation to obtain the mass.This is also the case for other NIKA and NIKA2 clusters with published SZ fluxes and masses, at higher masses and lower redshifts.When using direct HSE mass measurement, only XLSSC 102 deviates from the relation by about 2σ.However, as investigated in detail in Paper XLIV and Appendix D, this may be related to systematic uncertainties in the mass measurement due to the very complex morphology and dynamical state of this cluster or a very large hydrostatic mass bias.
Either way, we do not observe any significant deviation from the Y SZ − M scaling relation in the three high-redshift, lowmass XXL clusters.Moreover, our results implicitly show that the three clusters follow the Y SZ − Y X relation remarkably well.While the size of our sample does not allow us to infer statistical conclusions on the relation, this provides a first indication of these relations being stable down to M 500 ∼ 2 × 10 14 M ⊙ and z ∼ 1.

Discussion
The results presented in the paper may be affected by the analysis choices, which we discuss here.For instance, we used the XXL detection center as the reference for extracting the profiles and derived quantities.While not much freedom is available for XLSSC 072 given the agreement between the different cluster components on the center, this is not the case for XLSSC 100 and XLSSC 102.Paper XLIV explored the systematic uncertainty associated with this choice for the most perturbed cluster of our sample, XLSSC 102.Hence, this provides us with an upper limit on this uncertainty for the sample, which is in fact modest for the global quantities (M HSE,500 , Y SZ,500 ) -on the order of 1  2 σbut it can be as high as about 1σ for the profiles in the central region.Similarly, the morphological analysis showed that the clusters are not spherically symmetric.By performing the analysis in different sectors for XLSSC 102, Paper XLIV es-timated the corresponding dispersion to be ≳ 1σ on the profiles, but slightly smaller on global quantities.
Our analysis also relies on the modeling of the pressure profile or the HSE mass profile.Nevertheless, we tested that different methodologies relying on very different assumptions led to consistent results.Therefore, the systematic uncertainty associated with the modeling is expected to be much smaller than statistical uncertainties.The direct HSE masses that we derived rely on spherical symmetry and the HSE assumption and are more likely to be affected by systematic effects (cluster geometry, clumping, etc.; see Appendix D) than indirect methods, but they are the only direct measurement that can be used to test the Y SZ − M relation.On the other hand, the most robust and precise masses are likely to be the ones derived from the Y X − M relation, the Y SZ − M relation, or the UPP normalization fit, but they rely on the calibration of these relations at lower redshifts and higher masses, which is what we aimed to test here.Moreover, these methods do not generally propagate the intrinsic scatter associated with the scaling relation or pressure profile.
The main conclusion of our work is that the pressure profile and the SZ mass proxy are in line with standard extrapolation down to z ∼ 1 and M 500 ∼ 2 × 10 14 M ⊙ .However, this is based on a sample of only three clusters.Stronger conclusions would require increasing the size of the sample, but this might require a significant amount of observing time for such masses and redshifts, which is not straightforward to obtain.Despite the limitation of the present work, the results indicate that the physics that drives cluster formation is already in place in the regime that we explored.

Summary and conclusions
The SZ structure of the ICM gives us precious information about the thermal state of galaxy clusters and the astrophysical processes at play during their formation.This is reflected in the cluster thermal pressure profile and the scaling relation maintained by the SZ flux and their mass.Detailed investigations of these properties have been done at low redshift, and the effort is being put in at high redshifts for massive clusters.However, at high redshifts and low masses, where the largest deviations from self-similarity are expected, the investigation of the SZ structure with resolved data has remained nearly unexplored to date.
In this paper, we present the analysis of three XXL-selected clusters at z ∼ 1 and M 500 ∼ 2 × 10 14 M ⊙ observed with the NIKA2 camera via their SZ signal, at a resolution of about 18 arcsec.We investigated the dynamical state of the sources using SZ, X-ray, and optical data.We extracted their pressure profile and compared them to expectations from standard evolution.Complementary X-ray data were used to extract the gas density profile, which we combined with the pressure to measure the hydrostatic masses of the systems.We also estimated the masses using the UPP normalization fit to the data, Y SZ − M, and the Y X − M scaling relations.
The main conclusions of this work are listed here.
-The three clusters, XLSSC 072, XLSSC 100, and XLSSC 102, at z ∼ 1 and M 500 ∼ 2 × 10 14 M ⊙ are well detected with NIKA2 in about hours hours per source.The signal is extended and the peak S/N reaches −9.7, −9.2, and −6.9, respectively.These are among the first resolved SZ data available down to such low masses and high redshifts.-All three clusters present evidence for ongoing merging activity.This is shown by their disturbed morphologies that present deviation from a compact, spherically symmetric distribution.In the case of XLSSC 100 and XLSSC 102, this is further confirmed by the offset between the peak and centroid of the SZ, the X-ray, the galaxy density, and the BCGs.Assuming standard evolution, this is also confirmed by the flatness of their pressure profiles.-The pressure profile is well constrained up to 2R 500 , which is a remarkable achievement given the low masses and high redshifts of the clusters.-The pressure profile of the three clusters agrees with that of local dynamically disturbed systems, once rescaled according to standard evolution in mass and redshift.In the case of XLSSC 072, the data are in better agreement with expectations from dynamically disturbed systems, but they also agree with the UPP.-Despite their perturbed ICM, their low masses, and high redshifts, we do not find any significant deviation in the Y SZ − M scaling relation followed by our targets.-The comparison of the pressure profile and the Y SZ − M scaling relation to that of local samples is limited by uncertainties in the mass.This highlights the difficulty of obtaining accurate and robust mass estimates in this new regime.
Galaxy cluster formation is primarily driven by gravity, on top of which feedback processes help regulate cluster evolution.This includes shock heating, turbulent cascade of energy injected from large-scale structures accretion and mergers, and supernova and AGN feedback.These processes are expected to shape the radial thermodynamical profiles and scaling relations followed by galaxy clusters.The lack of significant nonstandard evolution in the pressure profile and the Y SZ − M relation of clusters when extrapolating those expected from lower redshift and more massive objects suggests that the dominant mechanisms that drive clusters' observational properties are already in place around z ∼ 1, down to M 500 ∼ 10 14 M ⊙ .
work is based in part on data products produced at Terapix available at the Canadian Astronomy Data Centre as part of the Canada-France-Hawaii Telescope Legacy Survey, a collaborative project of NRC and CNRS.This paper present data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by Subaru Telescope and Astronomy Data Center at National Astronomical Observatory of Japan.Data analysis was in part carried out with the cooperation of Center for Computational Astrophysics, National Astronomical Observatory of Japan.
the detection parameters, it is estimated to be 0.88 at 150 GHz and 0.93 at 260 GHz.Once the source catalogs are made in each band independently, we measure the flux of the counterpart band by fitting a source at the location of the detection.We match the catalogs with other bands to associate the detected sources using an aperture of one beam FWHM.We also list the possible matches within the same band of two nearby sources if they fall within a single beam.
The list of NIKA2 identified sources in both bands, within 5 arcmin of the cluster centers, are reported in Table B.1 for all clusters.In Figure B.1, the source positions are reported on the images.

B.2. Submillimeter contamination
We compute the mean 150 GHz to 260 GHz flux ratio for all the sources detected at 150 GHz, excluding radio sources, F 150 /F 260 = 0.221 ± 0.014.Using this reference value, we find that the expected S/N is nearly the same at 150 and 260 GHz, implying that any source that could significantly bias the SZ signal should be detected at 260 GHz.If instead we use the mean ratio for the sources detected at 260 GHz, F 150 /F 260 = 0.123 ± 0.008, the S/N should be two times larger at 260 GHz, and the potential bias would be reduced.Accordingly, and given the fact that no hint of a 260 GHz source is visible in the SZ region for the three clusters, we do not expect that any significant submillimeter source that is blended in the SZ signal could be missed and significantly bias the clusters analysis.
The field of XLSSC 072 was observed with Herschel/SPIRE (Griffin et al. 2010) at 500, 350, and 250µm (obsID 1342189031 and 1342190313), and we compare the SPIRE maps to NIKA2 images in Figure B.2.Although the image depth is relatively shallow, the brightest regions match the NIKA2 sources well.This confirms that no bright submillimeter source is missed by NIKA2 and that the SZ signal from XLSSC 072 is not contaminated.

B.3. Radio GMRT counterparts
While NIKA2 260 GHz data can be used to assess the contamination from submillimeter sources, radio data are necessary to address the contamination from radio galaxies.We use XXL/GMRT images and catalogs from Smolčić et al. (2018), hereafter XXL Paper XXIX, to do so.In Table B.2, we list all the GMRT sources identified in the 10 arcmin × 10 arcmin region around the three clusters.The GMRT positions uncertainties are at most 0.5 arcsec, which is negligible for our purpose.In addition, we provide estimates of the fluxes expected at 150 GHz assuming a power-law spectrum: . This is done by using a spectral index equal to the mean value of the sample selected in Paper XXIX, between 610 and 1400 MHz (or its mean value plus one standard deviation).This provides an upper limit of the expected flux since a steepening of the radio spectrum is common at higher frequencies.Some of the GMRT sources are also detected in the NVSS surveys, in which case a spectral index estimate is available for individual sources, which we use to compute a more reliable flux estimate at 150 GHz.Given these estimates, only XXL-GMRT J021511.4-034309 in the field of XLSSC 072 was expected to be detected and is indeed detected (NIKA2-150 J021511.5-034308, NIKA2-260 J021511.6-034309).Figure B.1 shows the GMRT maps and how they compare to the NIKA2 data.We also use the FIRST catalog, which we compare to the GMRT images.We note that it is very un-

Appendix D: Thermodynamic profile diagnosis
This appendix presents the temperature, entropy, and gas fraction profiles of our target clusters, which we use as a thermodynamic diagnosis of the dynamical state and to address systematic effects in the mass measurement.They are computed within the framework of the MINOT software (Adam et al. 2020) given the pressure and the density inferred from the SZ and X-ray data.As a reference, we use the pressure profile inferred from the gNFW model fit to the data.We note that detailed discussions regarding the recovered temperature, entropy and gas fraction were presented in Paper XLIV for XLSSC 102.The thermodynamic profiles reported here are in agreement with our previous analysis, but they differ slightly due to the updated density profile.
The entropy and temperature are given by The three entropy profiles are compared to the self-similar baseline, in the case where only gravitational effect are present (Voit et al. 2005), given by K(r) ∝ (r/R 500 ) 1.1 .The masses used for the comparison are obtained from the Y X proxy.We can observe that all clusters present a large excess entropy in the core, within r ≲ 300 kpc, of about 300 keV cm 2 .Such a feature is typical for disturbed systems (Pratt et al. 2010).This indicates that all three systems are dynamically disturbed, most likely because of the presence of merging events, in agreement with our imaging analysis (see Section 3).We note that XLSSC 102 agrees with a flat entropy profile at all scales.Beyond r ∼ 300 kpc, the profile is even lower than the self-similar baseline, although uncertainties are becoming very large.As the self-similar baseline corresponds to a minimal heat injection from gravitational collapse, such a feature is not expected.This could indicate an excess in the density profile caused by inhomogeneities in the gas, as observed for other nearby clusters with high-quality data (Tchernin et al. 2016), which would also bias low the HSE mass estimates.We note that this feature is reduced when using the direct HSE mass measurement since they are lower, but it does not entirely disappear.
The three temperature profiles agree with that of merging systems, with a profile decreasing from the core to the outskirt.We also report the projected temperature measured within 300 Given the uncertainties and the fact that the two measurements are not strictly comparable, good qualitative agreement is observed for XLSSC 100 and XLSSC 102.On the other hand, our SZ plus X-ray-derived temperature for XLSSC 072 is in qualitative agreement with the one from Paper XLVIII, higher than but still comparable to the one from Paper III, but in significant disagreement with the one from Paper XX.As the Y X -based masses are using Paper III temperatures, we conclude that while no significant issue is observed with XLSSC 102 and XLSSC 100, the Y X -derived mass for XLSSC 072 might be biased depending on the choice of the temperature measurement.For instance, using our SZ plus X-ray measurement would increase the mass within ∼ 2σ.
The gas fraction profiles increase from the center to the outskirts, in agreement with the expected baryon depletion generally expected in the center.As we can see from Equation D.3 and D.5, the gas fraction depends on the hydrostatic mass bias, which is set to b HSE = 0 here.On the other hand, the universal gas fraction at R 500 is given by f gas,univ (R 500 ) = Ω b Ω m Y b,500 − f ⋆,500 , (D.6) with Y b,500 ≃ 0.85 being the baryon depletion factor and where f ⋆,500 ≃ 0.015 accounts for the baryons that have condensed into stars (see Eckert et al. 2019 for details).Assuming that the pro-files should reach the universal gas fraction at R 500 , it is possible to estimate the value of b HSE .As we can see, the three profiles reach the cosmic baryon fraction at about R 500 .While XLSSC 072 and XLSSC 100 are in rough agreement with the expected value, XLSSC 102 presents a gas fraction that is higher by about 2σ.This might indicate that this system presents a high hydrostatic mass bias, in agreement with the fact that it is the most disturbed of our targets and already indicated thanks to the entropy profile.More quantitatively, a value of 1 − b HSE of 0.85, 1.0, and 0.5 would bring XLSSC 072, XLSSC 100, and XLSSC 102 to the expected universal gas fraction value at R 500 , respectively.
Given the entropy, temperature, and gas fraction profiles, we conclude that all three clusters are dynamically disturbed.Additionally, we find that the Y X -derived masses may be biased low for XLSSC 072 (within 2σ), and that direct HSEderived masses may be biased low by up to a factor of two for XLSSC 102.

Figure 1 .
Figure 1.NIKA2 images at 150 GHz (top) and 260 GHz (bottom) for XLSSC 072, XLSSC 100, and XLSSC 102, from left to right.The maps have been smoothed to an effective resolution of 18 and 27 arcsec at 260 and 150 GHz, respectively.S/N contours are shown with 1σ spacing starting at ±3σ.Data where the noise is greater than three times the value of the noise at the center of the map are masked.The effective beam size is shown in the bottom left corner of each panel.

Figure 2 .
Figure 2. Comparison of SZ, X-ray, and optical data for XLSSC 072 (left), XLSSC 100 (center), and XLSSC 102 (right).First row: Point-source-subtracted 150 GHz SZ surface brightness images with S/N contours.Second row: X-ray surface brightness images with S/N contours.Point sources from Paper XXVII have been masked.Third row: CFHTLS derived galaxy density images, Σ. Fourth row: HSC color images made by combining the R, I, and Z filters.In all panels, the cyan cross indicates the reference centers and the BCG positions are indicated as white hexagons.The gray dashed circles indicate the characteristic radii θ 500 estimated via the Y X − M scaling (see Section 5).In all panels, the black S/N (or S/B) contours start at 2σ and are separated by 1σ each.Magenta contours correspond to the SZ S/N, starting at 3σ and separated by 2σ each.We note that the XLSSC 102 data were already reported in Paper XLIV.

Figure 3 .
Figure 3. Probability distribution of signal peak (top) and centroid (bottom) location with respect to reference center in three wavelengths, for XLSSC 072, XLSSC 100, and XLSSC 102, from left to right.The BCG coordinates are indicated by the black stars.Contours give the 68% and 95% confidence interval.We note that we recover a posterior distribution that is in good agreement with that reported in Paper XLIV for XLSSC 102.
500 measured from the Y X − M relation; see Section 4.5) and Gaussian priors on the slope parameters based on Planck Collaboration et al. (2013): µ a,b,c = [1.33,4.13, 0.31] and σ a,b,c = [1.00,3.10, 0.23], corresponding to 75% of the mean value (i.e., 3 times larger than the values used in Paper XLIV, to ensure more freedom in the profile shape).

Figure 4 .
Figure 4. HSE mass profiles of XLSSC 072, XLSSC 100 and XLSSC 102.Left: Mass profile obtained using the gNFW pressure model together with the electron density profile.Right: Mass profile obtained by direct NFW mass modeling together with the electron density profile.For reference, the vertical dashed lines represent R 500 estimated using the Y X −M relation.The corresponding probability density functions for M 500 are shown as insets in the figures.The shaded region gives a 68% confidence interval.

Figure 8
Figure 8 compares the Y SZ − M relation followed by XLSSC 072, XLSSC 100, and XLSSC 102, to the Planck calibration sample used to derive cosmological constraints (Planck Collaboration et al. 2014).The Planck Collaboration et al. (2014) relation was obtained from a sample of z < 0.45 clusters, with a mean redshift ⟨z⟩ = 0.19.The masses were derived using the Y X − M relation from Arnaud et al. (2010), calibrated using an X-ray sample of

Figure 7 .Figure 8 .
Figure 7.Comparison of pressure profile as measured with different methods: gNFW modeling of the pressure, binned pressure profile, and NFW modeling of the mass with the joint use of the density profile.The vertical dashed lines give the location of R 500 as in Figure 5.

Figure B. 2 .
Figure B.2. Herschel/PACS images of XLSSC 072 at 500, 350, and 250 µm (from top to bottom).Black contours give the S/N in units of 1σ, starting at 3σ.The gray circles in the bottom left corner show the PACS beam FWHM in each band.White dashed contours show the NIKA2 150 GHz contours at -3, -5, -7, and -9 σ.Magenta and cyan contours show the 3, 4, and 5σ NIKA2 S/N at 150 and 260 GHz, respectively.The blue and red crosses indicate the point sources identified at 150 and 260 GHz in the NIKA2 data.
gas fraction is computed asf gas (r) = M gas (r) M tot (r) , (D.3)where the gas mass is obtained by integrating the electron density profile asM gas (r) = 4πr 0 µ e m p n e (r ′ )r ′2 dr ′ , (D.4) with the mean molecular weight µ e ≃ 1.15 and m p is the proton mass.The total mass is related to the hydrostatic mass via M tot (r) = M HSE (r) (1 − b HSE ) , (D.5) where b HSE is the hydrostatic mass bias.The entropy, temperature and gas fraction profiles are shown in Figure D.1 for the three clusters.

Figure C. 1 .
Figure C.1.Thermal electron density profile, from left to right, of XLSSC 072, XLSSC 100, and XLSSC 102.The solid lines provide the best profiles and the 68% confidence interval.The 1000 MC realizations are also provided via transparent markings to show the dispersion.

Table 1 .
Summary of the survey properties of the XXL targeted clusters.The masses obtained according to the mass-temperature relation are labeled M 500,MT .The value of M 500,MT based on Paper XX was computed from R 500 accounting for the different cosmological model.

Table 2 .
Observational summary of the NIKA2 observations after data selection.

Table 3 .
Coordinates of the candidate BCGs.

Table 4 .
Projected physical offsets between the measured peaks and best-fit centroids of different gas and galaxy tracers.The quantity Σ refers to the galaxy density.In the case of XLSSC 100 and XLSSC 102, the BCG 1 and BCG 2 are given in parentheses, respectively.

Table 5 .
Planck prior on the total flux.

Table 6 .
SZ fluxes and masses.The median value and the 68% confidence interval are reported.All the masses are assimilated to hydrostatic masses (see Section 4 for the different methodologies).Comparison between different mass measurements reported in Table Arnaud et al. (2010)asses.This indicates that the Y X − Y SZ relation followed by our targets is in excellent agreement with the one measured inArnaud et al. (2010).In the case of XLSSC 072, the observed difference vanishes when