© The Authors 2025

1. Introduction

Galaxy clusters, the largest gravitationally bound structures in the Universe, harbour a rich variety of physical phenomena. Radio observations have revealed that clusters often show diffuse synchrotron emission that can span megaparsec-sized regions, such as ‘radio halos’ (e.g. Bonafede et al. 2022) or ‘mega-halos’ (Cuciti et al. 2022), implying that clusters are filled with ultra-relativistic electrons and magnetic fields. The influence of the magnetic fields extends to particle acceleration models, radio synchrotron age estimates, the dynamics of the intracluster medium (ICM), and the transport of cosmic rays. Understanding the properties and origins of magnetic fields in clusters thus has broad importance (see Carilli & Taylor 2002; Govoni & Feretti 2004; Donnert et al. 2018, for reviews on magnetic fields in galaxy clusters).

The most promising tool to study magnetic fields is radio polarisation observations. A magnetised plasma such as the ICM causes a wavelength-dependent rotation of the polarisation angle (Faraday rotation) and depolarisation. In general, the Faraday depth of a source is defined as (Burn 1966; Brentjens & de Bruyn 2005)

$$\begin{array}{c}\hfill \varphi (\mathit{r})=812{\displaystyle {\int }_{\mathrm{LOS}}}{n}_{\mathrm{e}}\mathit{B}·\mathrm{d}\mathit{r}{\mathrm{rad}\mathrm{m}}^{-2},\end{array}$$(1)

where n_e is the electron density in parts per cm⁻³, B is the magnetic field in μGauss, and dr is the infinitesimal path length increment along the line of sight (LOS) in kiloparsecs, and we define ϕ(r) > 0 for the magnetic field pointing towards the observer. In the simple case of just one radio-emitting source along the LOS, the Faraday depth is equal to the rotation measure (RM). With a combination of RMs from radio observations and electron densities from X-ray observations, it is thus possible to study the magnetic field properties of galaxy clusters.

Such studies are best done at low redshifts to maximise the cluster angular size because polarised radio sources are relatively rare (e.g. Rudnick & Owen 2014). The most detailed analyses have been of the Coma cluster (Bonafede et al. 2010) and Abell 2345 (Stuardi et al. 2021), where seven polarised radio sources were detected per cluster. The Coma cluster magnetic field was found to agree with a Kolmogorov power spectrum with a central strength of 5 μG and a scaling of magnetic field energy density linearly proportional to the thermal gas density (B² ∝ n_e). The central magnetic field strength in Abell 2345 was found to be similar to the Coma cluster, but with a magnetic field energy density that scales more steeply instead (i.e. B² ∝ n_e²). Several other low-redshift clusters have been analysed in polarisation (Murgia et al. 2004; Govoni et al. 2006, 2010, 2017; Guidetti et al. 2008; Vacca et al. 2012), with typically less than five polarised radio galaxies per study, resulting in large uncertainties on the magnetic field estimates (see e.g. Johnson et al. 2020, for a detailed discussion).

Because cluster magnetic fields are thought to be generally turbulent and disordered, the observed Faraday rotation is the result of a random walk process and thus a random variable. Accurate magnetic field estimates, therefore, require a statistical analysis probing many independent sight lines. Another potential problem is that polarised radio galaxies are often embedded in the cluster, and the degree to which the observed RM variations are caused by local interaction of the lobes with the ICM is debated (Rudnick & Blundell 2003; Laing et al. 2008; Guidetti et al. 2012; Osinga et al. 2022). Such problems can be overcome by stacking clusters to increase the number of polarised radio sources located behind clusters and thus the number of independent sight lines through a cluster (Clarke et al. 2001; Bonafede et al. 2011; Böhringer et al. 2016; Stasyszyn & de los Rios 2019; Osinga et al. 2022). Although stacking experiments have limited ability to probe differences between clusters, they are useful for obtaining average cluster magnetic field properties and are currently the only way to study cluster magnetic fields beyond the few nearest clusters.

We recently published the largest homogeneous stacking experiment using Karl G. Jansky Very Large Array (VLA) observations of 124 galaxy clusters (Osinga et al. 2022) selected from the Planck 2nd Sunyaev-Zeldovich source catalog (PSZ2) (Planck Collaboration XXVII 2016). This study presented the first clear depolarisation trend tracing the radial profile of cluster magnetic fields using over 600 polarised radio sources. While depolarisation traces the smaller-scale structure of the magnetic fields on the scale of the resolving beam (usually kiloparsec), the larger-scale structure (i.e. along the LOS, usually megaparsec) can be inferred from the Faraday rotation of the radio sources. In this paper, we added the information from the Faraday rotation of the same sample of sources to study the large-scale properties of the magnetic fields in galaxy clusters. By jointly fitting both depolarisation and Faraday rotation, we aimed to constrain the average magnetic field strength, scaling with density, and power spectrum. Cosmological calculations were performed assuming a flat ΛCDM model with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7 and the R₅₀₀ radius is the radius enclosing an overdensity of 500 at the cluster redshift.

2. Chandra-Planck ESZ sample

The sample of galaxy clusters is a subset of 124 out of 165 clusters from the Chandra-Planck Legacy Program for Massive Clusters of Galaxies¹ (Andrade-Santos et al. 2021) that have VLA observations presented by Osinga et al. (2022). The full details on the data reduction, polarised source identification and association, and determination of the polarisation properties are presented in the aforementioned paper, but we briefly summarise the important points here and highlight some improvements to the catalogue. All clusters are at low redshift (z = 0 − 0.35) to maximise their apparent size. The analysis of the Chandra data and resulting electron density profiles is presented in Andrade-Santos et al. (2017).

Each of the 124 clusters was observed for ∼40 min in the VLA L-band (1−2 GHz), resulting in typically 20 − 30 μJy beam⁻¹ noise levels at a resolution of 6 − 7″ after data reduction. Polarised sources were identified using RM-synthesis (Brentjens & de Bruyn 2005) and matched to total intensity components and optical counterparts. Optical counterparts were used to determine cluster membership, with spectroscopic redshifts where available and photometric redshifts otherwise. The full details on redshift determination are given in Osinga et al. (2022), where Eqs. (25) and (26) were used to divide sources into background and cluster members. In total 6807 and 819 source components were detected in total and polarised intensity respectively. Throughout this paper we simply refer to these components as ‘sources’. We have fit the following model to the polarised intensity as a function of wavelength λ, which accounts for rotation and depolarisation

$$\begin{array}{c}\hfill P({\lambda }^2)=p_{0}Iexp(-2{\mathrm{\Sigma }}_{\mathrm{RM}}^2{\lambda }^4)exp[2i({\chi }_0+\varphi {\lambda }^2)],\end{array}$$(2)

where p₀ denotes the intrinsic polarisation and Σ_RM² the inferred variance of the RM distribution within the observing beam which models the depolarisation as a function of wavelength². χ₀ is the intrinsic polarisation angle, and I denotes the total intensity model, which was assumed to be a simple power-law of the form I(ν) = I₀ν^α where ν is the frequency and α the spectral index. Equation (2) assumes that (see Sokoloff et al. 1998, for details):

1. The dominant mechanism causing Faraday rotation is an external Faraday screen (i.e. a non-emitting foreground magnetised plasma, such as the ICM).

2. The foreground screen magnetic field changes on scales (at least) smaller than the observing beam, which causes Faraday dispersion inside the beam, resulting in depolarisation.

3. The magnetic field fluctuates as a Gaussian random field, implying its scatter can be completely captured in the parameter Σ_RM.

It has been shown (Osinga et al. 2022) that this model accurately fits 75% of the polarised radio sources detected in this sample, and we discard sources that do not show a proper fit (see below).

We have improved the Monte Carlo Markov chain (MCMC) fitter used by Osinga et al. (2022) to now properly take into account the circular nature of χ₀ during the fitting. Osinga et al. (2022) used a uniform prior on χ₀, ∼𝒰(0, π), which would cause the sampler in some cases to become trapped around the boundary values. We removed this prior on χ₀ and fold the chain back into the range [0, π) after the sampling is completed. We also calculate the mean and spread using circular statistics where we take into account the fact that the angles are distributed on the half-circle ∈[0, π). In this way, the mean denotes the angle of the average vector on the unit circle, and the standard deviation is the spread in angles around the average vector. This agrees with the definition of the simple arithmetic mean and standard deviation when the angles are distributed away from the edges of the domain (see e.g. Mardia & Jupp 2009, for more details on circular statistics).

This paper is accompanied by an updated table of polarised components, shown in Appendix C. We note that this update mainly corrects the quoted mean and uncertainty of the intrinsic polarisation angle χ₀ and most sources have similar best-fit RM (ϕ) and depolarisation (Σ_RM) parameters to within a few per cent. This thus does not significantly impact the results. Finally, we used the same criteria for identifying bad fits as Osinga et al. (2022). All sources with a best-fit χ² value that is > 5σ away from the theoretical distribution, and sources with low signal-to-noise polarised emission resulting in artificially large values of Σ_RM² were flagged. This led to 195 bad fits among the 819 polarised sources identified in 124 galaxy clusters. After excluding three clusters (G121.11+57.01, G115.71+17.52, G139.59+24.18) due to a lack of accurate redshift estimates, a total of 610 sources across 117 clusters were available for subsequent analysis. Table 1 provides an overview of the number counts for various sub-samples and the unique clusters in which these sources were found, aiding in understanding the selection function of the defined subsets in the following analysis.

3. Methods

3.1. Determining RM scatter

We first provide a brief note on nomenclature to ensure clarity regarding the terminology associated with various types of scatter. We adopt the following definitions throughout this paper:

• δRM for the measurement uncertainty on RM.

• σ_RM for the measurable scatter (i.e. standard deviation) in RM between sources.

• Σ_RM for the inferred RM fluctuation across sources associated with beam depolarisation.

• δσ_RM for uncertainty on the scatter in RM (σ_RM).

• δΣ_RM for uncertainty on Σ_RM.

The observed (noisy) RM of a source, RM${}_{\mathrm{obs}}\prime$, not only encompasses RM_cluster induced by the ICM but also contains the intrinsic RM of the source (RM_src), the RM picked up as the light travels through the intergalactic medium (IGM), the RM introduced by the Galactic foreground (RM_gal) and finally the Earth’s ionosphere (RM_ion). Taking into account cosmological redshift, we can define

$$\begin{array}{cc}\hfill {\mathrm{RM}}_{{\mathrm{obs}}^{\prime }}=& \frac{{\mathrm{RM}}_{\mathrm{src}}}{{(1+z_{\mathrm{intr}})}^2}+{\mathrm{RM}}_{\mathrm{IGM}}+\frac{{\mathrm{RM}}_{\mathrm{cluster}}}{{(1+z_{\mathrm{cluster}})}^2}\hfill \\ \hfill & +{\mathrm{RM}}_{\mathrm{gal}}+{\mathrm{RM}}_{\mathrm{ion}},\hfill \end{array}$$(3)

where RM_IGM is the integrated contribution of the IGM, which itself has a redshift dependence (cf. Amaral et al. 2021, Eq. (2)) that is however not relevant for this work as we regard it as a single nuisance parameter.

The scatter in RM between sources, or RM variance, analogously is the sum of the variance of these individual terms, as they are independent contributions. The measurement noise is known from the estimated uncertainties in RM from the MCMC spectral fitting (we use the 16th and 84th percentile of the posterior, see Osinga et al. 2022, for details). The Galactic foreground can be largely removed by correcting the RMs for the Galactic contribution, which has been mapped most recently by Hutschenreuter et al. (2022). The most difficult components to isolate are the variance caused by the IGM and the intrinsic RM variance of radio sources. Studies using low RM variance sources (i.e. away from clusters and groups) that attempt to isolate the extragalactic from the Galactic and intrinsic scatter find that the extragalactic RM (rms) scatter is between ∼1 and 10 rad m⁻² (Oppermann et al. 2015; Vernstrom et al. 2019; O’Sullivan et al. 2020; Pomakov et al. 2022).

In this study, we aim to separate the cluster RM scatter from the other sources of scatter. The cluster RM is often the largest term, as near the centre of clusters, RMs can be of the order of 10³ rad m⁻². Meanwhile, the Galactic contribution to the RM is expected to be small since the cluster sample is selected from the Planck 2nd Sunyaev-Zeldovich source catalog (PSZ2) survey and thus avoids the Galactic plane by design (Galactic latitude |b|> 14°). We subtracted the Galactic contribution (of the order of 10¹ rad m⁻²) using the recent map from Hutschenreuter et al. (2022), propagating the uncertainties in quadrature.

After subtraction of the Galactic contribution, we can define the residual RM as $\mathrm{RRM}={\mathrm{RM}}_{\mathrm{obs}\prime }-{\mathrm{RM}}_{\mathrm{gal}}$, and are thus left with

$$\begin{array}{c}\hfill \mathrm{RRM}=\frac{{\mathrm{RM}}_{\mathrm{cluster}}}{{(1+z_{\mathrm{cluster}})}^2}+{\mathrm{RM}}_{\mathrm{extr}}\end{array}$$(4)

where we defined the last term as the extrinsic RM, which captures the RM contribution that is not from the cluster or the measurement uncertainties (δRRM). The extrinsic RM consists of components from the IGM, the radio source itself, and an ionospheric component, which are difficult to separate. However, we can estimate the scatter caused by RM_extr through sources located far away from clusters. At large radii (r > 3 R₅₀₀), the scatter is expected to be unrelated to the ICM, and caused mainly by the extrinsic scatter between sources and noise scatter due to measurement uncertainties. We find that the standard deviation of the observed RMs corrected for the Galactic contribution is equal to 13 ± 2 rad m⁻² (see Fig. B.1), with the measurement uncertainties accounting for 9 rad m⁻². This gives an estimate of σ_RM, extr = 9 rad m⁻², which is of the same order as the observed dispersion of extragalactic RMs (e.g. Schnitzeler 2010; Vernstrom et al. 2019).

Using this estimate of the extrinsic scatter, we can define our estimator for the RM scatter induced by the ICM, in the observer frame. We define the RM scatter corrected for measurement errors and extrinsic scatter analogously to Dolag et al. (2001), as

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{RRM},{\mathrm{cor}}^{\prime }}=\sqrt{{\sigma }_{\mathrm{RRM}}^2-\frac{{\displaystyle \sum _i^N}{(\delta \mathrm{RRM})}_i^2}{N-1}-{\sigma }_{\mathrm{RM},\mathrm{extr}}^2},\end{array}$$(5)

where σ_RRM is any estimator of the standard deviation, and the sum is taken over all sources i = 1 to i = N to account for varying measurement errors, δRRM, per source, and σ_RM, extr = 9 rad m⁻².

Finally, the cluster rest-frame RM scatter is higher than the observed RM scatter because we artificially measure smaller RMs (by a factor (1 + z)⁻²) due to cosmological redshift. This effect is difficult to account for exactly in a stacking experiment but can be approximated by using the median redshift of the cluster sample, weighted by the number of polarised sources per cluster. Finally, we thus define the corrected RRM scatter (standard deviation) as

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{RRM},\mathrm{cor}}={\sigma }_{\mathrm{RM},{\mathrm{cor}}^{\prime }}{(1+\mathrm{med}[z_{\mathrm{cluster}}])}^2.\end{array}$$(6)

This estimator was calculated in radial bins to measure the RM scatter profile.

3.2. Comparison to models

In the simple scenario of random magnetic field directions in cells of size Λ_c kpc, which have some uniform magnetic field strength and uniform electron density, the observed RM is the result of a random walk process. Because of the central limit theorem, the distribution of RMs is then expected to be a Gaussian distribution with zero mean, and variance given by (e.g. Murgia et al. 2004)

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{RM}}^2(a/R_{500})={812}^2{\mathrm{\Lambda }}_c{\displaystyle {\int }_{\mathrm{LOS}}}{(n_{\mathrm{e}}{B}_{||})}^2\mathrm{d}l,\end{array}$$(7)

where dl is the infinitesimal path length increment along the line of sight (LOS) in kpc, n_e is measured in cm⁻³ and B_|| is the magnetic field strength parallel to the line of sight in μG.

In reality, the magnetic field structure will more closely resemble a random field with fluctuations on many spatial scales, and both the magnetic field strength and electron density will scale with radius. Thus comparing observations to more realistic scenarios requires simulated magnetic fields. We followed the approach explained in Section 4 of Osinga et al. (2022) to generate mock rotation measure and depolarisation images for all clusters in our sample which have X-ray observations available (99/124). We assumed that the clusters have spherical and smooth electron density profiles, with best-fits to the de-projected electron density as a function of radius taken from Andrade-Santos et al. (2017), which are shown in Fig. B.4. In reality, there are likely density fluctuations in the ICM. However, these are difficult to measure and of the order of only a few per cent (e.g. Schuecker et al. 2004; Churazov et al. 2012; Sanders & Fabian 2012). Magnetic field fluctuations are expected to be stronger. Therefore, we have modelled the magnetic field as a Gaussian random field.

An example of a mock RM and depolarisation image is shown in Figure 1 for a random cluster in our sample. To compare the mock images to observations, we sample a mock source for every observed source at the same projected radius. If the source was classified as a background source, the integration is done over the full grid, while for a cluster member, the integration is done assuming the source lies at the midplane of the cluster.

Fig. 1. Simulated Faraday rotation (left) and depolarisation (right) screen for a random cluster from our sample as it would appear at 1.5 GHz with perfect sampling. The parameters are B0 = 5.0 μG, n = 2, η = 0.5 and the dashed and solid circles show 0.5 R500 and 1.0 R500, respectively.

In these models, the magnetic field is assumed to be a three-dimensional Gaussian random field with a single-power law spectrum characterised by the following parameters: B₀, η, n, Λ_min and Λ_max. The first two denote the variables that parameterise the magnetic field strength, assumed to follow

$$\begin{array}{c}\hfill B(r)=B_0{\left(\frac{n_{\mathrm{e}}(r)}{⟨n_{\mathrm{e}}(0)⟩}\right)}^{\eta },\end{array}$$(8)

where ⟨n_th(0)⟩ is the average central electron density of the clusters in our sample, which is equal to 0.056 cm⁻³. A relation between B and n_e is predicted by MHD simulations of cluster magnetic fields (e.g. Dolag et al. 2005; Vazza et al. 2018; Marinacci et al. 2018), with the physical motivation for this relation being related to the fact that magnetic fields and thermal electrons are subject to the same external forces (e.g. turbulence and mergers) or micro-scale plasma processes (e.g. Kunz et al. 2011). Depending on the exact physical conditions in the galaxy cluster, this can give rise to various relations between B and n_e. For example, the compression of a tangled magnetic field under the magnetic flux-freezing condition, yields $B\propto n_{\mathrm{e}}^{2/3}$, while assuming that the magnetic energy remains at a constant ratio with respect to the turbulent energy results in B ∝ ne^1/2 (e.g. Johnson et al. 2020, and references therein). Observationally, this relation has also been found, with values of η roughly between 0 and 1 (e.g. Bonafede et al. 2010; Govoni et al. 2017; Stuardi et al. 2021). We normalise B to follow Eq. (8) after generating the Gaussian random field (see Osinga et al. 2022 for details), which technically breaks the null-divergence property of the cluster magnetic field. However, this minimally impacts the final results as found by Murgia et al. (2004) and has been a standard method of generating cluster magnetic fields in the literature (e.g. Murgia et al. 2004; Bonafede et al. 2010; Vacca et al. 2010; Govoni et al. 2017; De Rubeis et al. 2024).

The last three parameters encode the power spectrum of the magnetic field:

$$\begin{array}{c}\hfill |B_k{|}^2\propto k^{-n},\end{array}$$(9)

between minimum and maximum fluctuations scales that are denoted by Λ_min and Λ_max in image space, respectively. In the simplified picture of scale-by-scale equipartition between the energy density of magnetic fields and turbulent motions (kinetic energy), we expect Kolmogorov turbulence to result in a magnetic field energy spectrum with n = 11/3³ (e.g. Schekochihin et al. 2004). The assumption of scale-by-scale equipartition is not necessarily valid in all galaxy clusters, and observationally, the spectrum was found to take values between n = 1 and n = 4 (Murgia et al. 2004; Govoni et al. 2006, 2017; Guidetti et al. 2008; Bonafede et al. 2010; Vacca et al. 2010, 2012; Stuardi et al. 2021; Osinga et al. 2022). We have computed all models on 1024³ pixel grids, to simulate all clusters homogenously. Here, one pixel represents 3 kpc, and clusters are thus simulated out to about 1.5 R₅₀₀, with a minimum fluctuation scale of Λ_min = 6 kpc. These models will be compared to observations in various ways, as detailed in the next section.

4. Results

4.1. Average magnetic field strength

In Figure 2, we plot the RRM as a function of the projected distance to the nearest cluster centre (denoted as a throughout this work). A clear trend is visible, with the scatter in RRM decreasing with a larger distance to the cluster centre. We determined the corrected RRM scatter for all sources within 1 R₅₀₀, using Eq. (6), and find that the RRM scatter of cluster members and background sources are similar⁴, being 200 ± 33 and 219 ± 66 rad m⁻², respectively, as shown in Table 2. Assuming that clusters are spherically symmetrical on average, or at least have no bias towards our direction, and that cluster members are on average located at the mid-plane of the cluster, we expect the RRM scatter of background sources to be at least $\sqrt{2}$ times that of cluster members at the same projected radius. However, cluster members are found preferentially at smaller radii (e.g. Osinga et al. 2022, Fig. B.1), where the scatter in RM is larger due to generally larger magnetic field strength and electron densities (e.g. Bonafede et al. 2010).

Fig. 2. RRM as a function of the normalised distance to the nearest cluster centre. The median uncertainty in RM is 11 rad m−2 for both cluster members and background sources. The plot is shown on a partially logarithmic scale in Appendix B to show low RM sources more clearly.

The mean value of the RMs agrees with zero as a function of radius, as shown in Figure 3, consistent with random magnetic field orientations along the line of sight. Assuming the simple random walk scenario denoted by Equation (7), we find that the most rudimentary estimate of the line-of-sight magnetic field strength is given by

$$\begin{array}{cc}\hfill \left(\frac{B_{||}}{\mathrm{\mu }G}\right)=& 2.46\left(\frac{{\sigma }_{\mathrm{RM}}/200}{{\mathrm{rad}\mathrm{m}}^{-2}}\right){\left(\frac{n_{\mathrm{e}}/{10}^{-3}}{{\mathrm{cm}}^{-3}}\right)}^{-1}\hfill \\ \hfill & \times {\left(\frac{{\mathrm{\Lambda }}_c/10}{\mathrm{kpc}}\right)}^{-1/2}{\left(\frac{L/1000}{\mathrm{kpc}}\right)}^{-1/2}\hfill \end{array}$$(10)

Fig. 3. Mean of the RRM in bins of the projected normalised distance to the nearest cluster centre. The horizontal error bars indicate the bin edges. Cluster members are shown in blue, background sources in orange, and the combined bins in black.

where L indicates the line-of-sight column length, which will be on average twice as large for background sources as cluster members. If we assume that cells are ordered on scales of 10 kpc with an electron density of 10⁻³ cm⁻³ (e.g. Böhringer et al. 2016), this reduces to

$$\begin{array}{c}\hfill B_{||}=\frac{{\sigma }_{\mathrm{RM}}}{2.57\sqrt{L}}·\end{array}$$(11)

If we assume approximately L = 1000 kpc for cluster members and twice as large for background sources, we find magnetic field strengths averaged within R₅₀₀ equal to 2−3 μG. Although the simple scenario suffices to give an order of magnitude estimate for the average magnetic field in galaxy clusters, it is clear from Figure 2 that the product of the magnetic field strength and electron density is not constant as a function of radius.

As we have information on the radial shape of the electron densities, we plot in Figure 4 the observed rotation measures as a function of ICM electron column density. This plot is less populated, as we now only show sources that are detected within a projected radius a < 2 R₅₀₀, where the Chandra-derived column density values are reliable. The scatter in rotation measure significantly increases with increasing column density, and the preferential sampling of cluster members at high column densities (i.e. low radii) is clearly pronounced. Following Böhringer et al. (2016), we calculated the scatter in rotation measure in bins of column density, from which the average magnetic field strength along the line of sight can be calculated as

$$\begin{array}{c}\hfill \left(\frac{B_{||}}{\mathrm{\mu }\mathrm{G}}\right)=3.801\times {10}^{18}\left(\frac{{\sigma }_{\mathrm{RM}}}{{\mathrm{rad}\mathrm{m}}^{-2}}\right)\left(\frac{N_{\mathrm{e}}}{{\mathrm{cm}}^{-2}}\right){\left(\frac{L}{{\mathrm{\Lambda }}_c}\right)}^{1/2}.\end{array}$$(12)

Fig. 4. RRM as a function of the electron column density for sources detected at a < 2 R500. The plot is shown on a partially logarithmic scale in Appendix B to show low RM sources more clearly.

The RM scatter as a function of column density is shown in Figure 5. The bottom panel of Figure 5 shows the resulting magnetic field estimate from Equation (12). We find that the average magnetic field strength is 0.2 μG $\sqrt{L/{\mathrm{\Lambda }}_c}$, resulting in around 2 μG for typical values of L = 1000 kpc and Λ_c = 10 kpc (Böhringer et al. 2016).

Fig. 5. Corrected standard deviation of RMs in bins of electron column density, for sources detected at a < 2 R500. The bins for the full sample have equal frequency with 71 sources per bin, while the bins for the sub-samples are defined in logarithmically spaced bins to sample similar densities. The horizontal error bars indicate the bins, and points are plotted at the median Ne. The bottom panel shows the magnetic field estimate assuming Equation (12), with a correlation length that is a factor 100 smaller than the line-of-sight distance, resulting in a mean strength of ∼2 μG.

4.2. Radial profile

Though the average magnetic field strength is a useful property to constrain, cluster magnetic fields (and electron densities) show a radial decline. We calculated the scatter in rotation measure in a running window of 50 points as a function of radius. This was computed following Eq. (6), using the interquartile range (IQR) as a robust measure of σ_RRM (i.e. σ_RRM = IQR/1.349). A clear trend of the RM scatter being significantly higher closer to the cluster centre is present in Fig. 2, showing the effect of the turbulent magnetised ICM increasing in density, magnetic field strength or both⁵. Background sources and cluster members show similar scatter profiles, indicating that they are largely tracing the same effect. However, background sources do show an increasing scatter above a ∼ 0.7 R₅₀₀. We checked whether this can be attributed to a change in possible confounding variables in Appendix A (i.e. the uncertainty in RRM, galactic latitude of clusters, redshift of clusters, or polarisation fraction of radio sources), and found that this is the region where we start preferentially sampling high-mass clusters, in contrast to the inner region where low-mass clusters are preferentially sampled (see Fig. A.1). As high-mass clusters provide larger column densities, and generally occupy denser parts of the cosmic large-scale structure with higher chances of intervening groups or dense filaments along the line-of-sight, it is likely that this is the cause of the rising RM scatter as a function of the radius after a ∼ 0.7 R₅₀₀.

Theoretically, we would expect background sources to display a scatter in RM that is on average greater by a factor $\sqrt{2}$, since cluster members are on average located at the midplane of the cluster and polarised light thus travels through half the column that background sources probe. However, from Figure 6 it is clear that such an effect is not observed. To quantify this, we performed a z-test to test two null hypotheses in radial bins: (i) that the scatter of background sources and cluster members are the same, and (ii) that the scatter of background sources is $\sqrt{2}$ times the scatter of cluster members. We found that we could reject neither of the null hypotheses with 95% confidence. This implies that the uncertainties are too large to identify a statistically significant difference between the RM of background sources and cluster members. To increase the number statistics, we also binned all sources into a radial bin bounded by [0, R₅₀₀], and found marginal evidence (p = 0.02) to reject the null hypothesis that IQR(RM${}_{\mathrm{behind}})=\sqrt{2}$ ⋅ IQR(RM_inside). However, this is likely caused by the fact that the cluster members are preferentially detected at smaller radii (median radius 0.34 R₅₀₀) than background sources (median radius 0.54 R₅₀₀), where the scatter is expected to be larger.

Fig. 6. Corrected scatter of observed RRMs as a function of the projected distance to the nearest cluster centre. The running scatter is calculated using a fixed sliding window of 50 points, the horizontal error bar indicates the median width of the sliding bin.

In the simplified model of the turbulent ICM as laid out in Section 3, the variance in rotation measure as a function of the projected distance should be proportional to the line-of-sight integral given in Equation (7). We can thus determine the magnitude of the 3D magnetic field fluctuations and the scaling between magnetic field and thermal electron density from the σ_RRM, cor profile. We sampled sources from our mock RM images created with n = 11/3 (i.e. a Kolmogorov power spectrum), with the same radial sampling as the observed sources. We calculated the simulated scatter profile using the same sliding window, which is shown in Figure 7 for a magnetic field to electron density scaling of η = 0.5 and η = 0.0 (i.e. constant B). The canonical profile with η = 0.5 is not a good fit to the data, as it underestimates the scatter near the outskirts and overestimates the scatter near the centre. The observed profile agrees more with a constant magnetic field, although given the smoothing scale of the sliding window, it is difficult to constrain the value of η precisely. The best-fit central magnetic field strength likely lies between 1 and 10 μG, agreeing with the value from the depolarisation of radio sources (Osinga et al. 2022). We note that only using sources behind clusters gives similar results, as the scatter is similar as a function of radius. However, the relative flatness of the observed profile is caused by the cool-core cluster sample, as investigated in the next section.

Fig. 7. Corrected scatter of RRMs as a function of the projected distance to the nearest cluster centre for modelled clusters. The running scatter is calculated using a fixed sliding window of 50 points, with simulated profiles shown for η = 0.0 and η = 0.5.

4.3. Merging versus relaxed clusters

Differences between the magnetic field properties of relaxed and merging clusters have been tentatively observed in various studies (Bonafede et al. 2010; Stasyszyn & de los Rios 2019; Osinga et al. 2022), but not yet clearly quantified. Following Osinga et al. (2022), we split our sample into merging and relaxed clusters based on the presence of a cool-core (CC) or absence of one (NCC). The corrected scatter in RM as a function of the distance for the NCC and CC clusters is shown in Figure 8. There is a clear difference between the two samples, with the CC sample showing a relatively flat scatter profile, while the NCC sample shows significantly increasing scatter with lower projected radius. Similar behaviour was observed in the depolarisation of radio sources in the same sample, with CC clusters showing a flatter depolarisation profile (Osinga et al. 2022). Additionally, splitting the sample in background sources and cluster members (left panel and right panel) shows that both populations trace similar effects in NCC clusters whereas cluster sources show a significantly flatter profile in CC clusters. It is unlikely that this is an ICM effect, given that it would imply an increasing magnetic field as a function of radius. Given the anomalies in the CC sample, In Fig. 9, we compare the NCC clusters scatter profile to simulated scatter profiles. NCC clusters do show general agreement with a model with η = 0.5 or η = 0.0, although the drop in the central region a < 0.2 R₅₀₀ cannot be explained by our model. This drop is observed in both the background population and sources inside clusters. It is possible that high RRM sources near the projected centre of clusters are depolarised, causing an observational bias against detecting high RRM sources near the centre. We exclude this a < 0.2 R₅₀₀ region in the fit in Fig. B.5, where we find that η = 0.0 and η = 0.5 have comparable χ² values.

Fig. 8. Corrected scatter of RRMs as a function of the projected distance to the nearest cluster centre. Cool-core clusters are shown in blue, while non-cool-core clusters are shown in orange. The left panel shows background sources, while the right panel shows cluster members. The running scatter is calculated using a fixed sliding window of size 0.50, indicated by the horizontal errorbar.

Fig. 9. Corrected scatter of RRMs observed in NCC clusters as a function of the projected distance to the nearest cluster centre. The running scatter is calculated using a fixed sliding window of 50 points, the median width of the sliding window is indicated by the horizontal error bar. Best-fit simulated profiles are shown for η = 0.0 and η = 0.5.

4.4. Constraining the magnetic field power spectrum

Both the RM and depolarisation of sources scale similarly (although not identically) with the magnetic field strength, but differently with n and Λ_max, and therefore the ratio can be used to constrain n and Λ_max, relatively independent of the other parameters (e.g. Bonafede et al. 2010). We defined the ratio of RM to depolarisation as

$$\begin{array}{c}\hfill \frac{|\mathrm{RM}|}{1-\mathrm{DP}}·\end{array}$$(13)

This ratio is also relatively independent of radial distance, as the depolarisation will increase (i.e. take lower values of DP) as the RM increases towards the cluster centre. We verified this with both Spearman and Pearson tests, which showed no significant correlation.

The median observed ratio of |RM|/(1 − DP) across all sources was found to be 113 ± 12 rad m⁻², with the uncertainty calculated by 1000 bootstraps, with background sources and cluster members having similar values of 106 ± 17 rad m⁻² and 119 ± 14 rad m⁻², respectively. Figure 10 shows the median ratio compared with simulated values for sources sampled at similar positions in the simulated RM and depolarisation images. It is clear from the figure that models with most of the magnetic field energy on small scales (i.e. n < 3) cannot reproduce the observed RM ratio, mainly because those models result in too low values of |RM| due to the rapidly fluctuating magnetic field along the line of sight. Instead, models with n ≥ 3 provide a good fit for various values of the maximum correlation scale Λ_max. Lowering Λ_max has an analogous effect to lowering n, namely decreasing the coherence length of the magnetic field along the line of sight. This thus results in a smaller average |RM| while the effect on depolarisation is less significant, as this is measured on scales below the observing beam (less than typically 15 kpc, although dependent on cluster redshift). Thus, the observed data best matches n ≥ 3.

Fig. 10. RM ratio, as defined in Equation (13), as a function of the magnetic field power spectrum index n. The median observed value is indicated by the black line, while the simulated values are indicated by the coloured points, with different markers for different magnetic field strengths. The uncertainty on the simulated values reflect the 16th and 84th percentile of 10 random initialisations.

For a Kolmogorov spectrum (n = 3.67) with typical values of the central magnetic field strength between 1 − 10 μG, the data are consistent with Λ_max > 100 kpc. Lower values of Λ_max would require significantly higher central magnetic field strengths.

4.5. RM and depolarisation forward model

Finally, we find the best model that reproduces the data by following the approach used in previous works such as Murgia et al. (2004), Bonafede et al. (2010), Govoni et al. (2017), Stuardi et al. (2021), directly comparing the simulated RM and depolarisation images to the observed data. For every source, we sample an equivalent source from the simulated clusters, and compute the expected RM and depolarisation. We minimise the difference between the simulated and observed radial scatter in RMs, and the depolarisation as a function of radius. We define the function to be minimised, q, as follows

$$\begin{array}{cc}\hfill q& =q_{\mathrm{depol}}+q_{\mathrm{RM}}\hfill \\ \hfill & ={\displaystyle \sum _r}{\left(\frac{{\mathrm{DP}}_{\mathrm{obs}}(r)-⟨{\mathrm{DP}}_{\mathrm{sim}}(r)⟩}{\mathrm{err}({\mathrm{DP}}_{\mathrm{obs}})}\right)}^2+{\left(\frac{{\sigma }_{\mathrm{RRM},\mathrm{obs}}(r)-⟨{\sigma }_{\mathrm{RRM},\mathrm{sim}}(r)⟩}{\mathrm{err}({\sigma }_{\mathrm{RRM},\mathrm{obs}})}\right)}^2,\hfill \end{array}$$(14)

where the observables are calculated in running bins equivalent to Fig. 6. Simulated observables are averaged over 10 different random initialisations as denoted by the angle brackets. Because sources in this sample were found to be intrinsically depolarised with DP = 0.92 at large radii (Osinga et al. 2022), we incorporated this into the simulated depolarisation. We fix the power spectrum to the Kolmogorov value of n = 3.67 to reduce the computational burden, as Section 4.4 showed that the RM ratio is consistent only with models that have n > 3.

Given the atypical behaviour of the CC sample, we forward model only the NCC sample, taking both behind and cluster members since the profiles are statistically similar (see Sect. 4.3). This smaller sample size consists of 182 sources detected at a < 1500 kpc in 42 unique clusters. First, we investigate the best-fit models when fixing the magnetic field to electron density scaling to the typical value of η = 0.5. Figures 11a and 11b show the values of q_depol and q_RM respectively, as a function of B₀ and Λ_max. We find the best agreement for a model with B₀ = 10 μG and a maximum correlation scale equal to Λ_max = 100 kpc for the depolarisation profile and B₀ = 5 μG with Λ_max = 1536 kpc for the RRM scatter profile. The former is the best-fit model when combining the q values in a normalised form. Figure 12 shows the measured and simulated radial profiles for both best-fit models. Both models are unable to explain the data fully, with both the observed σ_RRM and the depolarisation profile being flatter than the simulated profiles.

Fig. 11. Values of qdepol (a) and qRM (b) as defined in Equation (14) for combinations of B0 and ΛRM. Models are simulated with a Kolmogorov power spectrum and η = 0.5. The best-fit model is marked by a cross.

Fig. 12. Comparison between observed and forward-modelled depolarisation and RM scatter for two models that minimise the q statistic as defined in Equation (14). (a) Model that minimises qdepol. (b) Model that minimises qRM.

Section 4.2 showed that the radial profile of σ_RRM preferred low values of η ≈ 0.0. We thus also forward modelled with η = 0.0, and show the q-values and best-fit forward models in Appendix B, Figures B.6 and B.7. Strictly by q-value, models with η = 0.0 provide a better fit to the data than the η = 0.5 models, but again none of the models provide a consistent fit to both the depolarisation and rotation measure scatter that we observed.

5. Discussion

Previous studies of samples of clusters have investigated either the statistical depolarisation (Bonafede et al. 2010; Osinga et al. 2022) or the rotation measure as a function of the projected radius (Clarke et al. 2001; Böhringer et al. 2016; Stasyszyn & de los Rios 2019). In this work, we have for the first time measured both the rotation measure and depolarisation in a consistent way for polarised radio sources along sight-lines of stacked galaxy clusters.

The previously largest statistical study of cluster magnetic fields was performed by Böhringer et al. (2016), who cross-matched 1773 clusters with the RM catalogue compiled by Kronberg & Newton-McGee (2011) to find 92 RM sight-lines within R₅₀₀. We significantly increase the number statistics, detecting 261 radio sources with well-defined RM and depolarisation within R₅₀₀ in only 124 clusters. Comparing Table 2 with Table 1 in Böhringer et al. (2016), we find significantly higher σ_RRM in all radial bins, except outside R₅₀₀. The reason for this is likely twofold. First, there is a difference in the RM sensitivity as the Böhringer et al. (2016) is likely based on radio data with larger channel widths, preventing detection of high RMs (they find a maximum |RM|< 700 rad m⁻²). Second, the cluster samples are different, as Böhringer et al. (2016) quote a mean cluster mass of 3 × 10¹⁴ M_⊙, while only 15 out of our 124 clusters have a mass below this value, and the mean mass of our sample is 5.7 × 10¹⁴ M_⊙. Taking the different sampling into account by calculating the scatter as a function of the electron column density in Section 4.1, we found an average magnetic field estimate of 3 μG in our sample of galaxy clusters, perfectly consistent with the findings of Böhringer et al. (2016). Outside the typical R₅₀₀ radius of clusters, we observed a corrected standard deviation σ_RRM = 28 ± 4 rad m⁻². This is lower than the value of 56 ± 8 rad m⁻² found by Böhringer et al. (2016) which did not correct for measurement uncertainties, but significantly higher than the expected intrinsic variation in the radio sources of 5 − 7 rad m⁻² (Schnitzeler 2010). We thus observe an enhanced RM scatter beyond R₅₀₀ as well.

Inside R₅₀₀, the shape of the scatter profile was used to determine the magnetic field strength and electron density scaling, as shown in Figure 7. The RM scatter inside R₅₀₀ showed a flat radial profile compared to expectations from a simple scaling of the magnetic field with the electron density of clusters. Plasma theories generally predict η ≥ 0.5, with magnetic flux freezing giving $\eta =\frac{2}{3}$, adiabatic compression giving η = 1.0 and dynamo models often predicting a constant magnetic energy density to thermal energy density ratio (i.e. η = 0.5). Observationally, the best determined magnetic field profile is that of the Coma cluster, for which η = 0.4 − 0.7 was found. Other studies using resolved (cluster) radio galaxies also find η ≥ 0.5 (Murgia et al. 2004; Vacca et al. 2012; Govoni et al. 2017; Stuardi et al. 2021).

In contrast, we found that the best-fit σ_RRM(a) model preferred η < 0.5, below the other experimental and theoretical values. However, statistical studies such as this one, that use RMs of unresolved radio sources to determine the RM scatter, suffer from a significant observational bias as sources with high RM values near the centre of clusters are likely to be depolarised and thus missed. The observed electron density profiles, shown in Fig. B.4, illustrate the problem as the electron densities rise strongly towards the core of clusters. Using Eq. (7), we expect σ_RRM to be ∼8 times larger at r = 0.15 R₅₀₀ than at r = 0.5 R₅₀₀. Since the scatter observed at 0.5 R₅₀₀ is of the order of 50 rad m⁻², we expect a scatter of the order of 400 rad m⁻² at r = 0.15 R₅₀₀. Sources with such high RM scatter are quickly depolarised at the observed frequencies (Osinga et al. 2022), presenting a missing data problem. This artificially decreases the steepness of the σ_RRM profile, requiring low values of η to match the data. In contrast, studies relying on resolved radio sources often probe σ_RRM on smaller scales (i.e. the size of the radio source) where fluctuations are expected to be significantly smaller and thus if a polarised radio source is detected, the scatter can be determined more accurately (e.g. Bonafede et al. 2010).

Other statistical studies probing many unresolved sightlines such as Böhringer et al. (2016), Stasyszyn & de los Rios (2019) have also not observed a strong increase of σ_RRM at low radii, implying the same observational bias. However, at higher frequencies, RM values as high as 10 000 rad m⁻² have been observed in the centre of some clusters⁶ (e.g. Taylor & Perley 1993; Baidoo et al. 2023). Higher frequency data that suffers less from depolarisation is thus needed to determine the value of η accurately. This will be provided by, for example, the VLA Sky Survey (VLASS; Lacy et al. 2020) Depolarisation does not suffer as strongly from this observational bias at low radii, as upper limits on depolarisation fractions can still be set on sources that are significantly depolarised, given that they are detected at high signal-to-noise in total intensity. However, the depolarisation curves that we observed also do not agree with any of the Gaussian random field models explored in this paper, being generally flatter than models, as also seen in Osinga et al. (2022).

Thus, we finally fit both the σ_RRM and depolarisation radial profile jointly by means of forward modelling the NCC clusters, since the CC clusters showed anomalous behaviour as explained below. None of the models could explain the RM and depolarisation data consistently, with observations finding a flatter profile than expected from the models. However, our models assume clusters are spherical over-densities without any surrounding large-scale structure, which is unrealistic since clusters will on average lie between one or two filaments (e.g. Pimbblet et al. 2004). The densities and possibly magnetic field strengths in the regions connecting the clusters to the filaments will be higher than our spherical β-like model, plausibly flattening the observed RM and depolarisation profiles. Additionally, our models do not account for any fluctuations in the electron density with a similar or different power spectrum than that of the magnetic field. Although there are indications from simulations that the magnetic field and thermal density power spectra could be similar at small scales (e.g. Domínguez-Fernández et al. 2019), this is difficult to constrain observationally and expected to be of the order of a few per cent (e.g. Churazov et al. 2012). Additional density fluctuations would generally increase the modelled RM scatter, and thus require slightly lower magnetic field strengths to fit observations. Finally, the assumption that the magnetic fields are Gaussian is in contrast with expectations if the fields are generated by the fluctuation dynamo mechanism, which predicts spatially intermittent fields with non-Gaussian statistics (e.g. Schekochihin et al. 2004; Vazza et al. 2018; Seta et al. 2020; Kopyev et al. 2022). Comparisons with cosmological magneto-hydrodynamical (MHD) simulations are an obvious next step, with realistic turbulent gas profiles, potential non-Gaussian magnetic fields and clusters that are embedded in the large-scale structure around them (e.g. Marinacci et al. 2018).

Keeping these caveats in mind, we fixed the magnetic field to electron density scaling to the observationally best-determined value of η = 0.5 (Bonafede et al. 2011) and found that the best-fit model had B₀ = 5 μG and Λ_max ∼ 100 kpc. The central magnetic field strength is consistent with previous single object studies as well as statistical studies (Murgia et al. 2004; Govoni et al. 2006, 2017; Guidetti et al. 2008; Bonafede et al. 2010; Vacca et al. 2010, 2012; Stuardi et al. 2021; Böhringer et al. 2016). However, the correlation scale is significantly larger than found in resolved cluster member studies that typically find values around 50 kpc (Guidetti et al. 2008; Bonafede et al. 2010; Vacca et al. 2012; Govoni et al. 2017), although the resolved sources are often limited in size to less than a few hundred kiloparsec. Studies that use the brightness fluctuations or possibly polarised emission of radio halos to constrain the magnetic field power spectrum found results consistent with outer magnetic field fluctuation scales of ∼400 kpc (Govoni et al. 2005; Vacca et al. 2010), agreeing with the 100 kpc estimate. Fluctuations on scales of more than a few hundred kpc are also expected theoretically, as the turbulent dynamo process thought to be responsible for magnetic field amplification in galaxy clusters is expected to occur on various scales, from less than a kiloparsec up to a megaparsec (Donnert et al. 2018).

In the forward modelling, we have decided to combine the information from cluster members and background sources to maximise the number statistics, as the σ_RRM profiles are statistically similar, and Osinga et al. (2022) showed that the depolarisation properties of cluster members and background sources are also similar. We checked for a biased contribution to σ_RRM from cluster members in Section 4.2 but found that, using the full sample, we could reject neither the null hypothesis that cluster members and background sources have similar scatter nor the null hypothesis that background sources show a scatter that is $\sqrt{2}$ times larger. This is in accordance with the interpretation from Osinga et al. (2022) that in large samples, cluster members and background sources trace similar ICM effects. However, it is also possible that we do not detect the $\sqrt{2}$ difference because cluster members have slightly enhanced RM scatter due to interaction with the ICM (e.g. Rudnick & Blundell 2003). Because we could not reject either null hypothesis, a larger number of polarised radio sources is needed to observe whether background sources indeed show an RM scatter that is $\sqrt{2}$ times that of cluster members.

Finally, we checked for qualitative differences between merging and relaxed clusters. During cluster mergers, turbulence is injected on large scales, which can amplify magnetic fields, and drive large fluctuation scales (e.g. Vacca et al. 2018). Relaxed clusters, on the other hand, are expected to have smaller fluctuation scales, as the energy injected by a previous merger has dissipated through cascades to smaller and smaller scales. We found a significant difference between the merging and relaxed cluster samples. Merging clusters show higher values of σ_RRM in the region 0.5 < r/R₅₀₀ than relaxed clusters, and sources inside relaxed clusters probe a flat σ_RRM profile as a function of radius. Higher RM scatter in merging clusters was also observed by Stasyszyn & de los Rios (2019), although with larger radial bin sizes. If a scaling between B and n_e is assumed, then CC clusters should show a steeper radial profile near the centre, as CC clusters often have higher central electron densities, as also shown in Fig. B.4. It is unclear what is causing the relative flatness of the σ_RRM profile in the sample of cluster members inside CC clusters, so CC clusters were not used for the full RM and depolarisation forward model. In the future, exploring a two-component magnetic field model for CC clusters might alleviate the inconsistency, analogous to two-component models that are used to fit CC cluster thermal electron density profiles (e.g. Vikhlinin et al. 2006).

6. Conclusion

This work has presented the continuation of the study presented by Osinga et al. (2022), where VLA 1−2 GHz polarisation observations of 124 massive Planck clusters were presented, and the depolarisation properties of polarised sources were investigated as a function of the projected radius to the nearest cluster centre. We have incorporated the additional information from the best-fit RM and constrained cluster magnetic field properties by combining depolarisation and RM in a sample of clusters for the first time. We summarise the results of this work as follows:

1. We have clearly detected the increase of the scatter in rotation measure as a function of decreasing projected radius or increasing electron column density. Averaging all 124 clusters, we find a corrected scatter within R₅₀₀ of σ_RM = 209 ± 37 rad m⁻². The scatter outside of R₅₀₀ was found to be 28 ± 5 rad m⁻², still significantly larger than our measurement uncertainties.

2. Assuming that magnetic fields fluctuate on a single characteristic length scale Λ_c with a constant strength, the observed RM scatter agrees with an average magnetic field strength within R₅₀₀ equal to 2 (Λ_c/10 kpc)^−0.5 μG.

3. Differences between the RM scatter of cluster members and background sources are expected from the differences in path lengths alone. However, using the full sample of sources, we could not reject the null hypothesis that cluster members sources show similar σ_RRM profiles as background sources, consistent with the result that background sources and cluster members also show similar depolarisation in the same sample of clusters (Osinga et al. 2022).

4. The profile of σ_RRM shows a significant decrease as a function of the projected radius. However, the profile is flatter than expected from models that have a scaling of $B\propto n_{\mathrm{e}}^{\eta }$ with η = 0.5. We found that CC clusters show a flatter profile than NCC clusters, particularly when traced by cluster members.

5. Based on the ratio of RRM and depolarisation, models with magnetic field power spectral indices n < 3 are strongly rejected by the data. This is consistent with expectations from Kolmogorov turbulence (n = 3.67).

6. Jointly modelling both the depolarisation and rotation measure of sources in a forward modelling approach, we find that the observations cannot be explained by the simple model of a Gaussian random magnetic field that follows B ∝ n_e^η. Although η = 0 is slightly preferred over η = 0.5, more advanced modelling of cluster magnetic fields is required.

In this work, we implicitly assumed that all clusters have the same magnetic field parameters, while in reality, this might be a function of dynamical state, mass, or redshift. The universality of cluster magnetic fields has not been thoroughly tested (e.g. Govoni et al. 2017) and might significantly affect our results. Future observations with more sensitive telescopes such as MeerKAT and the Square Kilometre Array could test this assumption by detecting enough polarised sightlines through single clusters such that stacking is not required. Such observations should take into account that Equation (2) might not hold for all sources, as diffuse radio emission from clusters (e.g. van Weeren et al. 2019) might induce complexity into the Faraday spectra.

Data availability

Table C.1 is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/694/A44

Acknowledgments

We thank the referee, Amit Seta, for a thorough report that has improved the quality of this manuscript. EO and RJvW acknowledge support from the VIDI research programme with project number 639.042.729, which is financed by the Netherlands Organisation for Scientific Research (NWO). RJvW acknowledges support from the ERC Starting Grant ClusterWeb 804208. EO thanks both Wout Goesaert and Joppe Swart for the careful looks at his code and Bryan Gaensler for providing comments on the manuscript. The Dunlap Institute is funded through an endowment established by the David Dunlap family and the University of Toronto. Basic research in radio astronomy at the U.S. Naval Research Laboratory is supported by 6.1 Base Funding. This paper has made use of observational material taken with an NRAO instrument. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation. The Legacy Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID 2014B-0404; PIs: David Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID 2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; Prop. ID 2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. The Legacy Surveys project is honoured to be permitted to conduct astronomical research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University This research has made use of the NASA/IPAC extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This research has made use of NASA’s Astrophysics Data System (ADS). This research made use of Astropy, a community-developed core Python package for Astronomy. This research has made use of ChatGPT to restructure text and code.

References

Appendix A: Confounding variables

The radial distribution of polarised sources has an important selection effect. The fact that clusters are observed with pointed observations of a finite primary beam size (15′) means that there is a strong correlation between the projected radius of polarised sources and host cluster redshift, as shown in Fig. A.1. This, combined with the fact that PSZ clusters, show a strong correlation between mass and redshift, with lower mass clusters only being found at low redshifts (e.g. Fig. 26 in Planck Collaboration XXVII 2016), means that at small radii, the stacked sample is dominated by low-mass clusters, while at large radii the sample is dominated by high-mass clusters. Thus, even if clusters are self-similar in their thermal properties, at larger radii we will be probing preferentially larger column densities and thus expect larger RM scatter and depolarisation. Therefore, a simple analytical evaluation of Eq. 7 with a fixed path length does not follow the observed trend. However, if the radial profiles are compared to simulated profiles by forward modelling every cluster, integrating along appropriate path lengths, and sampling a simulated radio source at the same radial position, this is taken into account.

Fig. A.1. Distribution of all polarised radio sources in the (foreground) cluster versus projected radius plane. The black line shows the FWHM of the primary beam of the VLA L-band observations for a cluster of R500 = 900 kpc. The grey line shows the location where background sources start preferentially sampling higher redshift (higher mass) clusters, coinciding with the rising profile in Fig. 6.

We also investigate whether other variables than the projected radius correlate with increased scatter in RM. We found that the smallest background sources show the largest values of |RRM|, as shown in Figure A.2. Particularly, for the 29 sources with $\sqrt{\mathrm{Maj}\times \mathrm{Min}}<5$ arcseconds, the |RRM| rises rapidly. These sources show a median |RRM| = 28 rad m⁻², while larger sources show a median of |RRM| = 14.7 rad m⁻². This is not fully unexpected, as smaller sources probe a smaller (turbulent) foreground screen, and are thus less affected by smoothing. However, the results are not dominated by these 29 sources, as they are only a small subset of the whole sample, and the radial profile of the scatter is similar when excluding the smallest sources.

There is a strong correlation between the frational polarisation of radio sources and |RRM|, as shown in Figure A.3. This effect is particularly strong for cluster sources, which show significantly higher |RRM| at low polarisation fractions. Additionally, there are statistically more low polarisation sources inside clusters than behind (2% KS probability of being the same). However, this effect is likely because the strongest depolarisation happens near the centre of the cluster (Osinga et al. 2022), where it is easier to detect cluster members than background sources. It is difficult to separate the effect of an intrinsic relation between low fractional polarisation of sources and higher rotation measures from the assumed effect of the ICM, but we can re-calculate the scatter profile using only sources with a fractional polarisation above 2%, where the relation between |RRM| and p_1.5GHz flattens. This is shown in Figure A.4. The high fractional polarisation sources still show decreasing scatter as a function of radius, with similar profiles for cluster members and background sources. It is only the lowest fractional polarisation sources that show a significant difference between cluster members and background sources. Sources with low fractional polarisation are only found inside clusters when they are located near the centre. However, high fractional polarisation sources dominate the analysis, with 502 sources found with p_1.5GHz > 0.02 and only 107 sources found with p_1.5GHz < 0.02. The profile shown in Fig. 6 is thus dominated by high fractional polarisation sources, which explains that cluster and background sources follow similar profiles.

Finally, there is a strong correlation between |RRM| and Stokes I flux density for cluster members, where brighter sources in total intensity preferentially show higher RRM. Osinga et al. (2022) showed that we preferentially sample brighter sources near the cluster centre (i.e. there is a correlation between total flux density and projected radius from the cluster centre for cluster radio sources), which explains why the most centrally located, and thus brightest sources, correlate positively with |RRM|. In conclusion, the aforementioned correlations are either not dominating the observed trend (i.e. in the case of small sources, or low fractional polarisation) or are likely a by-product of the observed ICM effect (i.e. in the case of the flux density correlation).

Fig. A.2. Running median of the absolute value of RRM as a function of yjr source size. The distribution of sources in this parameter plane is plotted in the background and on the histograms beside the axes.

Fig. A.3. Running median of the absolute value of RRM as a function of the polarisation fraction at 1.5 GHz. The distribution of sources in this parameter plane is plotted in the background and on the histograms beside the axes.

Fig. A.4. Corrected scatter in RRM as a function of the projected radius, separately for sources with low (p1.5GHz < 0.02, 107 sources) and high (p1.5GHz > 0.02, 502 sources) fractional polarisation.

Fig. A.5. Absolute value of RRM as a function of the Stokes I flux density.

Appendix B: Alternative visualisations

In Fig. B.1 we show a zoomed-in version of Fig. 2 beyond a = 3R₅₀₀, used to determine the scatter between sources far outside clusters. We find no significant correlation between RRM and projected radius for sources far outside the central regions, although statistically significant peaks in RRM have been found at large radii in a larger sample of RRMs around cosmic overdensities (Anderson et al. 2024).

To better visualise the low RM sources in our sample, we plot in Figures B.2 and B.3 the cluster-induced rotation measure on a logarithmic scale as a function of the projected radius and column density, respectively. The electron density profiles of the clusters in our sample, as determined in Andrade-Santos et al. (2017), are shown in Figure B.4. Finally, the q values of the full forward model for η = 0.0 are shown in B.6. The best-fit profiles for these fits are shown in Fig B.7.

Fig. B.1. RMobs′ as a function of the normalised distance to the nearest cluster centre for sources located far away from clusters.

Fig. B.2. RRM as a function of the normalised distance to the nearest cluster centre. The y-axis is logarithmic for visualising low RM sources better than in Fig. 6. Uncertainties that are consistent with zero RM are plotted as downward-facing arrows.

Fig. B.3. RRM as a function of the electron column density. The y-axis is logarithmic for visualising low RM sources better than in Fig. 2. Uncertainties that are consistent with zero RM are plotted as downward-facing arrows.

Fig. B.4. Average electron density profile as a function of the normalised 3D radius for the 92 clusters that have Chandra X-ray observations. The mean profile is shown in the black line.

Fig. B.5. Like Fig. 9 but fit using only the data at r > 0.2R500.

Fig. B.6. Values of qdepol (left) and qRM (right) as defined in Equation 14 for combinations of B0 and ΛRM. Models are simulated with a Kolmogorov power spectrum and η = 0.0. The best-fit model is marked by a cross.

Fig. B.7. Comparison between observed and forward-modelled depolarisation and RM scatter for the model with η = 0.0 that minimises the q statistics as defined in Equation 14.

Appendix C: Updated polarised source catalogue

All Tables

All Figures

	Fig. 1. Simulated Faraday rotation (left) and depolarisation (right) screen for a random cluster from our sample as it would appear at 1.5 GHz with perfect sampling. The parameters are B0 = 5.0 μG, n = 2, η = 0.5 and the dashed and solid circles show 0.5 R500 and 1.0 R500, respectively.
In the text

	Fig. 2. RRM as a function of the normalised distance to the nearest cluster centre. The median uncertainty in RM is 11 rad m−2 for both cluster members and background sources. The plot is shown on a partially logarithmic scale in Appendix B to show low RM sources more clearly.
In the text

	Fig. 3. Mean of the RRM in bins of the projected normalised distance to the nearest cluster centre. The horizontal error bars indicate the bin edges. Cluster members are shown in blue, background sources in orange, and the combined bins in black.
In the text

	Fig. 4. RRM as a function of the electron column density for sources detected at a < 2 R500. The plot is shown on a partially logarithmic scale in Appendix B to show low RM sources more clearly.
In the text

Fig. 5. Corrected standard deviation of RMs in bins of electron column density, for sources detected at a < 2 R500. The bins for the full sample have equal frequency with 71 sources per bin, while the bins for the sub-samples are defined in logarithmically spaced bins to sample similar densities. The horizontal error bars indicate the bins, and points are plotted at the median Ne. The bottom panel shows the magnetic field estimate assuming Equation (12), with a correlation length that is a factor 100 smaller than the line-of-sight distance, resulting in a mean strength of ∼2 μG.

In the text

	Fig. 6. Corrected scatter of observed RRMs as a function of the projected distance to the nearest cluster centre. The running scatter is calculated using a fixed sliding window of 50 points, the horizontal error bar indicates the median width of the sliding bin.
In the text

	Fig. 7. Corrected scatter of RRMs as a function of the projected distance to the nearest cluster centre for modelled clusters. The running scatter is calculated using a fixed sliding window of 50 points, with simulated profiles shown for η = 0.0 and η = 0.5.
In the text

	Fig. 8. Corrected scatter of RRMs as a function of the projected distance to the nearest cluster centre. Cool-core clusters are shown in blue, while non-cool-core clusters are shown in orange. The left panel shows background sources, while the right panel shows cluster members. The running scatter is calculated using a fixed sliding window of size 0.50, indicated by the horizontal errorbar.
In the text

	Fig. 9. Corrected scatter of RRMs observed in NCC clusters as a function of the projected distance to the nearest cluster centre. The running scatter is calculated using a fixed sliding window of 50 points, the median width of the sliding window is indicated by the horizontal error bar. Best-fit simulated profiles are shown for η = 0.0 and η = 0.5.
In the text

Fig. 10. RM ratio, as defined in Equation (13), as a function of the magnetic field power spectrum index n. The median observed value is indicated by the black line, while the simulated values are indicated by the coloured points, with different markers for different magnetic field strengths. The uncertainty on the simulated values reflect the 16th and 84th percentile of 10 random initialisations.

In the text

	Fig. 11. Values of qdepol (a) and qRM (b) as defined in Equation (14) for combinations of B0 and ΛRM. Models are simulated with a Kolmogorov power spectrum and η = 0.5. The best-fit model is marked by a cross.
In the text

	Fig. 12. Comparison between observed and forward-modelled depolarisation and RM scatter for two models that minimise the q statistic as defined in Equation (14). (a) Model that minimises qdepol. (b) Model that minimises qRM.
In the text

Fig. A.1. Distribution of all polarised radio sources in the (foreground) cluster versus projected radius plane. The black line shows the FWHM of the primary beam of the VLA L-band observations for a cluster of R500 = 900 kpc. The grey line shows the location where background sources start preferentially sampling higher redshift (higher mass) clusters, coinciding with the rising profile in Fig. 6.

In the text

	Fig. A.2. Running median of the absolute value of RRM as a function of yjr source size. The distribution of sources in this parameter plane is plotted in the background and on the histograms beside the axes.
In the text

	Fig. A.3. Running median of the absolute value of RRM as a function of the polarisation fraction at 1.5 GHz. The distribution of sources in this parameter plane is plotted in the background and on the histograms beside the axes.
In the text

	Fig. A.4. Corrected scatter in RRM as a function of the projected radius, separately for sources with low (p1.5GHz < 0.02, 107 sources) and high (p1.5GHz > 0.02, 502 sources) fractional polarisation.
In the text

	Fig. A.5. Absolute value of RRM as a function of the Stokes I flux density.
In the text

	Fig. B.1. RMobs′ as a function of the normalised distance to the nearest cluster centre for sources located far away from clusters.
In the text

	Fig. B.2. RRM as a function of the normalised distance to the nearest cluster centre. The y-axis is logarithmic for visualising low RM sources better than in Fig. 6. Uncertainties that are consistent with zero RM are plotted as downward-facing arrows.
In the text

	Fig. B.3. RRM as a function of the electron column density. The y-axis is logarithmic for visualising low RM sources better than in Fig. 2. Uncertainties that are consistent with zero RM are plotted as downward-facing arrows.
In the text

	Fig. B.4. Average electron density profile as a function of the normalised 3D radius for the 92 clusters that have Chandra X-ray observations. The mean profile is shown in the black line.
In the text

	Fig. B.5. Like Fig. 9 but fit using only the data at r > 0.2R500.
In the text

	Fig. B.6. Values of qdepol (left) and qRM (right) as defined in Equation 14 for combinations of B0 and ΛRM. Models are simulated with a Kolmogorov power spectrum and η = 0.0. The best-fit model is marked by a cross.
In the text

	Fig. B.7. Comparison between observed and forward-modelled depolarisation and RM scatter for the model with η = 0.0 that minimises the q statistics as defined in Equation 14.
In the text