An improved Compton parameter map of thermal Sunyaev-Zeldovich effect from Planck PR4 data

Taking advantage of the reduced levels of noise and systematics in the data of the latest Planck release (PR4, also known as NPIPE), we construct a new all-sky Compton-$y$ parameter map (hereafter, $y$-map) of the thermal Sunyaev-Zeldovich (SZ) effect from the Planck PR4 data. A tailored Needlet Internal Linear Combination (NILC) pipeline, first validated on detailed sky simulations, is applied to the nine single-frequency Planck PR4 sky maps, ranging from $30$ to $857$ GHz, to produce the PR4 $y$-map over 98% of the sky. Using map comparisons, angular power spectra and one-point statistics we show that the PR4 NILC $y$-map is of improved quality compared to that of the previous PR2 release. The new $y$-map shows reduced levels of large-scale striations associated with $1/f$ noise in the scan direction. Regions near the Galactic plane also show lower residual contamination by Galactic thermal dust emission. At small angular scales, the residual contamination by thermal noise and cosmic infrared background (CIB) emission is found to be reduced by around 7% and 34%, respectively, in the PR4 $y$-map. The PR4 NILC $y$-map is made publicly available for astrophysical and cosmological analyses of the thermal SZ effect.


INTRODUCTION
The cosmic microwave background (CMB) radiation undergoes spectral and spatial distortions while travelling from the last-scattering surface at the recombination epoch up to our instruments at the present time because of the scattering and the deflection of the CMB photons by the matter intervening along the line of sight (los).Such distortions induce secondary CMB temperature anisotropies (Aghanim, Majumdar & Silk 2008).Analysing them allows us to use the CMB as a backlight to probe the baryonic and dark matter distributions in the Universe (Basu et al. 2021).
The most prominent spectral distortion of the CMB anisotropies arises from the thermal Sunyaev-Zeldovich (SZ) effect (Zeldovich & Sunyaev 1969;Sunyaev & Zeldovich 1972): when CMB photons travel through a hot ionised gas of electrons, predominantly located in the potential wells of massive galaxy clusters, they get upscattered to higher energies by the electrons through an inverse Compton scattering process.This causes an overall shift of the CMB blackbody spectrum to higher frequencies as the total number of photons is conserved, thus leading to a characteristic spectral signature of the thermal SZ effect, with a decrement of CMB intensity at low frequency (< 217 GHz) in the direction of galaxy clusters and an increment of CMB intensity at high frequency (> 217 GHz).
The peculiar frequency dependence of the thermal SZ effect has allowed the detection of thousands of galaxy clusters from multifrequency observations of the microwave sky over the past decade ★ E-mail: chandran@ifca.unican.es† E-mail: remazeilles@ifca.unican.es‡ E-mail: barreiro@ifca.unican.es(Bleem et al. 2015;Planck Collaboration XXVII 2016;Hilton et al. 2021;Melin et al. 2021), but also the mapping of the thermal SZ Compton- parameter (hereafter, -map) from the entire hot gas all over the sky, which includes diffuse, unbound gas between clusters (Planck Collaboration XXII 2016;Aghanim et al. 2019;Madhavacheril et al. 2020;Tanimura et al. 2022;Bleem et al. 2022).With increasing sensitivity and resolution, next-generation CMB experiments are expected to release even larger cluster catalogues (Abazajian et al. 2019;Ade et al. 2019) and cleaner -maps of the hot gas (Hanany et al. 2019;LiteBIRD Collaboration et al. 2022) in the near future.
Being independent of the redshift, the thermal SZ effect serves as an important cosmological probe of the large-scale structure in the Universe (Birkinshaw 1999;Carlstrom, Holder & Reese 2002).Cluster number counts as a function of the redshift from current SZ catalogues provide cosmological constraints on the amplitude of dark matter fluctuations,  8 , the matter density, Ω m , and the dark energy equation-of-state parameter, , which are independent of the constraints from primary CMB anisotropies, exhibiting the first tensions with respect to ΛCDM model predictions from the high-redshift CMB probe (Planck Collaboration XX 2014;Planck Collaboration XXIV 2016).Unlike cluster catalogues which solely rely on the most massive clusters that can be detected individually, Compton -maps probe the full thermal SZ emission over the sky, including the fainter emission from low-mass clusters and the diffuse, unbound gas outside clusters which in fact contribute statistically to the signal.As such, Compton -maps provide another important and complementary cosmological probe through the angular power spectrum of the Compton- field (Komatsu & Kitayama 1999;Refregier et al. 2000;Komatsu & Seljak 2002;Planck Collaboration XXII 2016;Bolliet et al. 2018;Remazeilles et al. 2019;Rotti et al. 2021;Tanimura et al. 2022), higher-order statistics (Rubiño-Martín & Sunyaev 2003;Bhattacharya et al. 2012;Wilson et al. 2012;Hill & Sherwin 2013;Planck Collaboration XXII 2016;Remazeilles et al. 2019) and cross-correlations with other tracers of the large-scale structure (e.g.Hill & Spergel 2014).
However, extracting thermal SZ Compton- anisotropies out of microwave sky observations is challenging because the signal is faint compared to Galactic and extragalactic foreground emissions at submillimetre wavelengths.In addition, thermal noise and instrumental systematics add further contamination to the data.The most significant foreground to thermal SZ emission at small angular scales arises from cosmic infrared background (CIB) anisotropies due to the cumulated emission of dusty star-forming galaxies.At the current stage where the model of various foregrounds like the CIB is relatively poorly known, the use of blind (i.e.non-parametric) component separation methods is warranted for thermal SZ map reconstruction.Hence, the latest all-sky thermal SZ -maps which have been publicly released by the Planck Collaboration (Planck Collaboration XXII 2016) have been obtained using tailored versions of the blind Internal Linear Combination (ILC) method for the reconstruction of the thermal SZ effect (Remazeilles, Delabrouille & Cardoso 2011a;Hurier, Macías-Pérez & Hildebrandt 2013;Remazeilles, Aghanim & Douspis 2013).The two public Planck thermal SZ maps were named NILC -map and MILCA -map after the respective names of the two component separation methods that were used.Both methods are ILC techniques but employ different frameworks for localization in pixel and spherical harmonic domains.For technical details, we refer the reader to Planck Collaboration XXII (2016) and the references therein.
These latest public all-sky thermal SZ -maps date back from 2015 as a product of the second Planck PR2 data release (Planck Collaboration XXII 2016).However, the Planck mission had four data releases in total, and the latest so-called PR4 data release in 2020 had significant updates with reduced noise and better control of systematics and calibration thanks to the NPIPE processing pipeline (Planck Collaboration Int. LVII 2020).
In this paper, we reconstruct an updated and improved all-sky thermal SZ Compton -map over 98% of the sky from the Planck Release 4 (PR4) data using a Needlet Internal Linear Combination (NILC, Delabrouille et al. 2009) specifically tailored for thermal SZ component separation (Remazeilles, Delabrouille & Cardoso 2011a;Remazeilles, Aghanim & Douspis 2013).A similar update has recently been reported for the MILCA -map (Tanimura et al. 2022), but the -map is not public to our knowledge.Benefiting from the improved quality of the latest Planck PR4 data, our new PR4 NILC -map is made public to the community for astrophysical and cosmological SZ analyses and cross-correlation studies.This paper is organised as follows.In Section 2, we introduce the Planck PR4 data used to construct the new -map, as well as some external data sets used as foreground tracers to characterise residual contamination in the -map.In Section 3, we describe our implementation of the NILC component separation method for thermal SZ reconstruction, highlighting the differences of processing with respect to the PR2 analysis.We present our results in Section 4 with a visual inspection of the PR4 NILC -map and comparison with PR2 -maps, estimation of its angular power spectrum and one-point statistics.In Section 5, we estimate the levels of residual contamination due to foregrounds and noise in the PR4 NILC -map and compare these with those of the PR2 NILC -map.We present our conclusions in Section 6.

Planck PR4 data
The latest PR4 (NPIPE) data release from Planck,1 as described in Planck Collaboration Int. LVII (2020), is used in this work for thermal SZ Compton -map reconstruction.The NPIPE processing pipeline was used to streamline the conversion of both LFI (Low-Frequency Instrument) and HFI (High-Frequency Instrument) raw time-ordered data (TOD) into nine calibrated full-sky maps corresponding to the nine frequency channels of Planck.
The main differences of the PR4 data with respect to earlier PR2 data that could benefit the reconstructed thermal SZ -map are: (i) Reduced noise levels in the PR4 frequency maps due to adding 8% more data from the repointing manoeuvre.
(ii) A different fitting of 4 K lines, a better flagging of pixels, and a smoother glitch removal which also contributes to the decrease of noise and half-ring correlations.
(iii) A destriping of data done with Madam (Keihänen, Kurki-Suonio & Poutanen 2005) using extremely short baselines which reduces the stripes due to systematic effects in the scanning direction of the Planck satellite in the PR4 sky maps.
(iv) Differences in the frequency bandpass responses for the PR4 HFI channels due to the differences in calibration between the PR2 processing pipeline and the PR4 NPIPE processing pipeline.
(v) Calibration of LFI and HFI data performed in a coherent pipeline.
We use the nine single-frequency full-mission maps from PR4, ranging from 30 to 857 GHz, to reconstruct the thermal SZ -map.We also use the two half-ring (HR) data splits from PR4 in nine frequency channels, which correspond to the first and second half of each stable pointing period of Planck and thus have practically uncorrelated noise.The two data sets from each half-ring are called HR1 maps and HR2 maps from here on.The principal use of this data split is to characterise the statistics of the noise in the PR4 fullmission -map but also to produce additional HR1 and HR2 -maps with maximally uncorrelated noise for thermal SZ power spectrum estimation.
All resultant -maps are given in the HEALPix pixelation scheme2 (Górski et al. 2005) with a pixel resolution of  side = 2048.The input PR4 sky maps are of resolution  side = 1024 for the LFI frequency channels (30-70 GHz) and  side = 2048 for the HFI frequency channels (100-857 GHz).For each frequency channel map, there is a smoothing effect due to the finite resolution of the optical beam of the detectors.This is treated using an effective symmetric beam transfer function for each channel.The specific instrumental beam windows from PR4 are used for component separation instead of the approximate Gaussian beam windows used in the PR2 -map analysis (Planck Collaboration XXII 2016).The PR4 beam window at 353 GHz, for instance, deviates from the Gaussian approximation by approximately 2% on average across the multipole range ℓ = 1000-2048.This discrepancy thus occurs at small angular scales where the thermal SZ signal from galaxy clusters prevails.Beam modelling errors are similar to calibration errors in an ILC, for which previous studies have demonstrated that even a minor percentage error can degrade the signal reconstruction in the high signal-tonoise regimes (Dick, Remazeilles & Delabrouille 2010).Therefore, accurate beam deconvolution at ℓ > 1000 using PR4 instrumental beams, instead of relying on Gaussian approximations, is preferred for the reconstruction of the small-scale features in the thermal SZ signal from compact galaxy clusters.

Masks
Although the component separation process by NILC is fairly localised on the pixelated sphere by construction, it is not perfectly local, so that the few pixels with the strongest emission in the Galactic centre can create ringing effects when we perform spherical harmonic transforms during needlet decomposition.These ringing effects can result in an overestimation of the sky-RMS signal at higher Galactic latitudes, ultimately affecting the effectiveness of foreground cleaning in those areas.To prevent unwanted ringing effects, we have followed the same strategy used in the Planck PR2 analysis (Planck Collaboration XXII 2016), by masking only the brightest 2% of pixels at 857 GHz along the Galactic ridge in all PR4 frequency maps before passing them through the NILC pipeline.This small processing mask, called NILC-MASK hereafter, is shown in Fig. 1 as the white area.The resultant PR4 -map is thus delivered over a fraction  sky = 98% of the sky.
The statistical analysis of the PR4 -map, including estimation of the SZ power spectrum and one-point probability density function (1-PDF) of the -map, requires masking the brightest extragalactic sources and a larger portion of the Galactic region in the -map to mitigate the residual foreground contamination after component separation.We use the apodized Galactic mask released from the PR2 analysis (Planck Collaboration XXII 2016), hereafter called GAL-MASK, conserving about  sky = 60% of the -map for statistical analysis.
For masking extragalactic radio sources in the PR4 -map (see Section 5.3), we use Planck point-source masks specifically constructed for the PR4 data at each frequency channel using the Mexican Hat Wavelet 2 (López-Caniego et al. 2006;Planck Collaboration XXVI 2016) as part of the Sevem pipeline (Planck Collaboration Int.LVII 2020).For those frequency channels with higher resolution than 10 ′ (i.e.⩾ 100 GHz), the masks are further convolved with a Gaussian beam of 10 ′ and made binary again by setting a threshold of 0.75.This is because our PR4 -map, like the public PR2 -maps, has a resolution of 10 ′ , so we need to increase the hole size for the sources to match this resolution.To study the effect of point source residuals in the -map, three different combinations of point-source masks, including only LFI, 30-143 GHz or all frequency channels are considered in this work (see section 5.3).Our reference point-source mask corresponds to that constructed from the masks of channels 30 to 143 GHz, which we call PS-MASK.For power spectra computation, the PS-MASK has been apodized with 0.1 deg transition length using the C1 apodization scheme in NaMaster (Alonso et al. 2019).
The combined GAL-MASK and PS-MASK for statistical analysis is displayed in black in Fig. 1, retaining a total sky fraction of 56%.

Foreground tracers
Estimating the residual contamination left by Galactic and extragalactic foregrounds after component separation is an essential part of determining the quality of the thermal SZ -map.This can be done at the map level, by visual comparison of the -map and a foreground template in specific regions of the sky, or by computing the crosspower spectrum between the foreground template and the -map as long as the foreground template does not suffer from thermal SZ contamination.
The two major foreground contaminants in thermal SZ maps are the CIB, a diffuse extragalactic dust emission from early star-forming galaxies yielding significant power at small angular scales, and the thermal dust emission from our Galaxy which prevails at large angular scales (see Fig. A2 in Appendix A).
To assess residual Galactic dust contamination in the -maps, we use as dust template the Improved Reprocessing of the IRAS Survey (IRIS) 100-m map (Neugebauer et al. 1984;Miville-Deschênes & Lagache 2005).At a wavelength of 100 m (equivalently, a frequency of ∼ 3000 GHz), the thermal SZ effect is completely negligible in the IRIS 100-m map and the sky emission is dominated by Galactic thermal dust emission, making the IRIS 100 m map a reliable tracer of large-scale dust contamination in the -maps.
To assess residual CIB contamination at small angular scales in the -maps, we use two independent Planck-based templates of the CIB emission, both at 857 GHz because the Planck 857 GHz channel map is not used for the construction of both PR2 and PR4 -maps at multipoles ℓ > 300 (see Planck Collaboration XXII 2016, and Section 3), which prevents from unwanted noise correlations between the -map and the CIB template.The intensity of the thermal SZ emission is also mostly insignificant at 857 GHz in the CIB templates.As a first template, we use the Planck GNILC CIB map at 857 GHz (Planck Collaboration Int.XLVIII 2016), which was processed with the datadriven Generalized Neelet ILC (GNILC) method (Remazeilles, Delabrouille & Cardoso 2011b) to disentangle CIB from thermal dust emission.As a second, independent template, we use the CIB map at 857 GHz from Lenz, Doré & Lagache (2019), derived from Planck data using a model-dependent approach to subtract thermal dust contamination based on Hi gas column density.

METHODOLOGY
We apply mostly the same NILC algorithm to Planck PR4 data as the one used for the PR2 -map release in Planck Collaboration XXII (2016), with some nuances as the data sets have different characteristics.The major steps and specifications of the current implementation of NILC on the Planck PR4 data are described hereafter.

Signal modeling
The thermal SZ (tSZ) effect is a frequency-dependent anisotropic distortion of the CMB temperature resulting from inverse Compton scattering of CMB photons off a hot gas of free electrons (Zeldovich & Sunyaev 1969): Here,  () is the characteristic frequency-dependence of the distortion (Fig. 2).In the non-relativistic limit and in thermodynamic temperature units, it is given by the analytic form: where  is the dimensionless frequency defined as Here, ℎ is the Planck constant,  is the Boltzmann constant, and  CMB is the CMB blackbody temperature.The direction-dependent Compton parameter ( n) in equation (1) represents the amplitude of the distortion, which is proportional to the los integral of the electron gas pressure  e =  e  e through the Thomson scattering cross-section: (4) In the equation above,  e is the electron gas temperature,  e is the electron number density,  e  2 is the electron rest mass energy, and   is the Thomson scattering cross-section.
To get an accurate spectral response of the thermal SZ effect in Planck frequency bands, the frequency dependence (), hereafter spectral energy distribution (SED), must be integrated over Planck PR4 frequency bandpasses.The resulting thermal SZ SED coefficients across frequencies for PR4 are shown in Fig. 2 (yellow dots) and listed in Table 1 in thermodynamic temperature units.They slightly differ from the coefficients of the PR2 analysis in HFI channels (Planck Collaboration XXII 2016) due to slightly different HFI bandpasses from the PR4 data release.
Mirroring PR2 assumptions in Planck Collaboration XXII (2016) for the sake of comparison, we have neglected relativistic corrections to the thermal SZ SED (Challinor & Lasenby 1998;Itoh, Kohyama & Nozawa 1998) in the PR4 analysis.However, it has been shown by Remazeilles et al. (2019) that relativistic SZ corrections statistically have an impact on the amplitude of the measured Planck -map power spectrum and the skewness of the 1-PDF.Future work will consider the release of another -map from PR4 that accounts for the relativistic SZ effect, which arises due to the variable temperature of the electron gas across the sky.
The observed data (, ) across the frequencies  can thus be expressed for all pixels  like: where (, ) is the unparametrized nuisance term accounting for all possible foreground emissions and the instrumental noise.The approach is therefore blind to the foregrounds because we do not assume a specific spectral model for the foregrounds, whose exact spectral properties are much less known than for the thermal SZ signal.Equation ( 5) can be inverted using the blind NILC method described hereafter to recover the Compton- parameter in each pixel from the multi-frequency data, with minimum-variance residual contamination from foregrounds and noise.
Hence, here we outline the main ingredients of the method and the specificities of our implementation for PR4, highlighting where it differs from the PR2 NILC implementation.
The Planck PR4 frequency maps are provided with the dipole and the frequency-dependent dipole-induced quadrupole included (Planck Collaboration Int.LVII 2020).Therefore, we first subtracted the dipole and the frequency-dependent quadrupole from the PR4 frequency maps using the best-fit templates from Commander (Eriksen et al. 2008) that are available at NERSC. 3 The PR4 maps are also masked with the small NILC-MASK (Fig. 1) to discard the 2% most-contaminated region in the Galactic plane from the NILC analysis, so that the released PR4 -map effectively covers  sky = 98% of the sky.
In order to produce the PR4 -map at the same 10 ′ angular resolution as the public PR2 -maps, the PR4 frequency maps have to be deconvolved from their native instrumental beam and reconvolved with a common 10 ′ symmetric Gaussian beam.In contrast to the PR2 SZ analysis where approximate Gaussian beams were used for beam deconvolution (Planck Collaboration XXII 2016), here we use the specific PR4 beam transfer functions4 of each frequency channel to ensure an accurate reconstruction of the small-scale thermal SZ emission from compact clusters.This operation is done in harmonic space and summarised by the following scheme: where SHT stands for spherical harmonic transform,  PR4 is the 10 ′ Gaussian beam transfer function.
Instead of operating on the maps (, ) in pixel space (pixelbased ILC) or on the spherical harmonic coefficients   ℓ in harmonic-space (harmonic ILC), NILC operates on needlet coefficients.Needlets are a set of wavelets forming a tight frame on the sphere which provide simultaneous localization in pixel domain and in harmonic space (Narcowich, Petrushev & Ward 2006;Marinucci et al., 2008).Localization in the pixel domain is important for the reconstruction of spatially localized signals such as thermal SZ and for adjusting foreground cleaning depending on Galactic latitudes, as Galactic foreground contamination is mostly localized at low latitudes while CIB and noise dominate at high latitudes.Simultaneous localization in harmonic space also enables differentiating between Galactic foregrounds dominating at large angular scales and CIB and noise dominating at small angular scales for customised foreground cleaning.In Appendix B, we have implemented a harmonic ILC (HILC) on PR4 data to highlight how NILC clearly outperforms HILC on thermal SZ map reconstruction.
To perform the decomposition of the PR4 data into needlet coefficients, we define 10 Gaussian-shaped needlet window functions ℎ  ℓ in harmonic space as where are the SHT of Gaussian functions with full-width-at-half-maximum 600, 300, 120, 60, 30, 15, 10, 7.5, 5] in arcmin for  = 1, . . ., 9. The needlet window functions defined in equation ( 7) thus satisfy the condition which guarantees the conservation of the signal at all angular scales after forward and inverse needlet transformations.As shown in Fig. 3, the 10 needlet windows ℎ  ℓ operate as bandpasses in harmonic space, each selecting a range of angular scales to ensure localization in harmonic space for component separation.For each needlet scale , the needlet coefficients  , ( ) are computed from the spherical harmonic coefficients   ℓ of the PR4 maps (equation 6) bandpass-filtered with the needlet window ℎ  ℓ : where  ℓ ( ) are spherical harmonics.For each frequency channel , we thus obtain 10 needlet maps  , ( ), each of them displaying sky emission for a specific range of angular scales as selected by the needlet windows.
An estimate of the Compton- parameter map at needlet scale  is obtained from a weighted linear combination of the PR4 needlet maps across the frequency channels, using the NILC weights: These weights depend only on   , which are the SED coefficients of the thermal SZ effect across the frequency channels (Table 1), and  −1  ( ), which is the inverse of the empirical covariance matrix of the PR4 data at pixel  and needlet scale ,   ( ), whose elements for all pairs of frequency channels are estimated as The Gaussian convolution kernel   ( ,  ′ ) defines, for each needlet scale , the effective size of the pixel domain around pixel  over which the product of data      ′  for a pair of frequency channels is averaged.For each needlet scale , the FWHM of   ( ,  ′ ) is chosen small enough for localized estimates of the sky covariance in pixel space, but large enough to average as many pixels as possible to minimise empirical chance correlations between signal and contaminants and keep the so-called ILC bias (Delabrouille et al. 2009) under control.
By construction, the NILC weights (equation 12) give unit response to the thermal SZ component since       = 1, such that the NILC estimate   ( ) (equation 11) recovers the full thermal SZ signal  at needlet scale  without multiplicative error.There is only an additive error to the -estimate due to residual foreground and noise contamination, but this error is kept as small as possible since the NILC weights give, by construction, the minimum-variance solution at the needlet scale .
An important aspect of the analysis is that different subsets of frequency channels are selected depending on the needlet window  for the construction of the NILC weights (equation 12) and -map estimate (equation 11) at that needlet scale (see Table 2).All nine LFI and HFI channels (30-857 GHz) are used by NILC in the first three needlet windows for large angular scales.However, due to their low resolution and limited sensitivity at smaller angular scales, LFI channel maps are excluded from subsequent needlet windows, as they are not well-suited for probing sky emission at these smaller angular scales.Therefore, only the 6 HFI channels (100-857 GHz) are combined by NILC in the fourth up to the sixth needlet window.Following Planck Collaboration XXII (2016), the 857 GHz channel map is not used at multipoles ℓ > 300 to mitigate CIB and infrared source contamination, which can be substantial in this channel at small angular scales.Consequently, only 5 HFI channels (100-545 GHz) are combined by NILC in the seventh up to the last needlet band.
The reconstructed needlet -maps,   ( ) (equation 11), from each needlet scale  are finally transformed into spherical harmonic coefficients,  ℓ,  , and combined together to form the complete PR4 NILC -map, ( ): The same NILC weights, computed from the full PR4 data in Eqs. ( 12)-( 13), are also applied to the two sets of half-ring PR4 maps (HR1 and HR2 frequency maps), which went through the same needlet decomposition process.The resulting PR4 HR1 and HR2 -maps,  HR1 ( ) and  HR2 ( ), have same thermal SZ signal and residual foreground contamination but mostly-uncorrelated noise.The half-difference between the PR4 HR1 and HR2 -maps,  HR1 ( ) −  HR2 ( ) /2, cancels out any sky emission but not the instrumental noise, and thus serves as a noise map estimate whose statistical properties are the same as those of the actual residual noise in the full PR4 -map.In addition, because of least-correlation of noise between the two half-rings, the cross-power spectrum between the PR4 HR1 -map and the PR4 HR2 -map enables to estimate the recovered thermal SZ angular power spectrum corrected for the instrumental noise bias.

PRTHERMAL SZ 𝑦-MAP CHARACTERIZATION
In this section, we assess the quality of the PR4 NILC -map and compare it with that of the PR2 -maps in terms of noise, residual systematics and foreground contamination, through map inspection, one-point statistics and power spectrum analysis.

Maps inspection
Fig. 4 shows, in orthographic projection, the new PR4 NILC -map (top panel) produced in this work, along with the PR2 NILC -map (bottom panel) released by the Planck Collaboration in 2015 (Planck Collaboration XXII 2016) for comparison.The left-hand side shows the northern hemisphere and the right-hand side shows the southern hemisphere with respect to Galactic coordinates.Both -maps are at the same 10 ′ angular resolution and show relatively high consistency at the map level.However, some differences between the two maps are already visible.The PR2 -map shows more prominent noise (small-scale granular pattern) than the PR4 -map in the bottom part of the northern hemisphere.In addition, blue patchy patterns are visible in the bottom part of the southern hemisphere of the PR2 -map, while these residuals are absent from the new PR4 -map.Thermal SZ sources are visible as red spots in the given colour scale.Prominent galaxy clusters like Coma and Virgo are clearly visible near the north pole.Some residual infrared compact sources may also appear in red in the maps while residual radio sources appear in blue due to the sign of the NILC weights flipping from positive at high frequencies to negative at low frequencies to match the spectral response of the thermal SZ sources.Diffuse Galactic foreground contamination, dominated by thermal dust, is visible with significant power at both ends of the range near the Galactic plane (along the edges in Fig. 4).
The black outline in Fig. 4 traces the boundary of the GAL-MASK masking 40% of the sky that, together with those pixels contaminated by point sources, are excluded in the 1-PDF and power spectrum analysis (Sections 4.2 and 4.3).The regions around the Galactic plane are best to be excluded from statistical analysis due to significant power from residual Galactic emission.However, these regions are still clean enough to spot thermal SZ sources, hence the release of the -map over 98% of the sky.
There are large-scale stripes along the scan direction of the satellite in the Planck -maps due to residual systematics in the Planck map-making pipeline (see Section 2.1).Different methods used for baseline noise correction and destriping affect the morphology of the residual stripes in the -maps (see Planck Collaboration XXII 2016, section 4.1).Improved destriping in Planck PR4 data compared to Planck PR2 data is clearly visible in Fig. 5, which shows the PR4 NILC -map (top) versus the PR2 NILC -map (middle) and the PR2 MILCA -map (bottom) after bandpass-filtering in multipole range.All the maps in Fig. 5 have been filtered with a common, analytical bandpass filter which selects the range of multipoles between ℓ = 20 and ℓ = 500, where the thermal SZ signal dominates.
Both the PR2 MILCA and PR2 NILC -maps exhibit more residual striping than the PR4 NILC -map.The NILC -maps also show lower levels of residual Galactic foreground contamination around the Galactic plane compared to the PR2 MILCA -map.Although the MILCA -map was updated using PR4 data in Tanimura et al. (2022), it is not publicly available, hence not included in this comparison.While stripes are large-scale residual patterns, they cause the clusters and granular noise in their direction to appear brighter (negative or positive) than the rest, hence the importance of improved destriping in the PR4 -map.Fig. 6 further highlights the reduced level of destriping at intermediate Galactic latitude in the PR4 NILC -map compared to the PR2 NILC -map.
Galactic thermal dust emission is the dominant foreground contaminant of the -map at large angular scales and low Galactic lat-    -8) and CIB (see Section 5.2) more efficiently in regions and angular scales where they are dominant.This can make thermal SZ sources in the background more visible and identifiable.
Fig. 8 shows galaxy clusters in regions around the Galactic plane, for which the PR4 -map (left-hand panels) shows lower diffuse foreground contamination than the PR2 -map (right-hand panels).The identified clusters shown in the centre of these images have been observed by the Clusters in the Zone of Avoidance (CIZA) project, which is an X-ray survey for galaxy clusters hidden by the Milky Way (Ebeling et al. 2002;Kocevski et al. 2007).Their CIZA names are given in the caption of Fig. 8.These regions should be excluded for cosmological inference as the contamination from Galactic foregrounds is still strong enough to cause biases on the thermal SZ power spectrum.However, lower residuals from the diffuse and filamentary structure of the Galactic foreground emission in the PR4 -map could in principle help in identifying more thermal SZ sources in regions near the Galactic plane and in defining their morphology with better accuracy.

1-PDF analysis
When CMB photons pass through hot gas in galaxy clusters, they are more likely to gain energy by inverse Compton scattering than to lose it.Therefore, the Compton- parameter field has a highly non-Gaussian distribution with a positive skewness (Rubiño-Martín & Sunyaev 2003).Being mostly insensitive to Gaussian contaminants, the characteristic skewness of the thermal SZ emission can also be used, as an advantageous alternative to the power spectrum, to constrain the cosmological parameter  The one-point probability distribution function (1-PDF) of the thermal SZ Compton- field is computed from the histogram of the -maps over 56% of the sky after masking the Galactic plane with the GAL-MASK and the point-sources with the PS-MASK.The resulting 1-PDF of the PR2 (solid red) and PR4 (solid black) maps, normalised to unity at maximum, are plotted in Fig. 9, showing the positive, skewed tail of the distribution which is characteristic of thermal SZ emission.The positive tails of the PR2 and PR4 NILC -maps match together, ensuring consistent recovery of the thermal SZ signal.However, the widths of the distributions do not match, indicating a reduced variance by 17% due to the lower foreground and noise contamination in the PR4 NILC -map compared to the PR2 NILC -map.Reduced variance owing to lower noise in the PR4  NILC -map is also evident from the tighter 1-PDF of the noise (see Fig. 11 and Section 5.1), which was computed for PR2 and PR4 from the histogram of the half-difference of HR1 and HR2 -maps.The blue dashed line in Fig. 9 shows the 1-PDF of the -map reconstructed from PR4 data using a harmonic-domain ILC (HILC), without localization in the pixel domain.The HILC -map is clearly inferior to the NILC -map with a much larger variance of the dis-tribution.This highlights the importance of spatial localization to reconstruct thermal SZ signals, as NILC proceeds, while HILC is not a recommended component separation method for thermal SZ mapping.This aspect is further elaborated in Appendix B As can be seen from simulations in Fig. A1 (Appendix A), which shows the contributions of various residual foregrounds to the 1-PDF of the NILC -map, non-Gaussian foregrounds from Galactic emission and extragalactic infrared sources add reasonably low skewness to the positive tail of the -map 1-PDF, while extragalactic radio sources add significant negative skewness to the distribution if not masked.Hence, we chose to use point-source masks from frequency channels below 217 GHz in order to suppress the negative contribution from radio sources in the PR4 and PR2 -map 1-PDFs, while we kept infrared sources since masking them can cause some slight loss of power in the thermal SZ tail of the distribution.This choice of masking is justified as the infrared sources, which dominate at high frequencies, are not a major contaminant to the thermal SZ -map (Fig. A1).Furthermore, it is possible that a small number of these sources are actually unresolved SZ sources that have been mistakenly identified as infrared sources in catalogues.Point-source residuals and masking are further explored in Section 5.3.
To get a constraint on  8 , we followed the approach presented by Wilson et al. (2012) and computed the unnormalized skewness ⟨  3 ( )⟩ of the PR4 NILC -map using the same sky area and source mask as in Planck Collaboration XXII (2016).To get a rough estimate of the uncertainty, we also computed the skewness of the associated noise estimated from the half-difference of half-ring -maps.By applying the characteristic scaling relation of ⟨  3 ( )⟩ ∼  11 8 , we derived a value of  8 = 0.76 ± 0.02 for the PR4 NILC -map, which is consistent with the value reported for the PR2 NILC y-map ( 8 = 0.77 ± 0.02) in Planck Collaboration XXII (2016).It was anticipated, given the matching positive tails of the 1-PDF of the PR4 and PR2 NILC -maps (Fig. 9).

Angular power spectrum analysis
The angular power spectrum of the thermal SZ effect has long been recognised as an important astrophysical and cosmological probe (Komatsu & Kitayama 1999  2002), as it integrates all the thermal SZ emission in the sky, both from diffuse, unbound hot gas and compact clusters of any mass and redshift.All-sky -maps from the Planck survey allow for computing the thermal SZ angular power spectrum over a large fraction of the sky and a relatively broad range of angular scales.
We use NaMaster5 (Alonso et al. 2019) to compute the angular power spectra of the -maps.It allows us to reconstruct the thermal SZ power spectrum from the masked -maps using a pseudo- ℓ estimation, while also taking care of the mode-coupling due to the application of the mask, the beam convolution, the pixelization, and the multipole binning.The PR2 and PR4 -map power spectra are computed with the GAL-MASK and PS-MASK, leaving about 56% of the sky available after apodization of the PS-MASK.A custom binning scheme is defined in the plots with linear bins of width Δℓ = 3 from multipoles ℓ = 2 to 30 and logarithmic bins with Δ log(ℓ) = 0.05 from ℓ = 30 onwards.This binning gives us 95 band powers.All power spectra are computed up to ℓ = 2048 since the data at higher multipoles is consistent with noise due to the 10 ′ resolution of the -maps.
Once the point sources are masked, instrumental noise dominates the power at high multipoles.To correct for the noise bias in the estimated thermal SZ power spectrum, we compute the cross-power spectrum between the PR4 HR1 and HR2 -maps since these two maps have mostly-uncorrelated noise due to half-ring data split (Section 3.2): where  HR1 ℓ and  HR2 ℓ are the spherical harmonic coefficients of the PR4 HR1 and HR2 -maps.The same is done for PR2 using public half-ring -maps from Planck Collaboration XXII (2016).
Fig. 10 shows the auto-power spectra of the PR2 and PR4 NILC maps (dashed lines) along with the reconstructed thermal SZ power spectra obtained from the cross-spectra between HR1 and HR2 maps (solid lines) for both PR2 (red) and PR4 (black).As a reference, the thermal SZ power spectrum from the Planck Sky Model (PSM; Delabrouille et al. 2013), which is used for the NILC analysis on Planck simulations in Appendix A, is overplotted as a thin blue line in Fig. 10.
The power spectrum of the noise in the -maps is also estimated using the half-difference of the half-ring -maps, (HR1 -map − HR2 -map)/2, and plotted in Fig. 12 for PR2 (red) and PR4 (black).The PR4 NILC -map benefits from lower noise at all angular scales compared to the public PR2 NILC -map, as a consequence of the reduced levels of noise in the Planck PR4 data from the inclusion of 8% more data (see Section 5.1 for further discussion).As we will see, the reduced levels of noise in the PR4 data give the possibility to NILC to further minimize the variance of extragalactic foreground contamination (CIB and radio sources) at small angular scales.
As we can see in Fig. 10, instrumental noise dominates the power at small angular scales (dashed lines), thus biasing the recovered thermal SZ power spectrum.After correcting for the noise bias using cross-power spectrum between half-ring -maps (solid lines), there is still some remaining excess power at high multipoles in the reconstructed thermal SZ power spectrum which is associated mostly with residual CIB contamination in the -maps, as also confirmed by the analysis on Planck simulations in Fig. A2 (see Appendix A).In contrast, the excess power at low multipoles is due to residual Galactic foreground contamination at large angular scales.Clearly, the PR4 map power spectrum (solid black) shows much less power at low and high multipoles compared to the PR2 -map power spectrum (solid red), which indicates lower residual foreground contamination in the PR4 NILC -map (see Section 5.2 for further discussion).The PR2 and PR4 -map power spectra are more consistent at intermediate multipoles (ℓ ∼ 30-300) where the reconstructed thermal SZ signal dominates over the residual Galactic and extragalactic foreground contamination.Around the same range of multipoles (ℓ ∼ 50-300), the reconstructed signal is also close to the PSM model of thermal SZ emission.While the PR4 -map conclusively has significantly lower contamination from extragalactic foregrounds at small angular scales compared to the PR2 -map, the difference in power seen at low multipoles is not statistically significant because of the large cosmic variance expected from the non-Gaussian SZ signal (Cooray 2001;Komatsu & Seljak 2002;Bolliet et al. 2018).Further analysis of the residual contamination in the thermal SZ -maps is done in Section 5.

RESIDUAL FOREGROUND AND NOISE CONTAMINATION IN THE PR4 𝑌 -MAP
According to Fig. A2 in Appendix A, where the same NILC pipeline is applied to Planck simulations, the major contaminants to the reconstructed -map are the instrumental noise, the CIB and extragalactic radio sources at small angular scales and the Galactic foreground emission at large angular scales.In contrast, CMB and infrared sources are not major contaminants to the NILC -map.
While residual compact sources and Galactic foregrounds can be further mitigated by the application of appropriate masks, this is not possible for CIB whose emission is more diffuse and homogeneous over the sky.In this section, we compare and quantify the residual contamination of the PR2 and PR4 NILC -maps by the noise, the CIB and extragalactic compact sources.

Noise
The half-difference between the HR1 and HR2 PR4 -maps gives us a noise map whose statistical properties are those of the actual residual noise fluctuations in the full (ring) PR4 -map.Similarly, we get an estimate of the noise contamination in the PR2 -map from the half-difference of the public PR2 half-ring NILC -maps released by the Planck Collaboration (Planck Collaboration XXII 2016).The resulting -noise maps are used to compare the noise characteristics of the PR4 and PR2 NILC -maps.Fig. 11 shows the binned normalized histogram (1-PDF) of the residual noise for the PR4 NILC -map (black), the PR2 NILC map (red), and the HILC -map (blue).The best-fitting Gaussian PDF (dashed line) is also shown in each case.As we can see, the noise distribution is mostly Gaussian in all -maps.The PR4 NILC -map (best-fitting Gaussian noise standard deviation of  = 1.12 × 10 −6 ) has lower noise as compared to the public PR2 NILC -map (best-fit Gaussian noise standard deviation of  = 1.16 × 10 −6 ), with a 6.8% reduction of the noise variance.The HILC -map has comparatively much higher noise than any of the NILC -maps, with a best-fit Gaussian standard deviation of  = 2.27 × 10 −6 .Fig. 12 shows the angular power spectrum of residual noise for the PR4 NILC -map (black) and the PR2 NILC -map (red) over  sky = 56% of the sky, as well as the relative decrease of noise power in the PR4 -map with respect to the PR2 -map across the multipoles, i.e.
. A linear binning of Δℓ = 20 is used for this plot.As evident from the bottom panel of Fig. 12, the PR4 -map has consistently lower noise power than the PR2 -map at all angular scales.The mean percentage improvement of residual noise power over the multipole range ℓ = 30-2048 is 6.7%, which is consistent with the result obtained from the 1-PDF.This improvement is due to the overall lower noise level in the Planck PR4 data for the reasons listed in Section 2.1.

Cosmic infrared background (CIB)
The CIB is the most significant foreground contaminant to thermal SZ -maps at small angular scales (see e.g.Fig. A2 in Appendix A).To assess and compare the levels of residual CIB contamination in the PR4 and PR2 -maps, we compute the cross-power spectrum between these -maps and Planck CIB maps at 857 GHz.We choose 857 GHz as the frequency of the CIB templates because the 857 GHz channel of Planck is not used at high multipoles ℓ > 300 for the reconstruction of the PR4 NILC -map (see Table 2) and the PR2 NILC/MILCA -maps (see Planck Collaboration XXII 2016).This enables the exclusion of spurious correlations between the noise of the -maps and that of the 857 GHz CIB map, as the noise from different frequencies is uncorrelated, while the CIB is still highly correlated across frequencies (Planck Collaboration XXX 2014).Since the thermal SZ intensity is negligible at 857 GHz (about 3% of its maximum intensity value at 353 GHz; see Fig. 2 and Table 1), using CIB maps at 857 GHz also allows us to exclude spurious correlations that would be caused by residual SZ contamination in CIB templates of lower frequency.
We use two independent CIB templates for our analysis (see Section 2.3): the Planck GNILC CIB map at 857 GHz (Planck Collaboration Int.XLVIII 2016) covering 57% of the sky and the Planck-based CIB map at 857 GHz from Lenz et al. (2019) which covers 18% of the sky.Combining with the GAL-MASK and PS-MASK of the -maps, the CIB cross -map power spectrum is computed with NaMaster over 50% of the sky when using the Planck GNILC CIB map and over 15% of the sky when using the Lenz et al. (2019) CIB map.The results are shown in Fig. 13 for high multipoles ℓ > 600, with a linear binning of Δℓ = 20.For either CIB template, a consistent pattern emerges, showing a much stronger correlation with the PR2 -maps (red / yellow) than with the new PR4 -map (blue).The residual CIB contamination in the PR4 NILC -map is thus considerably lesser than the CIB contamination of the public PR2 -maps.This contributes to lesser the overall contamination of the -map power spectrum at high multipoles as observed in Fig. 10.
The lower sub-panels of Fig. 13 display the relative decrease of CIB contamination in the PR4 NILC -map with respect to the PR2 NILC -map over the multipoles, i.e.
. By averaging over the multipole range ℓ = 600-2048, we infer a 34.2% decrease in residual CIB power in the PR4 NILC -map compared to the PR2 NILC -map over 50% of the sky when using the Planck GNILC CIB template (top panel), and a 56.7% decrease over 15% of the sky when using the Lenz et al. (2019) CIB template (bottom panel).As a matter of fact, the lower noise variance in PR4 data compared to PR2 data allows the NILC pipeline to further minimize the variance of the foregrounds at high multipoles, such as CIB (Fig. 13) and extragalactic compact sources (Section 5.3).Since the variance of instrumental noise is much larger than the variance of CIB fluctuations in the Planck data, even a few per cent reduction of noise in the PR4 data can make a big difference in the CIB variance minimization by NILC for the -map.
CIB maps reconstructed from half-ring data splits being available for the Lenz et al. (2019) template, we obtained a CIB noise map from the half-difference of half-ring CIB maps which we cross-correlated with the -noise maps also obtained from the half-ring data split.As already anticipated and shown in the bottom panel of Fig. 13 (dashed-dotted lines), the noise in the 857 GHz CIB map is mainly uncorrelated with the noise in the -maps due to the Planck 857 GHz channel not being used at high multipoles for the construction of the -maps.Hence, the noise contribution to the cross-spectrum analysis is mostly negligible and cannot account for the difference seen in Fig. 13.
We might wonder about the possibility of a bias in the cross-spectra comparison due to the use of GNILC CIB maps from the Planck PR2 release, which is the same data used for one of the -maps.However, given that we used the 857 GHz GNILC CIB map and the Planck 857 GHz map was not used for CIB-relevant multipoles in either PR2 or PR4 -maps reduces the likelihood of such a bias.Additionally, our findings are supported by the fact that the Lenz et al. (2019) CIB map, based not on PR2 or PR4 but on Planck PR3 data, yields similar cross-spectra results.

Extragalactic compact sources
Radio and infrared (IR) compact sources are another important foreground contaminant to thermal SZ emission at small angular scales.They are usually handled by the use of point-source masks.Radio sources from active galactic nuclei (AGNs) emitting synchrotron radiation are mostly detected at low frequencies below 217 GHz, while IR sources from star-forming dusty galaxies are detected at high frequencies (Planck Collaboration XXVI 2016).Hence, by masking the point sources detected in either low or high Planck frequency channels, we can explore the residual contamination from radio and IR sources separately and determine the optimal masking strategy.
In Fig. 14, we explore the impact of radio and IR source contamination on the 1-PDF of the PR4 and PR2 -maps by masking them with three different point-source masks: (i) A point-source mask defined as the union of point-source masks from all Planck channels, which masks both IR and radio sources (red line).(ii) A point-source mask defined as the union of the point-source masks from the three Planck LFI channels (30, 44 and 70 GHz), which masks only radio sources (yellow line).
(iii) A point-source mask defined as the union of the point-source masks from Planck's low frequencies below 217 GHz, where the spectral response of the thermal SZ effect is negative.This means that all the point sources that end up as negative spots in the -map are masked (black line).
The NILC weights give unit response to the SED of the thermal SZ signal across Planck channels.As the frequency dependence of the thermal SZ effect is negative at frequencies below 217 GHz (Fig. 2), and the majority of the radio sources are detected in these frequencies, they end up as negative point sources in the reconstructed Compton -maps.The opposite is true for IR sources in general.As the distribution of extragalactic sources over the sky is also non-Gaussian, radio sources contribute as a negatively skewed tail in the 1-PDF of the reconstructed -map.In contrast, IR sources contribute to the positively skewed tail of the distribution (see Fig. A1 in Appendix A for evidence from simulations).This assertion is supported by Fig. 14, where the drop of the negative non-Gaussian tail in the 1-PDF of the PR4 -map is largely obtained by masking only radio sources detected in Planck LFI channels (yellow line versus purple line).
Masking all detected Planck sources in the -map (red line) does not significantly attenuate the negative tail of the distribution beyond what is achieved by masking only the sources detected at frequencies < 217 GHz (black line).These two cases also show only marginal improvements compared to masking only LFI radio sources (yellow line).This is evident from the variance of the 1-PDF, which is reduced by 5% when masking sources detected below 217 GHz (black line) compared to the unmasked distribution (purple line), while a further reduction of only 0.5% is achieved by masking all sources (red line).In addition, when all sources in the -map are masked (red line), there is only a negligible drop compared to the black line in the positive non-Gaussian tail characteristic of thermal SZ (skewness reduced by only 1%).This suggests that masking sources detected at frequencies ⩾ 217 GHz have negligible benefits, while there is also a risk of losing part of the non-Gaussian SZ signal due to erroneous masking of unresolved compact thermal SZ sources mistaken as IR sources.Therefore, the choice is made to mask out only sources below 217 GHz in the -map using the corresponding Planck compact source masks (black line).This approach allows for a higher sky fraction without significant contamination of the signal.
As visible from Fig. 14 before masking, the negative non-Gaussian tail in the 1-PDF is slightly larger for the PR2 -map (dashed black line) compared to the PR4 -map (purple line), which suggests lower residual contamination from radio sources in the PR4 -map.Again, the lower noise levels in the PR4 data allow the NILC pipeline to further reduce the extragalactic foreground contamination from CIB and radio sources in the -map.

CONCLUSION
A new thermal SZ Compton -map has been produced over 98% of the sky by implementing a tailored NILC pipeline on the nine Planck frequency maps ranging from 30 to 857 GHz, as provided by the Planck PR4 data release.The newly introduced -map represents a substantial improvement over the Planck PR2 -maps that were previously released by the Planck Collaboration (Planck Collaboration XXII 2016).The Planck PR4 data feature reduced levels of noise and systematics, which translates into lower levels of noise and foreground contamination in the new PR4 NILC -map, compared to the public Planck PR2 -maps.
Several tests have been conducted with map inspections, one-point and two-point statistics to validate the quality of the PR4 NILC map.These tests reveal that the noise has been reduced by about 7% in the PR4 NILC -map, while the residual contamination from CIB has decreased by more than 34% compared to the PR2 -maps.Moreover, the PR4 NILC -map exhibits lower levels of large-scale striping from residual 1/  noise compared to the public Planck PR2 -maps and reduced contamination from extragalactic radio sources.The constraint on the cosmological parameter  8 = 0.76 ± 0.02, obtained from the skewness analysis of the PR4 NILC -map, remains consistent with the constraint derived from the PR2 NILC -map analysis as reported in Planck Collaboration XXII (2016).
The Planck PR4 NILC -map, as well as the associated half-ring -maps and masks, are publicly available at https://doi.org/10.5281/zenodo.7940376and in the Planck Legacy Archive.
We are considering several extensions to the current analysis for future studies.One approach is to incorporate external full-sky data in the NILC pipeline, along with the Planck PR4 channel maps, to enhance foreground cleaning in the Compton -map.This idea, similar to that of Kusiak et al. (2023) for CMB, uses additional channels from external data that trace foreground emission.
Another extension is to release a "CIB-free" PR4 -map for cross-correlation studies with large-scale structure (LSS) tracers such as lensing maps.To achieve this, constrained ILC methods (Remazeilles, Delabrouille & Cardoso 2011a;Remazeilles, Rotti & Chluba 2021) can be used to deproject the spectral moments of the CIB (Chluba, Hill & Abitbol 2017) from the PR4 data.Such an approach may significantly reduce biases caused by residual CIB-LSS correlations in thermal SZ-LSS cross-correlation studies.Recently, deprojection of CIB moments was applied to cluster detection algorithms and demonstrated favourable outcomes in simulations (Zubeldia, Chluba & Battye 2023).Although the release of a CIB-free Compton -map is planned, it is expected to have higher overall noise variance compared to the Planck PR4 NILC -map due to the additional CIB constraints imposed on NILC.With greatly reduced CIB contamination and minimal overall variance, the PR4 NILC -map may already provide a reliable thermal SZ template for cross-correlation studies.
The last proposed extension is to account for relativistic corrections to the thermal SZ SED in the NILC pipeline by incorporating the average temperature of the Planck clusters, as suggested by Remazeilles et al. (2019).This extension is warranted for the Planck PR4 -map, because even though the relativistic SZ correction is faint, it is still expected to contribute statistically to the signal at Planck sensitivity.results are based on observations obtained with Planck, 6 an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada.This work made use of Python packages like pymaster (Alonso et al. 2019), astropy (Astropy Collaboration et al. 2022), scipy (Virtanen et al. 2020), and matplotlib (Hunter 2007).Some of the results in this paper have been derived using the healpy (Zonca et al. 2019) and HEALPix (Górski et al. 2005) packages.We acknowledge the use of the PSM, developed by the Component Separation Working Group (WG2) of the Planck Collaboration.matic bandpasses and do not include 1/  noise but homogeneous Gaussian white noise with Planck channel sensitivities.
The simulated Planck frequency maps went through the same NILC pipeline as described in Section 3.2 to derive NILC weights and reconstruct the thermal -map.Each individual foreground and noise component of the simulation is propagated with the same NILC weights to reconstruct their respective residual map.Characterizing the residual contamination of the reconstructed -map is then possible by computing the 1-PDF and the angular power spectra of the residual map of each contaminant.
Fig. A1 shows the 1-PDF of the simulated input thermal SZ signal (dashed red), as well as the 1-PDF of the thermal SZ signal reconstructed from NILC (solid black) and that of each residual contaminant (solid coloured lines).Only pixels outside the Galactic plane mask (GAL-MASK) have been considered.Point-sources were not masked because we want to assess the contribution of radio and IR sources separately.In agreement with what we saw in Fig. 14 for Planck PR4 data (Section 5.3), the radio sources (purple line in Fig. A1) contribute dominantly in the negative tail and sparingly to the positive region.This is because radio sources are significant at frequencies below 143 GHz, where the frequency response of the thermal SZ is negative.Hence, they end up as negative residual sources when propagated through NILC.The opposite is largely true for the IR source residuals (brown line); they contribute primarily to the positive tail of the distribution and sub-dominantly to the negative tail.Nevertheless, IR source residuals are low enough so that the characteristic positive non-Gaussian tail of the thermal SZ (Rubiño-Martín & Sunyaev 2003) is almost perfectly recovered as can be seen by the overlapping positive tails of the input and output thermal SZ PDFs.This supports our findings in Section 5.3 and our strategy to mask only radio sources detected in the frequency range where thermal SZ has a negative frequency response.
Residual Galactic foregrounds (orange line in Fig. A1) contribute equally to the negative and positive non-Gaussian tails of the 1-PDF of the reconstructed -map.Finally, excess variance arises mainly from residual noise (blue), CIB (green) and Galactic foregrounds (orange).
Fig. A2 shows the angular power spectra of the input thermal SZ (dashed red), the output NILC thermal SZ obtained from autospectra (solid black) and from cross-spectra (i.e., free from noise bias; dashed black), the noise (blue), and the residual foreground components.Consistently with the PR4 data analysis, we apply the GAL-MASK and the PS-MASK7 to compute the power spectra.As can be seen, all foregrounds are mitigated below the thermal SZ signal by NILC over a wide range of multipoles.The CIB (green) and the noise (blue) are the main residual contaminants at small angular scales (high multipoles) after point-source masking, while Galactic foregrounds (orange) add residual excess power mostly at large angular scales.Residual contamination from CMB (grey) is sub-dominant at all angular scales in the NILC -map.

APPENDIX B: HARMONIC ILC FOR THERMAL SZ
For the sake of comparison, we also implemented a Harmonic Internal Linear Combination (HILC) on the PR4 data.The HILC operates in a spherical harmonic domain with perfect localization but lacks spatial localization in the pixel domain in contrast to NILC.Even though the HILC may perform well in extracting homogeneous and isotropic signals like the CMB (e.g.Tegmark, de Oliveira-Costa & Hamilton 2003), it is definitely not the best option to reconstruct localized signals like the thermal SZ, as we show here.
Fig. B1 (top panel) shows the thermal SZ Compton -map reconstructed with the HILC method using the Planck PR4 data (HILC -map).Except for needlet decomposition, all parameters used are similar to those of the NILC pipeline as described in Section 3.2.The maximum multipole till which a frequency channel is used in the HILC is determined by a cutoff multipole at which the instrumental beam window function of that channel drops to 10 −3 .As can be seen from the comparison between Fig. B1 and Fig. 4, the HILC -map is noticeably noisier than the NILC -map.This is confirmed by the 1-PDF of the -maps in Fig. 9 and the 1-PDF of the noise in Fig. 11, both of which show significantly larger variance for the HILC than for the NILC.The unnormalized skewness of the HILC -map is also larger than that of the NILC -map even after masking with GAL-MASK and PS-MASK.This points towards greater IR source contamination in the former map due to insufficient spatial localization from the HILC.All of this results in the NILC -map having a larger signal-to-noise than the HILC -map for thermal SZ observations.
The power spectrum of the HILC -map is shown in the bottom panel of Fig. B1, along with the power spectrum of the noise associated with the HILC -map as obtained from PR4 half-ring data splits.The HILC -map power spectrum obtained from the crosscorrelation of half-ring maps is also shown (thick black line).In this plot, the excess due to the noise is clearly visible at high multipoles and low multipoles in comparison with the PR4 NILC -map power spectrum.When considering the case of the cross-spectrum of half-ring maps, which is free from noise bias, the HILC -map power spectrum still shows some excess power at high multipoles compared to NILC, which indicates larger contamination from extragalactic foregrounds.Similarly, at very large scales, an additional excess not accounted for by the noise is also visible, which indicates stronger large-scale residual contamination from Galactic foregrounds in the HILC -map than in the NILC -map.This again is due to the fact that HILC lacks the spatial localization to deal efficiently with inhomogeneous and anisotropic foregrounds.
The average signal-to-noise ratio (SNR) is computed over the multipole range ℓ = 30-1000 for both HILC and NILC -maps as This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Masks used in the analysis, including the small processing mask for component separation (NILC-MASK, white area) leaving  sky = 98% of the sky, the Galactic mask from Planck Collaboration XXII (2016) (GAL-MASK, black area) leaving  sky = 60% of the sky and the point-source mask (PS-MASK, black dots) for statistical analysis.
transfer functions of each channel  and  gauss ℓ

Figure 3 .
Figure 3. Gaussian-shaped needlet windows (black) for localization in angular scales.The output 10 ′ beam window of the -map is also shown (blue).

Figure 4 .
Figure 4. NILC thermal SZ Compton -map from Planck PR4 data (top; this work) and Planck PR2 NILC -map (bottom; Planck Collaboration XXII 2016).The left-hand and right-hand sides correspond to the northern and southern hemispheres, respectively, in Galactic coordinates.The black outline shows the boundary of the Galactic mask used for power spectrum analysis.

Figure 6 .
Figure 6.A region at intermediate Galactic latitude in the thermal SZ -maps from Planck PR4 (top) and PR2 (middle) data and their difference PR4 − PR2 (bottom).The improved destriping in the PR4 NILC -map compared to the PR2 NILC -map is visible.

Figure 10 .
Figure10.Thermal SZ angular power spectra from the PR2 NILC -map (red) and the PR4 NILC -map (black), before correction for the noise bias through the auto-power spectrum of the -maps (dashed lines) and after correction for the noise bias through the cross-power spectrum between the HR1 and HR2 -maps (solid lines).The PSM SZ model (thin blue line) is shown as a reference.
Figure 11.1-PDF of residual noise estimate from PR4 (black) and PR2 (red) NILC -maps.Dashed lines show respective Gaussian fits.The noise 1-PDF from the PR4 HILC -map (Appendix B) is also shown for comparison.

Figure 12 .
Figure 12.Angular power spectrum of residual noise in the PR4 and PR2 NILC -maps (top), and relative decrease of noise power in the PR4 -map with respect to the PR2 -map (bottom).
Figure14.1-PDF of PR4 NILC -map without point-source mask (purple), with LFI source masks (yellow), with 30 to 143 GHz source masks (black), and with all LFI and HFI source masks (red).The major source contamination comes from radio sources in LFI channels.The PR4 NILC -map (purple) has a lower level of radio-source contamination than the PR2 NILC -map (dashed black).
Figure A1.1-PDF of NILC -map (solid black) versus input -map (dashed red), and individual contributions from various residual foregrounds (coloured solid lines) for Planck PSM simulations.

Figure A2 .
Figure A2.NILC thermal SZ and residual foreground power spectra for Planck PSM simulations.

Figure B1 .
Figure B1.Results from ILC in harmonic domain (HILC): reconstructed thermal SZ -map (top) and power spectra (bottom).The HILC -map is much noisier compared to the NILC -map and shows more significant residual excess power at both low and high multipoles.

ℓ
is the cross-spectrum of HR1 and HR2 -maps and the uncertainty is derived from the auto-spectrum  ℓ as ℓ = √︄ 2 (2ℓ + 1)  sky Δℓ  ℓ (B2)to account for the sample variance of the signal, residual foregrounds and noise.The SNR is 61.2 for HILC and 178.2 for NILC.This indicates a factor of 3 improvement in SNR for the NILC -map compared to the HILC -map.

Table 1 .
Thermal SZ SED coefficients in both thermodynamic temperature units and intensity units across Planck LFI and HFI channels after Planck PR4 bandpass integration.

Table 2 .
Frequency channels used for each needlet band.
-maps by taking cross spectra with Planck CIB maps at 857 GHz.Top: results with Planck GNILC CIB map at 857 GHz (  sky = 50%).Bottom: results with Lenz et al. (2019) CIB map at 857 GHz (  sky = 15%).The cross-spectrum between CIB noise and -map noise estimates is also shown (dashed-dotted lines).