© The Authors 2023

1. Introduction

Radio halos and radio relics are the two main types of extended megaparsec-scale radio sources in galaxy clusters. They have different features in terms of location, morphology, polarization, and spectral property (see the review of van Weeren et al. 2019). The synchrotron nature of these radio sources indicates that relativistic cosmic rays (CRs) and magnetic fields permeate the intracluster medium (ICM). Of all proposed CR origins, turbulent and shock (re)accelerations are the most plausible in-situ mechanisms for radio halos and relics, respectively (see the review of Brunetti & Jones 2014), where galaxy cluster mergers play an important role in creating the shocks and turbulence.

While the connection between radio relics and ICM shocks has been well established (e.g., Finoguenov et al. 2010; Akamatsu et al. 2015; Urdampilleta et al. 2018), the role that ICM turbulence plays in particle acceleration is less well understood. Theoretical works showed that typical approaches to second-order Fermi acceleration in the ICM include reacceleration of primary and secondary particles by compressive turbulence via transit-time damping (TTD; e.g., Brunetti & Lazarian 2007, 2011; Miniati 2015; Nishiwaki & Asano 2022) or by incompressive turbulence via nonresonant mechanisms (Brunetti & Lazarian 2016; Brunetti & Vazza 2020). In observations, radio halos have been found to be associated with a number of cluster properties. Their radio power P_ν is correlated with cluster X-ray luminosities L_X, temperature, and mass (Giovannini et al. 1999; Kempner & Sarazin 2001; Cassano et al. 2013; Kale et al. 2015; Cuciti et al. 2021). The presence of radio halos is statistically higher in dynamically disturbed clusters (Cassano et al. 2010; Cuciti et al. 2021). Moreover, the dynamic state of clusters can partially explain the scatter in the P_ν–L_X diagram (Yuan et al. 2015; Cuciti et al. 2021). The key to understanding ICM turbulent acceleration is to map turbulent velocity dispersions in the ICM and search for their correlations with radio properties. The direct way of mapping ICM turbulent velocity fields in galaxy clusters uses X-ray emission line broadening (Zhuravleva et al. 2012), which requires high spectral resolution and is beyond the capability of current X-ray imaging spectrometers. The alternative way is using power spectra to measure density fluctuations as a proxy of the turbulent velocity dispersion (e.g., Churazov et al. 2012; Gaspari et al. 2014; Zhuravleva et al. 2014a). The first attempt of connecting turbulent velocity dispersion and radio halo properties was made by Eckert et al. (2017, hereafter E17), who used the power spectral method to measure the velocity dispersion σ_v for 51 galaxy clusters and studied the turbulent Mach number distribution, concluding that P_ν is strongly correlated with σ_v.

The ongoing LOw-Frequency ARray (LOFAR, van Haarlem et al. 2013) Two-metre Sky Survey (LoTSS, Shimwell et al. 2017) is a suitable for a systematic detection of radio halos in the northern hemisphere owing to its unprecedented sensitivity of 0.1 mJy beam⁻¹ at low frequencies (120–168 MHz). In the footprint of the second LoTSS Data Release (LoTSS-DR2 Shimwell et al. 2022), which covers 27% of the northern sky, we found 83 Planck-Sunyaev-Zeldovich DR2 (PSZ2) clusters (Planck Collaboration XXVII 2016) hosting radio halos (Botteon et al. 2022, hereafter Paper I). The LoTSS-DR2-PSZ2 sample provides an excellent opportunity to systematically study the properties of radio halos in a large galaxy cluster sample. In this work, we focus on the X-ray properties and their connections to the radio halo properties of the PSZ2 clusters in the LoTSS-DR2 footprint. The data analysis includes two main parts. In the first part, we compute two morphological parameters that indicate cluster dynamic states and discuss the discrepancy of measurements from different X-ray telescopes. The morphological parameters will be used in a statistical analysis of radio halos and radio relics in forthcoming works (Cassano et al. 2023; Cuciti et al., in prep.; Jones et al. 2023). In the second part, we compute large-scale surface brightness (SB) and ICM density fluctuations. Using the density fluctuations, we estimate the turbulent velocity dispersion and explore its connection with the radio halo power.

This work is organized as the follows. In Sect. 2 we introduce the X-ray sample of the PSZ2 clusters in the LoTSS-DR2 footprint. Section 3 describes the data reduction and spectral analysis methods. Section 4 presents the results and systematic study of morphological parameters. In Sects. 5 and 6, we present the power spectral analysis and compare the radio halo power with the turbulent dissipation rate. We discuss and conclude this work in Sects. 7 and 8. We adopt a Λ cold dark matter cosmology model with cosmological parameters Ω_M = 0.3, Ω_Λ = 0.7, and h₀ = 0.7.

2. X-ray sample

The LoTSS-DR2 footprint covers 309 PSZ2 clusters. We used data from the archival XMM-Newton European Photon Imaging Camera (EPIC) and Chandra Advanced CCD Imaging Spectrometer (ACIS) for X-ray analysis. There are 115 and 110 PSZ2 clusters with Chandra and XMM-Newton observations, respectively. The data availability of individual PSZ2 clusters is listed in Table 1 of Paper I, while image products are available on the project website¹. The locations of clusters with available data are plotted in Fig. 1, and the sample sizes for different analysis are summarized in Table 1.

Fig. 1. Footprint of the LoTSS-DR2 overlaid with locations of the PSZ2 clusters with available X-ray data.

2.1. Sample for the morphological analysis

The summary of the X-ray sample we used for the morphology analysis is described in Sect. 3.4 of Paper I. Briefly, we applied several criteria including field-of-view coverage, observation mode, and image data quality to select the subsample for our analysis. We derived morphological parameters for 140 clusters, 105 of which were observed by Chandra and 98 by XMM-Netwon. Some PSZ2 objects are composed of multiple separate subclusters in X-rays. Taking all extended X-ray sources into account, 107 and 109 subclusters have Chandra and XMM-Newton measurements, respectively. The total number of subclusters with morphological parameters that we present is 150.

2.2. Sample for the power spectral analysis

A power spectral analysis for SB fluctuations requires more counts than the calculation of morphological parameters. Therefore, we considered an additional threshold of > 10⁴ net X-ray counts in the annulus between 100 kpc and r₂₅₀₀ to select a subsample for the SB power spectral analysis. Sixty-nine out of the total 109 XMM-Newton (sub)clusters met the criterion. We excluded several objects from these that are in a complex merger state, which prevents us from modeling their SB profile using a typical double β-model (Cavaliere & Fusco-Femiano 1978). These clusters are PSZ2 G093.94−38.82 ES and EN, which are in a late premerger phase; PSZ2 G124.20−36.48 N and S (Abell 115), which is an offset major merger after first core passage. In addition, we excluded PSZ2 G160.83+81.66 from the analysis because of its high redshift of 0.88 and small angular size. We also checked Chandra archival data. Because Chandra has only one-third of the effective area of XMM-Newton, we searched for clusters with total ACIS-I exposure > 80 ks and found that all clusters that meet this criterion have available XMM-Newton observations. Because we only investigate surface brightness fluctuations on large scales, where the XMM-Newton point spread function (PSF) size is not an issue, we did not include the Chandra data for analysis. Therefore, we have a sample size of 64.

The cluster masses listed in Paper I were retrieved from Planck Collaboration XXVII (2016), where they were estimated from the Compton-y parameter of each PSZ2 object and are close to the total mass for systems with multiple subclusters that are not resolved by Planck. For systems with multiple components in the X-ray images, we searched for mass ratios in the literature to accurately obtain r₅₀₀ values for individual subclusters. For PSZ2 G058.29+18.55 (Lyra complex), we adopted the hydrostatic M₅₀₀s reported by Clavico et al. (2019), which are 3.5 × 10¹⁴ M_⊙ and 2.5 × 10¹⁴ M_⊙ for the E and W subclusters, respectively. PSZ2 G107.10+65.32 (Abell 1758) has a weak-lensing mass of M_500, N = 9.6 × 10¹⁴ M_⊙ and M_500, S = 3.7 × 10¹⁴ M_⊙ for the N and S subclusters, respectively (Monteiro-Oliveira et al. 2017). No mass estimation is available in the literature for PSZ2 G093.94−38.82 (Abell 2572), and only the W subcluster is detected by Planck. We therefore continue using the PSZ2 mass as the mass of the W subcluster.

3. Data reduction and spectral analysis

We used the XMM-Newton Science Analysis Software (SAS) v18.0.0 and Chandra Interactive Analysis of Observations (CIAO) v4.12 (Fruscione et al. 2006) for the data reduction and analysis. The detailed reduction, image processing, and point-source detection methods are described in Sect. 3.4 of Paper I. In this section, we describe our method of spectral analysis and XMM-Newton EPIC-pn non-X-ray background (NXB) scaling. We used the pn filter wheel closed (FWC) version 2019v1.

3.1. XMM-Newton EPIC spectral analysis

We extracted MOS and pn spectra with event selection criteria #XMMEA_EM&&PATTERN< = 12 and FLAG==0&&PATTERN< =4, respectively. Redistribution matrix files and auxiliary response files were generated by the tasks rmfgen and arfgen, respectively.

We used SPEX v3.06 (Kaastra et al. 1996, 2020) for the spectral analysis. Because most of our objects have a temperature k_BT > 2 keV based on the M − k_BT scaling relations (Mantz et al. 2016), we used the atomic database SPEXACT v2.07, which includes fewer lines for a fast calculation. We used the spectral model combination cie1 × red × abs + cie2 × abs + pow, where the two cies are collisional ionization equilibrium models for the ICM and the foreground Galactic halo, red is the redshift of the object, abs is the Galactic absorption, pow is the power law for the cosmic X-ray background (CXB). For cie1, the abundances of metal elements are coupled to Fe, and we set the lower limit to 0.3 proto-solar (Lodders et al. 2009). The temperature for cie2 was fixed to 0.2 keV (Snowden et al. 1998), and the normalization of abs was set to the value from the database nhtot² (Willingale et al. 2013). The photon index of pow was fixed to 1.41 (De Luca & Molendi 2004). We binned the spectra using the optimal binning algorithm (Kaastra & Bleeker 2016) and used the energy range 0.7–7.0 keV for the spectral fitting. The Cash statistics (Cash 1979) was adopted to calculate the likelihood when the parameters were optimized.

3.2. pn background scaling

XMM-Newton observations are strongly affected by soft proton flares. Therefore, we need unexposed regions on the detectors to evaluate the level of the instrumental background. Different from the two EPIC-MOS detectors, there is no clean out of field of view (OoFoV) area in the four corners of the detector (e.g., Zhang et al. 2020; Marelli et al. 2021), that is, the pn NXB level of each observation cannot be estimated using the OoFoV regions.

The particle backgrounds of both XMM-Newton and Chandra show long-term variation that is anticorrelated with solar activity (Gastaldello et al. 2022). We used Chandra ACIS-S3 long-term monitoring data³ as a reference to predict the NXB level of the pn detector for any given epoch. We first fit the ACIS-S3 light curve using a Gaussian process regression method (Ambikasaran et al. 2015) with the George 0.4.0 package⁴. We adopted the product of an exponential squared kernel and a cosine kernel to represent the short-term stochastic and long-term periodic variation. The light curve and the fitted model are plotted in the left panel of Fig. 2.

Fig. 2. EPIC-pn background scaling using the ACIS-S3 light curve. Left: long-term light curve of the Chandra ACIS-S3 particle background. Middle and right: EPIC-pn vs. ACIS-S3 NXB levels for the FF and EFF observation modes, respectively. For both modes, a quadratic model (red) fits the ratio better than a linear model (dashed black).

We compared the pn FWC background 12–14 keV count rate with the predicted ACIS-S3 background count rate at each epoch of the calibration observations. We used a linear model and a quadratic model to fit the diagrams, and χ² was used to evaluate the goodness of fit. We found that for both the full-frame (FF) and extended-full-frame (EFF) observation mode, the diagrams are somewhat better fit by quadratic models (see the middle and right panels of Fig. 2). We therefore applied the two quadratic models to the science observations. For each observation epoch, we first predicted the ACIS-S3 NXB rate using the best-fit Gaussian process regression model, then we calculated the corresponding pn NXB rate in either FF or EFF modes based on the two quadratic models. We list the best-fit parameters for the two quadratic models in Table 2.

We evaluated the uncertainty of this method by calculating the standard deviation of the residuals of the quadratic fitting. The standard deviations are 5.7% and 3.8% for the FF and EFF modes, respectively.

4. Morphological parameters

To investigate the connection between diffuse radio emission and cluster dynamic states in this series of papers, we adopted two X-ray morphological parameters. They are the concentration parameter (Santos et al. 2008),

$$\begin{array}{c}\hfill c=\frac{F(r<r_{\mathrm{core}})}{F(r<r_{\mathrm{ap}})},\end{array}$$(1)

where F is the X-ray photon flux after vignetting correction, r_core is the aperture of the core region, and r_ap is the outer aperture. The second parameter is the centroid shift (Mohr et al. 1993; Poole et al. 2006),

$$\begin{array}{c}\hfill w={\left[\frac{1}{N_{\mathrm{ap}}-1}{\displaystyle \sum _i}{({\mathrm{\Delta }}_i-\overline{\mathrm{\Delta }})}^2\right]}^{1/2}\frac{1}{r_{\mathrm{ap}}},\end{array}$$(2)

where N_ap is the number of apertures, Δ_i is the centroid for the ith aperture, and $\overline{\mathrm{\Delta }}$ is the average centroid.

Following the convention of Cassano et al. (2010), we set r_core = 100 kpc and r_ap = 500 kpc. To determine the centers of the analysis apertures, we smoothed both XMM-Newton and Chandra images and used the maximum-intensity pixel after point-source subtraction as the center of the analysis aperture.

To calculate the parameters, we input σ = 30 kpc Gaussian smoothed Chandra images but unsmoothed XMM-Newton images. The Chandra flux images were generated by subtracting the blank-sky backgrounds that include CXB emission, while the background maps used to generate XMM-Newton flux images were NXB maps. Therefore, we subtracted from the XMM-Newton images a universal constant as the CXB before calculating the morphological parameters. The universal value S_CXB = 2.3 × 10⁻⁶ cts s⁻¹ cm⁻² arcmin⁻² is the mean value of the cluster-free regions beyond r₂₀₀ in the images of z > 0.3 clusters. The scatter of the CXB value in logarithmic space is 0.23 dex. We note that the scatter is not only contributed by the cosmic variance, but also due to imperfect NXB subtraction, and the scatter of the point-source detection limits is due to the different exposure time.

4.1. Individual and combined measurements

The individual XMM-Newton and Chandra measurements of c and w are provided in Table A.1. Observations are available from both the telescopes for 65 objects, for which we adopted a combined value of the two measurements. The combined values are used to investigate the correlations between cluster dynamic states and radio relic and radio halo properties in this series (Jones et al., in prep.; Cassano et al., in prep.; Cuciti et al., in prep.). For both c and w, the combined parameter

$$\begin{array}{c}\hfill {\mathbf{P}}_{\mathrm{comb}}=\frac{1}{2}({\mathbf{P}}_{\mathit{XMM}}+{\mathbf{P}}_{\mathrm{Chandra}}).\end{array}$$(3)

The uncertainty of the combined parameter contains two parts,

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathbf{P}}_{\mathrm{comb}}={(\mathrm{\Delta }{\mathbf{P}}_{\mathrm{comb},\mathrm{stat}}^2+\mathrm{\Delta }{\mathbf{P}}_{\mathrm{comb},\mathrm{sys}}^2)}^{1/2},\end{array}$$(4)

where the statistic uncertainty

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathrm{\Delta }{\mathbf{P}}_{\mathrm{comb},\mathrm{stat}}=\frac{1}{2}\times & (\mathrm{\Delta }{\mathbf{P}}_{\mathit{XMM}}^2+\mathrm{\Delta }{\mathbf{P}}_{\mathrm{Chandra}}^2\hfill \\ \hfill & {+2\mathrm{\Delta }{\mathbf{P}}_{\mathit{XMM}}\mathrm{\Delta }{\mathbf{P}}_{\mathrm{Chandra}})}^{1/2}\hfill \end{array}\end{array}$$(5)

and the systematic uncertainty

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathbf{P}}_{\mathrm{comb},\mathrm{sys}}=\frac{1}{2}\times |{\mathbf{P}}_{\mathit{XMM}}-{\mathbf{P}}_{\mathrm{Chandra}}|\end{array}$$(6)

is half of the measurement discrepancy. The combined measurements as well as the uncertainties have been presented in Table 2 of Paper I.

4.2. Discrepancy in morphological parameters

The uncertainty of the combined measurements in Eq. (4) is dominated by the systematic term, that is, the discrepancy of the measurements by the two telescopes. The Chandra vs. XMM-Newton measurements of c and w of the 65 objects are plotted in Fig. 3, where the insets in each panel illustrate the discrepancy of the measurements.

Fig. 3. Comparison of Chandra and XMM-Newton morphological parameter measurements. Top left and top right: Chandra vs. XMM-Newton measurements of c and w, respectively. The black lines are the diagonal. In each panel, the subplot is the histogram of the ratio of the measurements from the two telescopes. Dashed vertical red and blue lines indicate the mean and 1σ of the distributions, respectively. The large error bars of two Chandrac measurements are due to low count numbers. Bottom left and bottom right: discrepancy of c and w vs. redshift, respectively. In each panel, the two horizontal orange lines as well as the bands indicate the mean and scatter of the discrepancy of the low- and high-redshift populations.

The c measurements from the two telescopes agree well with each other, with a mean deviation of 7% and a scatter of 11%. When we divide the sample into two different redshift ranges, c from Chandra measurements are higher by 15.3% and 5.3% overall than the XMM-Newton measurements, with 1σ scatters of 12.4% and 10.5% for the high- (z > 0.3) and low-redshift (z < 0.3) populations, respectively (see the bottom left panel of Fig. 3). The centroid shift measurements show a larger discrepancy between the two telescopes. The distribution of w_Chandra/w_XMM has a mean of 0.03 dex and a 1σ scatter of 0.34 dex. We did not find a redshift dependence of the ratio. The small discrepancy of c and the large scatter of w agree with a recent study by Yuan et al. (2022), where c and w of clusters in the full archival XMM-Newton and Chandra data were reported. A detailed investigation of the cross-instrument systematic uncertainty is beyond the scope of this work. Nevertheless, we examine a few possible origins of the systematics in Appendix B.

4.3. Relaxation score

Recently, Ghirardini et al. (2022) proposed a novel method for combining the measurements of different morphological parameters into a new parameter, the relaxation score ℛ. The method calculates the joint cumulative probability function in a multidimensional parameter space. In our case, the joint cumulative distribution function in the space of c and w is

$$\begin{array}{c}\hfill \mathcal{R}(c,w)={\displaystyle {\int }_{-\infty }^C{\int }_W^{\infty }{f}_{c,w}(c\le C,w\ge W)\mathrm{d}w\mathrm{d}c,}\end{array}$$(7)

where f_c, w is the joint probability density function. Using this method, we are able to compare the degree of relaxation of clusters within our sample. We use this parameter in the next section to explore the correlation between SB fluctuation and cluster dynamic state.

5. ICM density fluctuations on large scales

5.1. Calculation of 2D surface brightness fluctuations

The evaluation of the SB fluctuations, especially on large scales, is sensitive to the underlying SB model, as illustrated in Zhuravleva et al. (2015) and Bonafede et al. (2018). For many clusters in our sample, the morphologies are clearly eccentric, which means that a spherically symmetric β-model will overestimate the SB fluctuations. Therefore we used an elliptical β-model to fit the SB on large scales. For clusters with bright cool cores, we additionally used a second β-model to fit the core. For all clusters, we also added a constant model for the CXB during the fit. We fit the parameters directly in the 2D plane. The combination of the SB models can be written as

$$\begin{array}{cc}\hfill S_{\mathrm{model}}(x,y)& =\mathrm{Beta}(x,y,x_1,y_1,s_1,r_1,{\beta }_1)\hfill \\ \hfill & +\mathrm{EBeta}(x,y,x_2,y_2,s_2,r_2,{\beta }_2,{\theta }_2,e_2)\hfill \\ \hfill & +C(s_3),\hfill \end{array}$$(8)

where Beta is the 2D β-model describing the cores, EBeta is the 2D elliptical β-model describing the bulk ICM, C is the constant model describing the sky background, x_i, y_i is the center of the ith model in image coordinates, s_i is the SB normalization, r_i is the β-model core radius, and β_i is the β-model slope. The residual map is calculated as

$$\begin{array}{c}\hfill \delta S(x,y)=\frac{[N_{\mathrm{obs}}(x,y)-B(x,y)]/E(x,y)-C}{S_{\mathrm{model}}(x,y)-C},\end{array}$$(9)

where N_obs is the observed count image, E is the vignetting corrected exposure map, and B is the NXB map. To measure the fluctuations contributed by Poisson noise, we simulated Poisson randomizations of model count images and converted them into flux. The simulated noise residual maps can be expressed as

$$\begin{array}{c}\hfill \delta S_{\mathrm{noise}}(x,y)=\frac{[N_{\mathrm{rand}}(x,y)-B(x,y)]/E(x,y)-C}{S_{\mathrm{model}}(x,y)-C},\end{array}$$(10)

where

$$\begin{array}{c}\hfill N_{\mathrm{rand}}\sim \mathrm{Pois}(\lambda =S_{\mathrm{model}}\times E+B_{\mathrm{smoothed}}),\end{array}$$(11)

is the Poisson randomization of the model count image, where B_smoothed is the smoothed NXB map. To minimize the uncertainty from the background, we chose 0.4r₅₀₀ as the outer boundary for the analysis, at which radius the flux from the ICM is approximately a factor of 2 higher than the sum of CXB and NXB.

We used a modified Δ-variance method (Arévalo et al. 2012) to calculate the 2D power spectra of the residual flux maps. This method cleanly compensates for data gaps and allows us to mask out regions of point sources and substructures of mergers. For each cluster, we obtained the power spectrum of the SB fluctuation component P_2D(k)⁵ by subtracting the noise power spectrum from the power spectrum of the residual map, where we used a Monte Carlo approach to simulate 100 noise maps using Eq. (10). At large wavenumbers, the total power spectrum is dominated by the noise component. Therefore, we set a cutoff at the wavenumber where the power of the fluctuation component is twice that of the noise component. The noise-removed SB fluctuation power spectra were converted into 2D amplitude spectra using the equation

$$\begin{array}{c}\hfill A_{2\mathrm{D}}(k)=\sqrt{P_{2\mathrm{D}}(k)2\pi k^2}.\end{array}$$(12)

5.2. A_2D spectra and correlations with other parameters

E17 adopted a fixed scale of 660 kpc to calculate A_2D. To have at least two independent resolved components in the analysis aperture, we adopted scales of 0.4 × r₅₀₀, which covers the range of physical sizes from 450 kpc to 600 kpc and is close to 660 kpc for massive clusters. After applying a wavenumber cut for each cluster, the A_2D spectra of 36 cover cover the wavenumber of (0.4 × r₅₀₀)⁻¹. We provide the results of power spectral analysis in Appendix C, and the results of A_2D at k = (0.4 × r₅₀₀)⁻¹ are listed in the third column of Table C.1. The radio classifications of the 36 clusters are listed in Table 3. Eleven of the 36 objects have extended radio emission that is identified as a radio halo. In 7 of the remaining 25 objects, diffuse emission in the forms of radio relics or sources of uncertain nature is detected. We therefore did not consider them in the following analysis. Even though there are 10 nondetections of diffuse emission clusters, meaningful radio upper limits can be determined for only 3 clusters that are not contaminated by extended radio galaxies or residuals of the subtraction of discrete sources (see the detailed discussion in Bruno et al. 2023).

We compared A_2D at k = (0.4r₅₀₀)⁻¹ with morphological parameters (see Fig. 4). We calculated the Pearson correlation coefficients and corresponding p-values for A_2D–c, A_2D–w, and A_2D–ℛ in logarithmic space. We found that A_2D is marginally anticorrelated with the concentration parameter c with a p-value of 0.021, whereas the p-value of A_2D − w is 0.18, suggesting no correlation. As the combination of c and w, the relaxation score ℛ is also marginally anticorrelated with A_2D, where the p-value 0.057 is mostly driven by the weak anticorrelation between c and A_2D. We conclude that for our sample, the ICM dynamic state is marginally correlated with SB fluctuations at a scale of 0.4 × r₅₀₀, implying that more relaxed clusters tend to have weaker SB fluctuations on large scales.

Fig. 4. Morphological parameters c (left), w (middle), and relaxation score ℛ (right) vs. A2D. The Pearson correlation coefficient and the corresponding p-value are labeled in each panel.

We also explored the correlations between A_2D, radio halo power P_150 MHz, and cluster mass M₅₀₀ (see Fig. 5). The upper limits of the radio power were obtained from Bruno et al. 2023. The p-values of the two pairs are 0.73 and 0.44, respectively, which means that at least in our sample, A_2D is independent of the radio halo power and cluster mass.

Fig. 5. Radio halo power (left) and cluster mass (right) vs. A2D. The marker size in the left panel indicates the cluster mass. The red points in the left panel indicate the upper limits of the radio halo power.

5.3. Turbulent velocity dispersion

Theoretical work illustrated that modest ICM turbulent motions excite isobaric perturbations, when the density fluctuation is proportional to the turbulent Mach number, that is, δρ/ρ₀ ≃ ηℳ_1D (Gaspari et al. 2014). We estimated the turbulent velocity dispersion based on the following assumptions: (1) all surface brightness fluctuations are contributed by turbulent motions, (2) the triggered perturbations are isobaric, (3) the proportionality coefficient η ≃ 1 (Zhuravleva et al. 2014b) holds for both relaxed and merging clusters⁶, and (4) the ICM can be approximated as isothermal in the radius of calculation, that is, we used an average temperature to calculate the sound speed.

We used pyproffit⁷ (Eckert et al. 2020) to recover the 3D density fluctuations from 2D SB fluctuations. The process is the same as described in E17. In short, we constructed an ellipsoid for the 3D density distribution using the elliptical β-model in Eq. (8) and then computed the power spectrum of the normalized emissivity distribution along the line of sight to convert P_2D into P_3D (Churazov et al. 2012). The final A_3D spectrum was converted as

$$\begin{array}{c}\hfill A_{3\mathrm{D}}(k)=\sqrt{P_{3\mathrm{D}}(k)4\pi k^3}.\end{array}$$(13)

The recovered A_3D spectra for the clusters hosting a radio halo are plotted in Fig. 6. Similar to A_2D, we took the value on the scale of k = (0.4 × r₅₀₀)⁻¹. The value of A_3D of each cluster is listed in the fourth column of Table C.1.

Fig. 6. Recovered A3D spectra for clusters with radio halo detection. The dashed line indicates the slope of the Kolmogorov turbulent cascade.

For each cluster, the temperature is measured from a circular region with radius of 0.4 × r₅₀₀ centered at the X-ray centroid and without point sources and the center core-component. We excluded the MOS1 detector from the spectral analysis because it might not cover the full region due to two missing chips. The measured temperatures were obtained following Sect. 3.1 and are listed in the sixth column of Table C.1. We calculated the average sound speed within the region of analysis from the measured k_BT. The average ICM sound speed is $c_{\mathrm{s}}=\sqrt{\gamma k_{\mathrm{B}}T/\mu m_{\mathrm{p}}}\simeq 507.3\times \sqrt{k_{\mathrm{B}}T/\mathrm{keV}}$ km s⁻¹. The 1D Mach number ℳ_1D on the scale 1/k is identical to A_3D(k) assuming η = 1. The 3D velocity dispersion is ${\sigma }_{\mathrm{v},3\mathrm{D}}=\sqrt{3}{\sigma }_{\mathrm{v},1\mathrm{D}}=\sqrt{3}{\mathcal{M}}_{1\mathrm{D}}{c}_{\mathrm{s}}$. The calculated σ_v, 3D values at k = (0.4 × r₅₀₀)⁻¹ are listed in the fifth column of Table C.1. We note that the A_3D values are linearly correlated with the A_2D values, which means that the relations of A_2D we obtained in Sect. 5.2 stand for A_3D and ℳ_1D as well. The scatter in Figs. 4 and 5 is propagated to the relations with A_3D and ℳ_1D. We discuss the scatter of ℳ_1D due to the systematic uncertainties in Sect. 7.1.

By using the estimated 3D turbulent velocity dispersion, similar to E17, we explored its correlation to radio halo power (see Fig. 7). The p-value of the Pearson correlation coefficient is 0.22, suggesting no correlation between radio power and turbulent velocity dispersion for our sample. Moreover, the velocity dispersions of the only three clusters with reliable radio upper limits are not at the lower end of the distribution. In the next section, we further explore the connection between radio halo power and ICM properties in the scenario of turbulent acceleration.

Fig. 7. Radio halo power at 150 MHz vs. turbulent velocity dispersion. The marker size indicates the cluster mass. The upper radio halo limits are presented in red.

6. Connecting radio halo power to turbulent acceleration

By using the radio halo radius r_RH and total mass within the radio halo M_tot(r_RH), Cassano et al. (2007) derived velocity dispersions σ_RH from the gravitational potential and adopted the quantity $M_{\text{tot}}(r_{\text{RH}}){\sigma }_{\text{RH}}^3$ as the approximated turbulent injection rate. In this section, instead of using the approximated turbulent injection power from the cluster total mass, we estimate the turbulent dissipation rate using quantities including the ICM temperature k_BT, ICM mass M_gas, and turbulent velocity dispersion σ_v, as we further test different turbulent reacceleration models.

We assumed a quasi-steady scenario of turbulent acceleration, which means that the total amount of energy loss including synchrotron and inverse Compton scattering of the cosmic microwave background (CMB) is balanced by the energy injection from acceleration. The turbulent dissipation rate per volume is

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{turb}}=C_{ϵ}{\rho }_{\mathrm{gas}}{\sigma }_{\mathrm{v},k}^{3}k,\end{array}$$(14)

where the coefficient C_ϵ is calculated from the Kolmogorov constant, although its value has been found not to be universal (e.g., Sreenivasan 1995). We adopted C_ϵ ≃ 5 (Zhuravleva et al. 2014a). The dissipation rate itself is the total flux of kinetic energy loss, where kinetic energy can be converted into heat, magnetic energy, and relativistic particles. When the turbulence is of a Kolmogorov nature, σ_v, k ∝ k^−1/3 and the term ${\sigma }_{\text{v},k}^{3}k$ is constant when k is in the inertial range of the turbulent cascade. For the A_3D spectra of our radio halo sample (see Fig. 6), the slope is close to -1/3, therefore we assumed the Kolmogorov nature and used the measurements at k = (0.4r₅₀₀)⁻¹ to estimate the dissipation rate. The total turbulent dissipation power in the volume of the radio halo is

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathrm{turb}}={\displaystyle {\int }_{V_{\mathrm{RH}}}{C}_{ϵ}{\rho }_{\mathrm{gas}}{\sigma }_{\mathrm{v},k}^{3}k\mathrm{d}V.}\end{array}$$(15)

Assuming the coefficient C_ϵ and ${\sigma }_{\text{v},k}^{3}k$ are invariant throughout the volume of the radio halo, we can write Eq. (15) as

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathrm{turb}}=C_{ϵ}{\sigma }_{\mathrm{v},k}^{3}k{M}_{\mathrm{gas}}(r_{\mathrm{RH}}),\end{array}$$(16)

where r_RH is the radius of the radio halo. The turbulent dissipation power is then proportional to the gas mass inside the volume of radius r_RH. Because only < 10% of the total turbulent flux goes into particle acceleration, we introduced a coefficient C_acc to denote the proportion of dissipation to particle acceleration, which is also assumed to be invariant throughout the radio halo volume. Therefore, the injected turbulent power for particle acceleration is

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathrm{inj}}=C_{\mathrm{acc}}{\mathcal{P}}_{\mathrm{turb}}.\end{array}$$(17)

Because of the energy loss of CMB in inverse Compton and the redshift dependence of the CMB luminosity, the energy that goes into synchrotron emission is

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathrm{syn}}=\frac{B^2}{B^2+B_{\mathrm{CMB}}^2}\times {\mathcal{P}}_{\mathrm{inj}},\end{array}$$(18)

where B is the ICM magnetic field strength, and B_CMB = 3.2(1 + z)² μG is the equivalent magnetic strength of the CMB inverse Compton.

We calculated the gas mass inside the volume of the radio halos in our sample by modeling the SB profiles. We extracted the SB profile for each cluster and fit it using a two-component projected density β-model (Cavaliere & Fusco-Femiano 1978), where the hydrogen number density n_H and projected surface brightness S_X are expressed as

$$\begin{array}{cc}& n_{\mathrm{H}}(r)={\displaystyle \sum _i^{n=2}{n}_{\mathrm{H},i}{[1+{\left(\frac{r}{r_{\mathrm{c},i}}\right)}^2]}^{-3{\beta }_i/2},}\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill & S_X(r)=2\times {\displaystyle {\int }_0^{\infty }{n}_{\mathrm{H}}{\left(\sqrt{l^2+r^2}\right)}^2\mathrm{\Lambda }\mathrm{d}l,}\hfill \end{array}$$(20)

where Λ is the cooling function and is approximately a constant for k_BT ≳ 2.5 keV gas in the 0.5–2.0 keV band. The k-correction of each cluster was calculated using calc_kcorr in Sherpa. The gas density can be converted from the hydrogen number density as ρ_gas ≃ 2.3n_Hμm_H, where μ ≃ 0.6 is the mean molecular weight. We integrated the gas mass using the best-fit density profile up to the radius of r_RH. Following the convention of this series, we used three e-folding radii as r_RH, where the e-folding radii of all radio halos were presented in Table 3 of Paper I. Of the 64 objects with deep XMM-Newton exposures, 21 have radio halo detections. The estimated M_gas(r_RH) of the 21 clusters is listed in the last column in Table C.1.

We first considered a fiducial scenario in which a constant fraction of the turbulent energy flux goes into particle acceleration. Because we did not compare the absolute values of turbulent power going into synchrotron emission but only compared the best-fit slopes in logarithmic scale, the value of C_acc is not important and was set to 0.05. We adopted a universal magnetic field strength B = 5 μG to calculate 𝒫_syn for our sample. The impact of the magnetic field strength is investigated in Sect. 7.3. The resulting P_150MHz–𝒫_syn diagram is plotted in the left panel of Fig. 8. Different from the result in the P_ν–σ_v diagram, the corresponding p-value of the Pearson coefficient is improved to 0.020 in the P_ν–𝒫_syn plane, which shows that the radio monochromatic power at 150 MHz is marginally correlated with the injected power from turbulent dissipation. We used the code BCES⁸, which uses the method taking bivariate correlated errors and intrinsic scatter (BCES) into account (Akritas & Bershady 1996), to calculate the slope for our sample. The slope from the BCES bisector method is 1.09 ± 0.42, which is close to unity. Alternatively, the BCES orthogonal method returns a slope of 1.31 ± 0.98, whose uncertainty is much larger than the BCES bisector method. This result of the unity slope agrees with the value of 1.24 ± 0.19 from Cassano et al. (2007), although we used a more detailed calculation and radio data at a different frequencies.

Fig. 8. Radio halo power at 150 MHz vs injection power from turbulent acceleration with constant energy input (left), TTD (middle), and ASA (right), respectively. Solid and dashed lines represent the best-fit results of BCES bisector and BCES orthogonal, respectively. The orange band shows the 1σ confidence band of the BCES bisector fit. The marker size indicates the cluster mass.

In addition to the constant acceleration coefficient, we further considered two different turbulent acceleration mechanisms, that is, TTD (Brunetti & Lazarian 2007) and adiabatic stochastic acceleration (ASA; Brunetti & Lazarian 2016; Brunetti & Vazza 2020). Assuming the ratio of the energy densities of the CR and the thermal ICM is constant across the cluster sample, the acceleration coefficients of the two models are dependent on sound speed and turbulent Mach number, which are C_acc, TTD ∝ c_s × ℳ_1D and $C_{\text{acc},\text{ASA}}\propto c_{\text{s}}\times {\mathcal{M}}_{\text{1D}}^{-1}$ for the TTD and ASA scenario, respectively (see Appendix E). For TTD, we therefore modified the acceleration coefficient as

$$\begin{array}{c}\hfill C_{\mathrm{acc},\mathrm{TTD}}=\left(\frac{c_{\mathrm{s}}}{c_{\mathrm{s},\mathrm{norm}}}\right)\left(\frac{{\mathcal{M}}_{1\mathrm{D}}}{{\mathcal{M}}_{1\mathrm{D},\mathrm{norm}}}\right)\times C_{\mathrm{acc}}^{\prime },\end{array}$$(21)

and for adiabatic stochastic compression, it is

$$\begin{array}{c}\hfill C_{\mathrm{acc},\mathrm{ASA}}=\left(\frac{c_{\mathrm{s}}}{c_{\mathrm{s},\mathrm{norm}}}\right){\left(\frac{{\mathcal{M}}_{1\mathrm{D}}}{{\mathcal{M}}_{1\mathrm{D},\mathrm{norm}}}\right)}^{-1}\times C_{\mathrm{acc}}^{″}.\end{array}$$(22)

We note that in Eqs. (21) and (22), the turbulent Mach number is at a fixed scale, which can be scaled from 0.4r₅₀₀ assuming a Kolmogorov slope. Meanwhile, the normalization terms c_s, norm and ℳ_1D, norm and the two constants $C_{\text{acc}}^{\prime }$ and $C_{\text{acc}}^{\prime }$ in the two equations are arbitrary because the exact values are nontrivial to calculate. We adopted the mean values of the 11 clusters for c_s, norm and ℳ_1D, norm and fixed the values of $C_{\text{acc}}^{\prime }$ and $C_{\text{acc}}^{\prime }$ to 0.05. The P_150 MHz − 𝒫_syn diagrams of the two different acceleration models are plotted in the middle and right panels of Fig. 8. The BCES bisector slopes of TTD and ASA are 0.93 ± 0.31 and 1.04 ± 0.29, respectively. Both slopes are close to unity within the uncertainty. In other words, we cannot distinguish the two different acceleration models with our sample. We also note that the smaller scatter of the ASA scenario compared to the TTD scenario arises because it contains less systematic uncertainties from the ℳ_1D measurement. It does not mean that the data favor the ASA model.

By comparing the two best-fit slopes using Eqs. (21) and (22), we find that the slope does not change due to the large scatter of ℳ_1D. When we substitute C_acc in Eq. (17) with either Eqs. (21) or (22), the turbulent flux that is tunneled into particle acceleration can be written as

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathrm{inj}}\propto c_{\mathrm{s}}^4{M}_{\mathrm{gas}}(r_{\mathrm{RH}})f({\mathcal{M}}_{1\mathrm{D}}),\end{array}$$(23)

where f(ℳ_1D) is an ℳ_1D dependent function, which is different in the two scenarios. When we ignore f(ℳ_1D), that is, when we assume it as a constant, Eq. (23) can further be simplified as

$$\begin{array}{cc}\hfill {\mathcal{P}}_{\mathrm{inj}}& \propto c_{\mathrm{s}}^4{M}_{\mathrm{gas}}(r_{\mathrm{RH}})\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill & \propto {\left[k_{\mathrm{B}}T\right]}^2{M}_{\mathrm{gas}}(r_{\mathrm{RH}}).\hfill \end{array}$$(25)

The new quantity [k_BT]²M_gas(r_RH) suggests that the ICM sound speed (temperature) and mass within radio halo volume are the two main factors behind the turbulent power for particle acceleration. This quantity can be also written as [k_BT ⋅ Y_X]_rRH, which is a product of temperature and the well-known mass proxy Y_X (Kravtsov et al. 2006) within the radio halo radius. We used all 21 clusters in Table C.1 with radio halo detections to calculate [k_BT ⋅ Y_X]_rRH. For the temperature measurements, we directly adopted the k_BT_0.4r500 measurements in Table C.1, whose measurement radii are close to r_RHs. The difference of the measured temperatures due to the different radius adoptions is only at the percent level (Mantz et al. 2016); see also Appendix D. The diagram of P_150 MHz vs. [k_BT ⋅ Y_X]_rRH is plotted in Fig. 9, where the scatter is much smaller than in Fig. 8, and the p-value of the Pearson coefficient is 1.05 × 10⁻⁶. The best-fit BCES bisector slope is 1.18 ± 0.18, which agrees with the unity slope.

Fig. 9. Radio halo power vs. quantity [kBT ⋅ YX]rRH. The marker color indicates cluster M500, and the marker sizes denotes the radio halo radius. Solid and dashed lines represent the best-fit results of BCES bisector and BCES orthogonal, respectively. The orange band is the 1σ confidence band of the BCES bisector fit.

7. Discussion

7.1. Systematic uncertainties and scatter of ℳ_1D

The systematic uncertainty of the turbulent Mach number estimation using the method of SB fluctuation power spectrum has two main origins. The first origin lies in the assumption that all SB fluctuations on top of the underlying model are from turbulent motions, where the fluctuations, especially on a large scale, are determined by the choice of the underlying SB model. A simple circular β-model will overestimate the density fluctuation, whereas a “patched” model (Zhuravleva et al. 2015) might underestimate the density fluctuation. The analyses in Zhuravleva et al. (2015) and Bonafede et al. (2018) demonstrated that the systematic uncertainty due to model choice might be even larger than 50%. Although we adopt the elliptical β-model to fit the global underlying SB, which is intermediate compared to the circular β-model and the patchy model, it is still possible for a given cluster that the density fluctuation is either overestimated or underestimated. This might be the reason why we only find marginal correlations between A_2D and morphological parameters, and no correlation between A_2D and P_150 MHz.

The second origin is the assumption that the density fluctuations is proportional to ℳ_1D with a unity slope. The exact slope may be different from scale to scale and from system to system. Zhuravleva et al. (2014b) reported a scatter along the scale k of 30%. The scatter of the slope at a fixed scale of eight simulated clusters in Simonte et al. (2022) is about 16%. This systematic uncertainty additionally increases the scatter in our P_150 MHz–[k_BT ⋅ Y_X]_rRH plots, and artificially results in a better correlation in the ASA scenario than in the TTD scenario. The tight correlation in the P_150 MHz–[k_BT ⋅ Y_X]_rRH plot implies that the scatter of the true values of ℳ_1D may be much smaller.

Because the scatter of ℳ_1D due to the systematics is large for the subsample of clusters that host radio halos, we cannot reject the null hypothesis that the radio halo power P_150 MHz is independent of ℳ_1D. The consistency of the best-fit slopes in the diagrams of P_150 MHz − [k_BT ⋅ Y_X]_rRH, P_150 MHz − 𝒫_syn of the ASA scenario and the TTD scenario also suggests that P_150 MHz is independent of ℳ_1D, or has at most a weak dependence on f(ℳ_1D).

7.2. Comparison with previous studies

E17 first applied the SB power spectral analysis to investigate the connection between turbulent velocity dispersion and radio halo properties. In this section, we compare our work to that of E17 in terms of sample properties and results.

The radio halo sample used in E17 was adopted from Cassano et al. (2013), where the mass range was M₅₀₀ > 6 × 10¹⁴ M_⊙ and the radio observation frequency was 1.4 GHz. This work uses a Planck-SZ selected sample with the mass range extended to ∼3 × 10¹⁴ M_⊙ and radio observations at 150 MHz. In addition to the different mass ranges, the two samples have different ranges of the radio halo power. The LoTSS-DR2 radio halo sample has a median radio halo power of 1.5 × 10²⁵ W Hz⁻¹. When we assume a typical radio halo spectral index −1.3, the expected median value at 1.4 GHz is 8.1 × 10²³ W Hz⁻¹, which is lower by about a factor of three than the median value of the sample used in E17. Similarly, the median radio halo power of the 11 clusters with both A_2D measurement and radio halo detection in this work is 1.3 × 10²⁵ W Hz⁻¹, whose expected median radio halo power at 1.4 GHz is also lower by about a factor of three than that in E17.

The analysis of E17 has two main results. First, the A_2D distribution shows a bimodality, in which clusters with radio halos have higher SB fluctuations than clusters with only upper limits of a radio halo detection. Second, the radio halo power at 1.4 GHz is correlated with the turbulent velocity dispersion with a best-fit slope of 3.3 ± 0.7.

Our analysis cannot reproduce the first result directly because of the 36 objects with A_2D measurements, only 11 have radio halo detections and 3 have a sufficient radio-image quality to estimate upper limits. We cannot place reasonable upper limits on the other 22 objects either because the radio-image quality is poor or because the emission is significantly affected by calibration artifacts. Nevertheless, the anticorrelation we found in the c–A_2D plane (see Sect. 5.2) indirectly proves that clusters hosting radio halos have higher A_2D than clusters without a radio halo. Statistical studies showed that the occurrence of radio halos is significantly higher in clusters with low concentration parameters (e.g., Cassano et al. 2010; Cuciti et al. 2015). Clusters with higher SB fluctuations are more likely to host less dense cores and therefore have a higher probability to host radio halos. However, the bimodality shown in E17 might be due to the nonuniformity of the sample, which includes clusters from the two flux-limited X-ray selected samples REFLEX (Böhringer et al. 2004) and eBCS (Ebeling et al. 2000) and clusters with radio halos reported in the literature (see Cassano et al. 2013 for a sample description).

Different from the second conclusion of E17, our analysis does not find a correlation between P_150 MHz and σ_v. This is simply due to the small sample size of 11 and the large scatter of ℳ_1D (or A_2D). If the radio halo power P_ν is independent of ℳ_1D, the underlying M − k_BT and M − P_ν scaling relations will result in a weak σ_v − P_ν correlation. On the other hand, the strong correlation reported by E17 might be amplified by the two most luminous radio halos, that is, the Bullet cluster and MACSJ0717, whose X-ray morphologies are extremely disturbed, and the corresponding σ_v could be overestimated. When we exclude the two clusters, the p-value of the Pearson correlation coefficient drops from the original 2 × 10⁻⁶ in E17 to 0.05. Moreover, the observation frequency and selection function of the LoTSS-DR2 radio halo sample are different from those of Cassano et al. (2013). The sample of Cassano et al. (2013) contains more luminous radio halos, and our sample is likely to contain more ultra-steep spectrum radio halos (USSRHs), which are hard to detect at higher frequencies. The properties of the most luminous radio halos might be different from USSRHs in terms of the P_ν − σ_v relation.

7.3. Unity slope of the P_150 MHz − 𝒫_syn relation

Both this work and Cassano et al. (2007) investigated the slope in the diagram of radio halo power vs. turbulent flux. The quantity that Cassano et al. (2007) used to denote turbulent flux only takes the gravitational potential and the radio halo size into account. From a macroscopic view of energy conservation, it is clearly expected that the gravitational potential energy is eventually converted to heat, magnetic energy, and relativistic particles. The detailed astrophysical processes and channels that convert the gravitational potential energy need to be investigated to interpret the observed phenomena, however. For this reason, our study went one step deeper and focused on the baryonic contents within the ICM. The quantity we used, 𝒫_syn, was calculated using mass, turbulent Mach number, and sound speed of the ICM based on turbulent dissipation and detailed acceleration models. Although we cannot distinguish the TTD and ASA scenarios, the unity slope we obtained between the expected synchrotron emission that is originally from the turbulent acceleration and the observed radio halo power further supports the theory of turbulent (re)acceleration.

When we calculated 𝒫_syn, we only tentatively used a fixed magnetic field strength of 5 μG. We additionally tried B = 1 μG and B = 10 μG for the constant energy input scenario, and the corresponding slopes are 1.10 ± 0.51 and 1.08 ± 0.40, respectively, which means that a different choice of magnetic field strength will not significantly affect the result due to the large scatter of the ℳ_1D. Even when we ignore the dependence on B and on M_1D, the slope in the P_150 MHz − [kT_B ⋅ Y_X]_rRH diagram remains unity, which implies that for our sample, whose redshift median is ∼0.2, the parameter of B does not affect the result. Future studies of high-redshift samples could shed light on the impact of the magnetic field strength.

8. Conclusion

We analyzed archival XMM-Newton and Chandra X-ray data of 140 PSZ2 clusters in the footprint of LoTSS-DR2. We computed two morphological parameters, the concentration parameter and the centroid shift. For 36 clusters that were observed with deep exposures we used the power spectral analysis and measured the amplitudes of surface brightness and density fluctuations at the scale of 0.4r₅₀₀. We also estimated the turbulent velocity dispersion at the same scale. Using the turbulent velocity dispersion, we calculated the turbulent dissipation rate, investigated the relation between turbulent flux and radio halo power, and tested different acceleration models. Our results are summarized below.

1. The measurements of the concentration parameter obtained with the two telescopes agree well with each other with a global discrepancy of 7 ± 11%. In contrast, the discrepancy of the centroid shifts from the two telescopes is large, with an rms of 0.34 dex.

2. We found a marginal correlation between the surface brightness amplitude A_2D and concentration parameter. However, we did not find correlations between A_2D and cluster mass and radio halo power, which further implies that the turbulent Mach number could be independent of the cluster mass and radio halo power.

3. The flux of turbulent acceleration that goes into synchrotron radiation is well correlated with the radio halo power with a unity slope. The two acceleration mechanisms, transit-time damping and nonresonant adiabatic stochastic acceleration, cannot be distinguished within the uncertainties of the slopes.

4. We introduced a new quantity [k_BT ⋅ Y_X]_rRH, which is easy to calculate and denotes the turbulent acceleration flux assuming a constant turbulent Mach number. The quantity [k_BT ⋅ Y_X]_rRH is well correlated with radio halo power, where the slope is also unity. This quantity can be applied to different samples in the future to verify whether the slope and scatter remain consistent with the tight direct proportionality reported here.

The purpose of this project was to explore the connection between the nonthermal properties of radio halos and the gas dynamics as well as the thermal contents of galaxy clusters. Future high spectral resolution X-ray observations that directly measure turbulent velocity dispersions using microcalorimeters and radio observations that cover a wide frequency range will deepen our understanding of the particle acceleration and radio halo formation.

Acknowledgments

The authors thank the anonymous referee providing useful comments that improved the paper. XZ thanks Marco Simonte for sharing the simulation results and acknowledges the support from China Scholarship Council. AS is supported by the Women In Science Excel (WISE) programme of the Netherlands Organisation for Scientific Research (NWO), and acknowledges the Kavli IPMU for the continued hospitality. SRON Netherlands Institute for Space Research is supported financially by NWO. FG and MR acknowledges support from INAF mainstream project ‘Galaxy Clusters Science with LOFAR’ 1.05.01.86.05. AB acknowledges support from the VIDI research programme with project number 639.042.729, which is financed by NWO, and from the ERC-StG DRANOEL n. 714245. RJvW acknowledges support from the ERC Starting Grant ClusterWeb 804208. MB and FdG acknowledge support by the Deutsche Forschungsgemeinschaft under Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306. AD acknowledges support by the BMBF Verbundforschung under the grant 05A20STA. This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. This research has made use of data obtained from the Chandra Data Archive and software provided by the Chandra X-ray Center (CXC) in the application packages CIAO and Sherpa. LOFAR data products were provided by the LOFAR Surveys Key Science project (LSKSP; https://lofar-surveys.org/) and were derived from observations with the International LOFAR Telescope (ILT). LOFAR (van Haarlem et al. 2013) is the Low Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, which are owned by various parties (each with their own funding sources), and which are collectively operated by the ILT foundation under a joint scientific policy. The efforts of the LSKSP have benefited from funding from the European Research Council, NOVA, NWO, CNRS-INSU, the SURF Co-operative, the UK Science and Technology Funding Council and the Jülich Supercomputing Centre. This research made use of the Dutch national e-infrastructure with support of the SURF Cooperative (e-infra 180169) and the LOFAR e-infra group. The Jülich LOFAR Long Term Archive and the German LOFAR network are both coordinated and operated by the Jülich Supercomputing Centre (JSC), and computing resources on the supercomputer JUWELS at JSC were provided by the Gauss Centre for Supercomputing e.V. (grant CHTB00) through the John von Neumann Institute for Computing (NIC). This research made use of the University of Hertfordshire high-performance computing facility and the LOFAR-UK computing facility located at the University of Hertfordshire and supported by STFC [ST/P000096/1], and of the Italian LOFAR IT computing infrastructure supported and operated by INAF, and by the Physics Department of Turin university (under an agreement with Consorzio Interuniversitario per la Fisica Spaziale) at the C3S Supercomputing Centre, Italy. This research made use of Astropy, a community-developed core Python package for Astronomy http://www.astropy.org (Astropy Collaboration 2013, 2018).

References

Appendix A: Morphological parameters of the sample

We list the morphological parameters c and w of the 150 individual clusters in Table A.1.

Appendix B: Systematic uncertainties of the morphological parameters

B.1. Discrepancy in concentration parameter

The PSF of the telescopes is one of the main origins of the discrepancy in c, especially for distant cool core clusters, that is, a large PSF smooths the core and leads to an underestimation of c. The result of the high-redshift population agrees with this explanation. This discrepancy can be corrected when c is recovered from a surface brightness profile that takes the instrumental PSF into account (e.g., Lovisari et al. 2017). However, for low-redshift objects, the effect of the PSF is not expected to be important. In our analysis, we already smoothed the Chandra image with a 30 kpc kernel before the calculation, which was not applied to the XMM-Newton image. This approach makes the smoothness of the Chandra images comparable to the XMM-Newton images at z ∼ 0.3 and even higher for objects at lower redshifts, which means that the PSF is not the only effect that adds to the observed discrepancy. Therefore, we additionally checked the systematic uncertainty due to CXB subtraction for low-zXMM-Newton clusters. We examined the discrepancy when CXB levels of 170% and 60%, respectively, were used, which corresponds to the 0.23 dex scatter of the CXB values of the high-redshift population (see Sect. 4). The corresponding discrepancies are plotted in Fig. B.1. Universal 60% or 170% CXB levels can decrease or increase the measured c with median shifts of 2.0% and 3.5%. This analysis suggests that for our low-zXMM-Newton subsample, the CXB level could be globally higher than the universal value we used, which is obtained from the high-z subsample. This might be due to the large angular sizes of the low-z clusters, where more point sources are hidden behind the ICM emission and are not detected. This effect will be stronger for XMM-Newton observations because its PFS size is one order of magnitude larger than that of Chandra, and therefore its sensitivity to point sources in a cluster field is reduced.

Fig. B.1. Impact of the CXB level adoption on the c measurement. Top: Comparisons between c calculated using 170% (orange) and 60% (blue) CXB levels and the original level for XMM-Newton clusters. The solid line is the diagonal. Bottom: Histograms of the discrepancy of the c values with 170% (orange) and 60% (blue) CXB levels. The dashed lines denote the median values of the two distributions.

B.2. Discrepancy in centroid shift

As shown by Nurgaliev et al. (2013), a low count number (< 2000) may result in systematically overestimating the centroid shift. We investigated this possible bias by plotting the discrepancy in logarithmic space versus the Chandra net count number in the analysis aperture in Fig. B.2. We adopted the Chandra count number because it is lower overall than the XMM-Newton count number. The Pearson correlation coefficient of 0.04 and the corresponding p-value of 0.75 suggest no trend of w bias as a function of count number. We further selected the sources with the largest discrepancy to investigate the origin of the difference. Five sources have discrepancies larger than 2σ_sys. Of these, G187.53+21.92 and G192.18+56.12 have much larger w_XMM, while G172.74+65.30, G092.69+59.92, and G066.41+27.03 have much larger w_Chandra.

Fig. B.2. Discrepancy of w in the logarithmic space vs Chandra net count number in the analysis aperture.

For clusters with w_Chandra ≫ w_XMM, G092.69+59.92 is faint in the shallow Chandra image, which could lead to a large uncertainty. For the remaining two objects, we checked the coordinates of the aperture centers in maps of the two instruments and found large distances between them (see the left and middle panel of Fig. B.3). The two clusters do not host bright cool cores, which means that the uncertainty of the maximum intensity pixel is based on the count number. In addition, the count numbers of the XMM-Newton images are much larger than the Chandra images, suggesting that the X-ray peaks of the Chandra images have a large uncertainty, which leads to overestimates of w. We note that due to the flat morphology of the two clusters, the measurements of c agree with each other within a 10% level even though the X-ray peaks determined by the two telescopes are different.

Fig. B.3. XMM-Newton images of examples with wChandra ≫ wXMM (left and middle) and wChandra ≪ wXMM (right). Triangles and points indicate the maximum intensity coordinates of the Chandra and XMM-Newton images, respectively. The dashed circle masks the residual of a point source near the bright core.

For the two w_Chandra ≪ w_XMM objects, we found that PSZ2 G187.53+21.92 has a peaked morphology and there is a residual of a point source that is not perfectly removed near the core in the XMM-Newton image (see the right panel of Fig. B.3). Because the PSF of XMM-Newton is significant, the traditional point source removing process⁹ cannot work perfectly due to the large gradient of the ICM emission when a bright point source is near the cool core. However, we have no clear explanation for the discrepancy of G192.18+56.12.

Appendix C: Results of the power spectral analysis of the sample

We list the results of the power spectral analysis, temperature, and gas mass of the sample in Table C.1.

Appendix D: Temperature measurements of the sample

We plot the mass versus temperature in Fig. D.1. Although our spectral extraction region is 0.4r₅₀₀, the measurements are close to the M₅₀₀ − k_BT₅₀₀ scaling relation (e.g., Lovisari et al. 2015; Mantz et al. 2016). The typical radio halo radius is in the range of 0.4 to 1.0 r₁₀₀, which means that the k_BT_0.4r500 measurements can be used as the emission-weighted temperatures within r_RHs.

Fig. D.1. Mass vs temperature of our sample. The overplotted lines are the scaling relations of Lovisari et al. (2015) (red) and Mantz et al. (2016) (green).

Appendix E: Fractions of different turbulent acceleration models

In reacceleration models, a fraction C_acc (efficiency) of the turbulent energy flux $F~\rho {\sigma }_v^3{L}^{-1}$ is assumed to be converted into (re)acceleration of relativistic electrons and positrons in the ICM (e.g., Brunetti & Lazarian 2007),

$$\begin{array}{c}\hfill C_{\mathrm{acc},\mathrm{e}}F\sim {\displaystyle \int }{\mathrm{d}}^{3}pE\frac{\partial f_{\mathrm{e}}(p)}{\partial t},\end{array}$$(E.1)

where f_e(p) is the relativistic-electron distribution in momentum space. Under the assumption that the isotropy of the pitch-angle distribution of relativistic particles is preserved during the (re)acceleration, we can use the Fokker-Planck equation to link the right side of E.1 to the coefficient of particle diffusion in momentum space, D_pp (e.g., Brunetti & Lazarian 2007),

$$\begin{array}{c}\hfill C_{\mathrm{acc},\mathrm{e}}\sim F^{-1}{\displaystyle \int }{\mathrm{d}}^{3}p\frac{E}{p^2}\frac{\partial }{\partial p}\left(p^2{D}_{\mathit{pp}}\frac{\partial f_{\mathrm{e}}}{\partial p}\right),\end{array}$$(E.2)

where in the case of TTD (Brunetti & Lazarian 2007; Miniati 2015),

$$\begin{array}{c}\hfill \frac{D_{\mathit{pp}}}{p^2}\propto \frac{c_{\mathrm{s}}^2{\mathcal{M}}_{\mathrm{turb}}^4}{L},\end{array}$$(E.3)

and in the case of nonresonant (re)acceleration with incompressive turbulence assuming a fixed energy flux of magnetic hydrodynamic turbulence is channeled into magnetic field (Brunetti & Lazarian 2016; Brunetti & Vazza 2020),

$$\begin{array}{c}\hfill \frac{D_{\mathit{pp}}}{p^2}\propto \frac{c_{\mathrm{s}}^2{\mathcal{M}}_{\mathrm{turb}}^2}{L}.\end{array}$$(E.4)

Combining Eq. E.2 with Eqs. E.3 and E.4, we can estimate how the efficiency scales with the relevant physical quantities in the two models,

$$\begin{array}{c}\hfill C_{\mathrm{acc},\mathrm{e},\mathrm{TTD}}\propto {\sigma }_{v,k}\left(\frac{U_{\mathrm{e}}}{\rho c_{\mathrm{s}}^2}\right),\end{array}$$(E.5)

$$\begin{array}{c}\hfill C_{\mathrm{acc},\mathrm{e},\mathrm{ASA}}\propto \frac{c_s^2}{{\sigma }_{v,k}}\left(\frac{U_{\mathrm{e}}}{\rho c_{\mathrm{s}}^2}\right),\end{array}$$(E.6)

where $U_{\text{e}}/(\rho c_s^2)$ is essentially the ratio of the energy densities of relativistic electrons and thermal plasma in the ICM. Assuming it is a constant, we can simplify the two equations and obtain

$$\begin{array}{cc}\hfill C_{\mathrm{acc},\mathrm{e},\mathrm{TTD}}& \propto c_{\mathrm{s}}{\mathcal{M}}_{\mathrm{turb}},\hfill \end{array}$$(E.7)

$$\begin{array}{cc}\hfill C_{\mathrm{acc},\mathrm{e},\mathrm{ASA}}& \propto \frac{c_{\mathrm{s}}}{{\mathcal{M}}_{\mathrm{turb}}}.\hfill \end{array}$$(E.8)

All Tables

All Figures

	Fig. 1. Footprint of the LoTSS-DR2 overlaid with locations of the PSZ2 clusters with available X-ray data.
In the text

	Fig. 2. EPIC-pn background scaling using the ACIS-S3 light curve. Left: long-term light curve of the Chandra ACIS-S3 particle background. Middle and right: EPIC-pn vs. ACIS-S3 NXB levels for the FF and EFF observation modes, respectively. For both modes, a quadratic model (red) fits the ratio better than a linear model (dashed black).
In the text

Fig. 3. Comparison of Chandra and XMM-Newton morphological parameter measurements. Top left and top right: Chandra vs. XMM-Newton measurements of c and w, respectively. The black lines are the diagonal. In each panel, the subplot is the histogram of the ratio of the measurements from the two telescopes. Dashed vertical red and blue lines indicate the mean and 1σ of the distributions, respectively. The large error bars of two Chandrac measurements are due to low count numbers. Bottom left and bottom right: discrepancy of c and w vs. redshift, respectively. In each panel, the two horizontal orange lines as well as the bands indicate the mean and scatter of the discrepancy of the low- and high-redshift populations.

In the text

	Fig. 4. Morphological parameters c (left), w (middle), and relaxation score ℛ (right) vs. A2D. The Pearson correlation coefficient and the corresponding p-value are labeled in each panel.
In the text

	Fig. 5. Radio halo power (left) and cluster mass (right) vs. A2D. The marker size in the left panel indicates the cluster mass. The red points in the left panel indicate the upper limits of the radio halo power.
In the text

	Fig. 6. Recovered A3D spectra for clusters with radio halo detection. The dashed line indicates the slope of the Kolmogorov turbulent cascade.
In the text

	Fig. 7. Radio halo power at 150 MHz vs. turbulent velocity dispersion. The marker size indicates the cluster mass. The upper radio halo limits are presented in red.
In the text

	Fig. 8. Radio halo power at 150 MHz vs injection power from turbulent acceleration with constant energy input (left), TTD (middle), and ASA (right), respectively. Solid and dashed lines represent the best-fit results of BCES bisector and BCES orthogonal, respectively. The orange band shows the 1σ confidence band of the BCES bisector fit. The marker size indicates the cluster mass.
In the text

	Fig. 9. Radio halo power vs. quantity [kBT ⋅ YX]rRH. The marker color indicates cluster M500, and the marker sizes denotes the radio halo radius. Solid and dashed lines represent the best-fit results of BCES bisector and BCES orthogonal, respectively. The orange band is the 1σ confidence band of the BCES bisector fit.
In the text

	Fig. B.1. Impact of the CXB level adoption on the c measurement. Top: Comparisons between c calculated using 170% (orange) and 60% (blue) CXB levels and the original level for XMM-Newton clusters. The solid line is the diagonal. Bottom: Histograms of the discrepancy of the c values with 170% (orange) and 60% (blue) CXB levels. The dashed lines denote the median values of the two distributions.
In the text

	Fig. B.2. Discrepancy of w in the logarithmic space vs Chandra net count number in the analysis aperture.
In the text

	Fig. B.3. XMM-Newton images of examples with wChandra ≫ wXMM (left and middle) and wChandra ≪ wXMM (right). Triangles and points indicate the maximum intensity coordinates of the Chandra and XMM-Newton images, respectively. The dashed circle masks the residual of a point source near the bright core.
In the text

	Fig. D.1. Mass vs temperature of our sample. The overplotted lines are the scaling relations of Lovisari et al. (2015) (red) and Mantz et al. (2016) (green).
In the text