© E. Osinga et al. 2022

1. Introduction

Through observations of diffuse synchrotron emission such as radio halos (e.g., van Weeren et al. 2019, for a recent review) and Faraday rotation measures (RMs) of polarised radio sources (e.g., Akahori et al. 2018, for a recent review), it has been proven that galaxy clusters have magnetic fields. These fields play a key role in many astrophysical processes such as heat conduction, gas mixing, and cosmic ray propagation, but the exact properties and origin of these magnetic fields are still uncertain (see Carilli & Taylor 2002; Donnert et al. 2018, for reviews on magnetic fields in galaxy clusters). Estimates of the magnetic field strength from observations of diffuse synchrotron emission (i.e. radio halos) place the magnetic field strengths of galaxy clusters around the μG level (Ferrari et al. 2008). Recently, observations of high-redshift radio halos have revealed that clusters at z > 0.6 might have similar magnetic field strengths to local galaxy clusters (Di Gennaro et al. 2021a), implying that magnetic field amplification should happen fast during cluster formation. However, estimates of the magnetic field strength from diffuse synchrotron emission require various assumptions as to the energy spectrum and distribution of relativistic particles (e.g., equipartition or minimum energy; Beck & Krause 2005).

The most promising method to derive magnetic field properties in clusters is through Faraday rotation of polarised radio emission (see Govoni & Feretti 2004, for a review). Various studies have constrained the magnetic field strength and structure of nearby galaxy clusters using the RM of extended radio sources (e.g., Murgia et al. 2004; Govoni et al. 2006, 2017; Guidetti et al. 2008, 2010; Laing et al. 2008; Bonafede et al. 2010; Vacca et al. 2012). These studies have found central magnetic field strengths of the order of 1–10 μG, and a magnetic field power spectrum index between n = 2 and n = 4.

The depolarising effect of Faraday rotation can also be used to constrain magnetic field properties (e.g., Taylor et al. 2006; Bonafede et al. 2011; O’Sullivan et al. 2019; Stuardi et al. 2020; Sebokolodi et al. 2020; Di Gennaro et al. 2021b; de Gasperin et al. 2022; Rajpurohit et al. 2022). Since we observe radio sources with a finite spatial resolution, the differential Faraday rotation between different lines of sight within a single beam reduces the observed degree of polarisation. This beam depolarisation effect depends on the correlation scales of the magnetic field and the magnetic field strength. In this way, the average properties of magnetic fields in clusters can be investigated, and differences can be studied between various cluster properties, such as the presence or absence of a cool core (Bonafede et al. 2011). The advantage of using fractional polarisation over the RM of radio sources is that unpolarised sources can also be taken into account, as upper limits on the polarisation fraction can be estimated.

A drawback in most studies of Faraday rotation and the resulting depolarisation is that the polarised radio sources are often cluster members. This introduces a small uncertainty because the location of the radio sources inside the cluster cannot be determined accurately, but a larger uncertainty is introduced by the gas in the intracluster medium (ICM), whose properties are usually not known in detail. Often, it is assumed that the interaction between the ICM gas around the radio source and the radio plasma is negligible. However, it is debated to what extent this assumption is true, with some studies showing evidence for local Faraday rotation being induced in radio lobes (e.g., Rudnick & Blundell 2003) and other studies finding no evidence for this (e.g., Ensslin et al. 2003). The ICM could be locally compressed around cluster radio sources, causing higher densities and thus also higher depolarisation, potentially biasing results. Additionally, bent-tailed radio galaxies which are often seen in clusters might not have the same intrinsic polarisation as classic double-lobed radio galaxies (Feretti et al. 1998).

In this paper, we aim to alleviate these problems through a study of the polarisation properties of sources inside and behind clusters. Because the number of polarised radio sources (behind clusters) is typically low (e.g., Rudnick & Owen 2014), such a study will be most often statistical. Although polarisation properties of sources behind clusters have been investigated for some single clusters (e.g., Bonafede et al. 2010), this has not yet been studied thoroughly in a sample of clusters. This paper focuses on the beam depolarisation effect and considers only the implications of the fractional polarisation measurements of the radio sources. In a follow-up paper, we extensively study the Faraday RM of the polarised sources.

Samples of galaxy clusters can be selected relatively unbiased through the thermal Sunyaev-Zel’dovich effect, which imprints a redshift-independent distortion on the spectrum of the cosmic microwave background (Sunyaev & Zeldovich 1970, 1972). The sample we use as a starting point for this work is the Planck Early Sunyaev Zel’dovich sample (Planck Collaboration VIII 2011). This provides mass-selected samples of galaxy clusters up to high redshifts. We obtained observations with the Karl G. Jansky Very Large Array (VLA), which are detailed in Sect. 2, to study the linear polarisation properties of radio sources located inside and behind ESZ clusters. The source finding, determination of the polarisation properties, host galaxy identification and redshift estimation process is explained in Sect. 3. Theoretical depolarisation expectations are derived through modelling of the magnetic fields as Gaussian random fields in Sect. 4 and results are shown in Sects. 5 and 6. Finally, we conclude with a discussion and summary in Sects. 7 and 8. Possible biases are discussed in Appendix A. Throughout this paper, we assume a flat ΛCDM model with H₀ = 70 kms⁻¹ Mpc⁻¹, Ω_m = 0.3 and Ω_Λ = 0.7. We refer to the intensity of linearly polarised light simply as the polarised intensity.

2. Data

2.1. Chandra-Planck ESZ sample

The Planck Early Sunyaev-Zel’dovich (ESZ) results presented 189 cluster candidates all-sky (Planck Collaboration VIII 2011), of which 163 clusters are at a redshift of z < 0.35. The Chandra-Planck Legacy Program for Massive Clusters of Galaxies¹ observed all 163 clusters with sufficient exposure time to collect at least 10 000 source counts per cluster. This makes it (one of) the largest relatively unbiased samples of galaxy clusters with high-quality X-ray data available. The Chandra observations of 147 clusters from the ESZ sample are presented in Andrade-Santos et al. (2017, 2021), where the sample has been reduced by 16 because six clusters are too close to point sources, nine clusters are classed as multiple objects and one system was too large to allow for a reliable background estimate in the Chandra field of view. High-quality X-ray data is particularly important for polarisation studies to be able to break the degeneracy between electron density and magnetic field. In this paper, we used the thermal electron density profiles which were calculated from the fitting procedure detailed in Andrade-Santos et al. (2017).

2.2. Observations and data reduction

We have obtained VLA L-band (1–2 GHz) observations of 126 Planck clusters at z < 0.35 and Dec > −40° (VLA project code 15A-270). The redshift cut is made because the angular size of higher redshift clusters on the sky becomes too small to find a significant number of polarised background sources. Out of these 126 clusters, 102 are from the ESZ catalogue and 24 clusters are new detections in the PSZ1 (Planck Collaboration XXXII 2015) and PSZ2 (Planck Collaboration XXVII 2016) catalogues that have been added to the sample. The observations were taken in the B(nA) array configuration, with the BnA configuration employed for targets located at Dec < −15° or Dec > +75° to match the resolution of targets observed in favourable declination ranges. Targets are observed for ∼40 min each. The full L-band comprises 16 spectral windows, each consisting of 64 channels before frequency averaging.

The calibration of the radio data was done using the Common Astronomy Software Application (CASA; McMullin et al. 2007) and proceeded in the following fashion. For each observing run, the initial data calibration was done per spectral window, such that bad spectral windows could be identified and flagged. We Hanning smoothed the spectral axis to reduce the effect of Gibbs ringing due to strong radio frequency interference (RFI) in the L-band. Shadowed antennas were flagged and the initial flagging of RFI was done with the CASATFCrop algorithm. The effect of the elevation on the antenna gain and efficiency was calculated and antenna position corrections were applied. The flux scale was set to the Perley & Butler (2017) scale. We calculated initial-bandpass calibration solutions using a large solution interval and initially calibrated the complex gains with the central channels of the spectral window. The antenna delay terms were then calculated and applied, after which the final-bandpass solutions could be calculated. A polarised calibrator (either 3C138 or 3C286) was used to solve for a global cross-hand delay and an unpolarised calibrator (3C147) was used to calibrate on-axis polarisation leakage. Subsequently, the polarised calibrator was then used to calibrate the polarisation angle. Off-axis polarisation leakage due to a time, frequency, and polarisation-dependent primary beam becomes important as the distance from the pointing centre increases but is known to be less important in Stokes Q and U than in Stokes V (Uson & Cotton 2008). Typically for VLA L-band observations, the leakage from Stokes I into Q and U is around 1% at the primary beam full-width half maximum (Jagannathan et al. 2017). While this effect can mimic depolarisation due to the frequency dependence of the primary beam, we do not consider it to be a major issue for this study because all clusters are observed near the pointing centre. We discuss off-axis leakage in more detail in Sect. 7.4. Finally, the antenna-based complex gain solutions were calculated using the calibrator sources, and another round of automatic flagging was performed using the CASATFcrop and Rflag algorithms. All spectral windows were then combined and the resulting data were averaged to 8 MHz channels and 6-second timesteps. Leftover RFI was then flagged with the AOflagger (Offringa et al. 2012) and a custom strategy to flag RFI in the cross-hand correlation (rl,lr) plane was employed. Spectral window 8 was fully lost to RFI in every observing run, resulting in a total of 90 frequency channels after initial calibration.

To remove residual amplitude and phase errors in the direction of the target fields and increase the quality of the final images, we further performed six rounds of self-calibration, automatically calculating the solution interval based on the mean flux density in each field. This was done to ensure enough signal-to-noise during the calibration steps, with larger solution intervals used for fields with fainter sources. The imaging and cleaning were done using WSclean version 2.7.3 with the options -join-polarizations and -squared-channel-joining for Stokes Q and U imaging (Offringa et al. 2014). The six rounds of self-calibration involve three phase-only calibration rounds and three amplitude and phase rounds, decreasing the solution interval each round. For the majority of targets, this automatic self-calibration pipeline proved sufficient to obtain high-quality images of the target fields. A small number of target fields needed manual tweaking of parameters or flagging of RFI. For those clusters, one or two additional rounds of self-calibration were performed after the pipeline.

Each 8 MHz channel was corrected for the VLA primary beam attenuation, and all channels were smoothed to a circular Gaussian restoring beam at the resolution of the lowest frequency channel, to ensure that all channels have the same angular resolution. This resulted typically in a synthesised beam size of 6–7 arcsec. The distribution of central root-mean-square (RMS) noise in the full-band Stokes I images is given in Fig. 1. Most targets have an RMS noise of around 20–30 μJy beam⁻¹ in the centre of the field. Two clusters, G033.46−48.43 and G226.17−21.91, have been removed from the sample. Calibration artefacts from a bright radio source were completely dominating the G033.46−48.43 field and during observations of G226.17−21.91 most of the data was lost to interference by a thunderstorm.

Fig. 1. Central RMS noise in Stokes I in the 124 observed target fields.

We found significant flux density variations between the spectral windows in all observations, also noted by Di Gennaro et al. (2021b), probably related to bandpass calibration or deconvolution uncertainties. To mitigate this problem as much as possible, we aligned the flux scale per observing run by fitting a simple power-law model to all bright Stokes I sources with at least a signal-to-noise ratio of 100,

$$\begin{array}{c}\hfill I_{\nu }=I_0{\nu }^{\alpha },\end{array}$$(1)

where α represents the spectral index and I_ν is the Stokes I intensity. Correction factors for each spectral window were determined per observing run by averaging the correction factors of individual sources. These correction factors were usually of the order of 5–10%. The corrections were applied to the Stokes I, Q and U fluxes.

The final 124 calibrated radio images are shown as a mosaic in Fig. 2. The five fields with RMS noise higher than 60 μJy beam⁻¹ in Fig. 1 are caused by calibration artefacts from bright sources at the edge of the fields, and in one case in the centre of the field. Direction-dependent calibration (e.g., Tasse 2014) could improve the quality of the images affected by bright off-axis sources, but these few fields should not significantly affect the results presented here. We decided to keep all 124 fields for our analysis because even in the five fields with bright artefacts 27 polarised radio sources were still relatively unaffected by those artefacts and could be used in the analysis.

Fig. 2. VLA 1–2 GHz primary beam corrected total intensity images of the 124 Planck clusters. The images are smoothed to the resolution of the lowest frequency channel (typically 6″) and the size is equal to the field-of-view at 2 GHz (0.35  ×  0.35 deg2). The colour scale is logarithmic with the scale range determined individually per cluster for visualisation purposes. The order of clusters follows the order in the Table E.1, in row-major order.

3. Methods

This section details the source finding of both polarised and unpolarised radio sources and the determination of their polarisation properties. Finally, we explain the optical counterpart identification and estimation of the redshift of the sources, such that we can classify the host galaxies as background sources, cluster members or foreground sources.

3.1. Polarised source finding

Linear polarisation can be expressed as a complex quantity by a combination of Stokes Q and Stokes U or written as a complex vector

$$\begin{array}{c}\hfill P({\lambda }^2)=Q+iU=p_{0}Iexp(2i\chi ),\end{array}$$(2)

where λ indicates the observed wavelength, p₀ the polarisation fraction, I refers to the Stokes I intensity and

$$\begin{array}{c}\hfill \chi ({\lambda }^2)=\frac{1}{2}arctan\left(\frac{U}{Q}\right)\end{array}$$(3)

is the polarisation angle. Faraday rotation introduces a wavelength-dependent rotation of the polarisation angle χ. In the general case, the Faraday depth of a source is defined as (Burn 1966; Brentjens & de Bruyn 2005)

$$\begin{array}{c}\hfill \varphi (\mathbf{r})=0.81{\displaystyle \int }{n}_{\mathrm{e}}\mathbf{B}·\mathrm{d}\mathbf{r}\left[\mathrm{rad}{\mathrm{m}}^{-2}\right],\end{array}$$(4)

where n_e is the electron density in parts per cm⁻³, B is the magnetic field in μGauss and dr the infinitesimal path length increment along the line of sight in parsecs. We adhere to the definition that ϕ(r) > 0 implies that the magnetic field is pointing towards the observer.

In the simplest case, where only one source is present along the line of sight without internal Faraday rotation, the Faraday depth ϕ is equal to the RM of a source, and the observed rotation can be expressed as

$$\begin{array}{c}\hfill \chi ({\lambda }^2)={\chi }_0+\varphi {\lambda }^2.\end{array}$$(5)

Faraday rotation may cause polarised sources to be undetected in the wide-band Stokes Q and U images or in the linearly polarised intensity ($\sqrt{Q^2+U^2}$) images. This is because the Stokes Q and U intensities can be both positive and negative, resulting in averaging out the frequency integrated signal if the RM is significant. To solve this problem, we used the Faraday RM-synthesis (Brentjens & de Bruyn 2005) technique. RM-synthesis aims to approximate the Faraday dispersion function F(ϕ) by Fourier inversion of the following equation

$$\begin{array}{c}\hfill P({\lambda }^2)={\displaystyle {\int }_{-\infty }^{+\infty }}F(\varphi )e^{2i\varphi {\lambda }^2}\mathrm{d}\varphi ,\end{array}$$(6)

where P(λ²) is the complex polarised surface brightness (Eq. (2)) as a function of the observing wavelength (squared) and ϕ is the Faraday depth of the source (Eq. (4)). Calculating the Faraday dispersion function F(ϕ) essentially corresponds to de-rotating polarisation vectors to their position at an arbitrary wavelength ${\lambda }_0^2$. However, we note that RM-synthesis only approximates the Faraday dispersion function because we cannot sample all wavelengths. The limitations of our frequency setup can be expressed with the three following quantities (Brentjens & de Bruyn 2005). The maximum Faraday depth to which we have more than 50% sensitivity is given by the channel width: δλ²

$$\begin{array}{c}\hfill ||{\varphi }_{max}||\approx \frac{\sqrt{3}}{\delta {\lambda }^2}\approx 1200\left[\mathrm{rad}{\mathrm{m}}^{-2}\right].\end{array}$$(7)

The resolution in ϕ space is determined by our wavelength coverage, with the full-width half-maximum given by

$$\begin{array}{c}\hfill \delta \varphi \approx \frac{2\sqrt{3}}{\mathrm{\Delta }{\lambda }^2}\approx 52\left[\mathrm{rad}{\mathrm{m}}^{-2}\right].\end{array}$$(8)

The maximum scale we can resolve in ϕ space (analogous to resolving-out extended radio sources in synthesis imaging) is given by the shortest observable wavelength

$$\begin{array}{c}\hfill \mathrm{maximum}\mathrm{scale}\approx \frac{\pi }{{\lambda }_{min}^2}\approx 140\left[\mathrm{rad}{\mathrm{m}}^{-2}\right].\end{array}$$(9)

Because the resolution in ϕ space is smaller than the maximum scale we can resolve, we are technically able to detect slightly extended sources in Faraday space (i.e. Faraday thick sources). Typical values of RM found in clusters are usually of the order of 10² rad m⁻², going up to 10³ rad m⁻² in dense cool-core clusters (e.g., Abell 780 and Cygnus A; Taylor et al. 1990; Sebokolodi et al. 2020). Thus with the current frequency setup, we are sensitive to the typical amount of Faraday rotation in clusters.

We performed RM-synthesis using the pyrmsynth² module, weighting by the inverse RMS noise of the channels and ignoring bad channels. The result is an ‘RM-cube’ with two spatial axes and a Faraday depth ϕ axis, that contains the polarised intensity at each pixel location as a function of the Faraday depth, sampled from ϕ = −2000 to ϕ = 2000 rad m⁻² in steps of 10 rad m⁻². The peak polarised intensity map is then made by taking the maximum value along the ϕ axis. The peak polarised intensity map for all clusters is shown as a mosaic in Fig. 3.

Fig. 3. Same as Fig. 2, but for peak polarisation intensity. The colour scale used here is arcsinh. We note that the calibration artefacts visible in cluster position 107 (zero-based row-major index Beck & Krause 2005; Bonafede et al. 2011) were not used in the analysis, but the field was kept as two polarised sources at the edge of the primary beam were relatively unaffected by the artefacts.

To find polarised source candidates automatically, we used the source finder program PyBDSF (Mohan & Rafferty 2015) on both the Stokes I images and the peak polarised intensity maps after RM-synthesis. We set the parameters thresh_pix = 5.0 and thresh_isl = 3.0, meaning that a five-sigma threshold was used for the source detection and a three-sigma threshold was used during the fitting of the total intensity source properties. The background noise was calculated over the image in a box with a size of 3 arcmin in steps of 1 arcmin to account for the varying background noise to the primary beam.

We found that PyBDSF performed better when de-correcting the peak polarised intensity maps for the primary beam, such that an approximately flat-noise image was used for the source finding. More involved methods for polarised source finding were considered (e.g., moment analysis; Farnes et al. 2018), but our simple method was found to be sufficient given the still relatively small data size which allowed for visual inspection of the polarised source candidates. All polarised source candidates were cross-matched with sources found in the Stokes I images and source candidates that lie inside the extent of the Stokes I source were retained as real polarised sources. We defined the extent of the sources in the Stokes I map as twice the full width at half maximum (FWHM) of the fitted Gaussian, which is empirically found to be a good estimate of source sizes (e.g., Hardcastle et al. 2019). This source finding process proved successful in most of the observations, but all fields were also visually inspected and some manual intervention was needed for rare cases, such as clear polarised source candidates that were positioned just outside of the extent of the source in the Stokes I map. In total, PyBDSF found 6 807 source candidates in Stokes I and 819 source candidates in polarisation over the 124 target fields.

3.2. Fractional polarisation measurement

To determine the polarisation properties such as the intrinsic polarisation fraction p₀ and the Faraday depth ϕ of polarised source candidates, we can model the polarised emission as

$$\begin{array}{c}\hfill P({\lambda }^2)=p_{0}Iexp[2i({\chi }_0+\varphi {\lambda }^2)].\end{array}$$(10)

However, if a source emits at different Faraday depths along the same line of sight it suffers from depolarisation due to the differential Faraday rotation causing the emission from the far side of the source to be rotated more than emission from the nearby side of the source. This internal depolarisation can be modelled as (see Sokoloff et al. 1998, for details)

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

where ${\mathrm{\Sigma }}_{\text{RM}}^2$ represents the amount of depolarisation. A similar effect happens because we observe the sources with a finite spatial resolution. If the magnetic field in an external Faraday screen (e.g., the ICM) changes on scales smaller than the restoring beam sources are partly depolarised by beam depolarisation. This is an external depolarisation effect and can be modelled as (see Sokoloff et al. 1998, for details)

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

where ${\sigma }_{\text{RM}}^2$ models the amount of depolarisation. Finally, if the polarisation angle rotates significantly in a single frequency channel, bandwidth depolarisation occurs. This limits the maximum observable RM, as is given in Eq. (7).

Distinguishing between internal and external depolarisation effects can be done by measuring the spectral index of the polarised emission at lower frequencies with high resolution because external depolarisation effects are stronger at low frequencies (Arshakian & Beck 2011). In reality, there are probably both internal and external Faraday effects at play and a combination of the models could be used to fit the data. However, for this study distinguishing exactly between polarisation mechanisms is not important, as we are only interested in the polarisation fraction trend. The internal depolarisation of radio sources should not affect the general trend and can be found from the depolarisation ratio of sources at cluster outskirts (see Sect. 5). Therefore we decided to fit only the external depolarisation model given by Eq. (12). We fitted this model to the Stokes Q and U channels simultaneously and the total intensity (Stokes I) spectrum was modelled as a simple power-law (Eq. (1)). Fitting the Stokes I, Q and U channels directly has the advantage that we can assume Gaussian likelihoods because these channels have Gaussian noise properties, unlike the polarised intensity maps, whose distribution is Ricean. We fitted the integrated Q and U flux densities of each polarised source candidate, where the integration was performed over the extent of the polarised source as defined in Sect. 3.1. This means that separate polarised components of the same physical source (e.g., two polarised lobes of a single radio galaxy) were treated as separate sources during fitting. The uncertainty in the integrated flux density per channel was calculated as

$$\begin{array}{c}\hfill {\sigma }_i=\sqrt{{({\sigma }_{\mathrm{rms}}\times \sqrt{N})}^2+{({\delta }_{\mathrm{cal}}\times f_i)}^2},\end{array}$$(13)

where N is the number of beams covered by the source and σ_rms the background RMS noise in the corresponding channel. The second term accounts for the flux density variations explained in Sect. 2 by assuming a δ_cal = 5% error on the measured flux densities per channel, denoted by f_i.

The fitting was done using a Monte Carlo Markov chain (MCMC) fitting code developed by Di Gennaro et al. (2021b) to sample the posterior probability. The following uniform priors were assumed:

$$\begin{array}{c}\hfill \{\begin{array}{c}{I}_0\sim \mathcal{U}(0,\infty )\hfill \\ \alpha \sim \mathcal{U}(-\infty ,\infty )\hfill \\ p_0\sim \mathcal{U}(0,1)\hfill \\ {\chi }_0\sim \mathcal{U}(0,\pi )\hfill \\ \varphi \sim \mathcal{U}(-2000,2000)\hfill \\ {\sigma }_{\mathrm{RM}}^2\sim \mathcal{U}(0,\infty ).\hfill \end{array}\end{array}$$(14)

The initial values for the parameters were found through a least-square fit, using the ϕ as obtained from the RM-synthesis method as the initial guess for the Faraday depth. The posterior was sampled with 200 walkers for 1000 steps and a burn-in period of 200 steps was removed from each chain. The one-sigma uncertainties on the best-fit parameters are given by the 16th and 84th percentile of the chain. An example of the results on a polarised source with a good signal-to-noise ratio is shown in Fig. 4.

Fig. 4. Results of the MCMC fitting of the external depolarisation model on the Stokes I, Q and U channels to a polarised source with good signal-to-noise ratio. Panel a shows the measured flux densities per channel in black and the best-fit model in blue. Panel b shows the posterior distribution of the model parameters visualised as one and two-dimensional projections in a corner plot. Panel c shows the inferred polarisation fraction from the data and the best-fit model together with the residuals.

To judge whether the model m, given by Eq. (12), is a good fit to the data points y_i (i.e. the Stokes Q and U flux densities), we inspected the normalised residuals,

$$\begin{array}{c}\hfill R_i=\frac{y_i-m({\lambda }^2,\theta )}{{\sigma }_i},\end{array}$$(15)

where σ_i is the uncertainty on the polarised flux. It is not possible to determine analytically the number of degrees of freedom k in the external depolarisation model because it is a non-linear model, which is why we are not able to determine the reduced chi-squared value. In Fig. 5 we plot the distribution of the sum of the squared residuals (i.e. the χ² value) of the best fitting external depolarisation model to each polarised component. Most polarised sources have around 84–89 data points (i.e. channels) after masking the bad channels. This would give 80–85 degrees of freedom if the model was linear with 4 parameters. For comparison, we show also the theoretical χ² distribution with 80 degrees of freedom. The main peak of the sum of the squared residuals shows good agreement with the theoretical χ² distribution, indicating that most sources have acceptable fits. There is however a long tail of large χ² values, mainly caused by bright Stokes I sources, where residual calibration artefacts are more noticeable. To automatically discard bad fits, we decided to cut all fits with a χ² value that is 5σ away from the theoretical distribution, indicated by the dashed line in the Figure. This cut removed 148 polarised components. A few of these components are possibly Faraday complex sources, for which the simple model with a single RM component is not sufficient. We note that there are no apparent correlations between the χ² parameter and the derived best-fitting parameters p₀, χ0, ϕ, ${\sigma }_{\text{RM}}^2$, or the projected radius to the cluster centre, so we are not biasing our analysis by removing these sources. Additionally, for sources with low signal-to-noise polarised emission, a good fit (according to the χ² parameter) can be found by artificially large values of σ_RM. For these sources, the best-fit σ_RM is basically unconstrained, with large error bars. Therefore we decided to also cut sources where the fractional uncertainty on the best-fit σ_RM is larger than two. This cuts 45 additional sources, so in total 193 polarised components have bad fits. These components are indicated in the polarised source Table C.1 by the column ‘Flagged’, where we have also flagged 11 components that are part of a radio relic.

Fig. 5. Distribution of the sum of the squared residuals (χ2) of the model fits to all polarised components. The plot is truncated at χ2 = 400 for visibility.

To calculate the upper limit on the fractional polarisation of sources detected only in Stokes I, we followed a method similar to Bonafede et al. (2011). We randomly sampled empty regions with a size of 10″, a bit larger than the synthesised beam. For these ‘noise sources’ we computed the polarised surface brightness and compared the distribution of the polarised surface brightness of the ‘noise sources’ to the distribution of real sources, taking into account the varying background noise level due to the primary beam of the VLA. We put the noise-dependent threshold of the surface brightness where noise sources constitute 10% of the real sources. The resulting threshold (P_t) as a function of the RMS noise is 0.04, 0.05, 0.06 and 0.08 mJy beam⁻¹ for sources with background RMS values 0–29, 29–36, 36–46 and 46+ μJy beam⁻¹ respectively, where the background noise bins are chosen such that in each bin there are an equal amount of simulated noise sources. For comparison, the threshold P_t calculated independently of the background noise level gives a value of P_t = 0.06 mJy beam⁻¹. The one-sigma upper limit on the fractional polarisation value is then calculated as

$$\begin{array}{c}\hfill F_p\le \sqrt{\frac{P_t^2}{I^2-{\sigma }_I^2}},\end{array}$$(16)

where I is the surface brightness of the unpolarised source and σ_I is the background RMS noise. This method gives conservative upper limits for extended unpolarised sources because the Stokes I surface brightness is computed over the entire extent of the source (rather than e.g., per lobe).

3.3. Optical counterparts

To determine the redshift of the radio sources, each radio source needs to be associated with an optical counterpart. PyBDSF is known to occasionally split up components of a single physical radio source (e.g., Williams et al. 2019). This particularly happens for large and extended sources, and often when the source has multiple disconnected patches of emission. To group PyBDSF Stokes I source candidates into single physical sources and to identify the optical counterpart, radio-optical overlays were created and every source was visually inspected. We used the g, r, z filters from the Legacy Survey (Dey et al. 2019) where available and used the Pan-STARRS survey (Chambers et al. 2016) for the 32 fields outside of the Legacy survey sky coverage.

As a first guess of the optical host galaxy, the nearest optical neighbour to the radio source was marked. This proved a good guess in 5806 out of 6807 total intensity source candidates. For the remaining 1001 sources, the best candidate optical counterpart position was manually marked from visual inspection of radio-optical overlays.

The source association was done in the same visual inspection step as the host galaxy identification. Out of the 6807 total intensity source candidates, 411 candidates were components of another source, leaving 6396 physical sources detected in total intensity. This indicates that PyBDSF in most cases correctly identified the total extent of the Stokes I source. We did not perform the source association step for the polarised components. Because different parts of a radio source can have different RM determinations and thus polarised intensities, we decided to treat separate polarised components as separate polarised sources, as for example was also done in Böhringer et al. (2016).

3.4. Redshift estimation

With the best-estimated location of the host galaxies determined, we employed different methods to estimate the source redshift. First, we checked whether a source has a spectroscopic redshift measurement available by cross-matching the host galaxy positions to the NASA/IPAC Extragalactic Database (NED³) and the Sloan Digital Sky Survey (SDSS DR16; Ahumada et al. 2020) with a matching radius of 0.5 and 3 arcsec respectively. If a spectroscopic redshift was found, but no uncertainty was quoted the redshift uncertainty is set to 0 in the catalogue.

If no spectroscopic redshift was found, a photometric redshift estimation was done from the Legacy Imaging Surveys Data Release 8 (Dey et al. 2019), which is the most sensitive optical survey covering the majority of our clusters. This approach is detailed fully in Duncan (2022) and provides high quality redshifts for galaxies at z < 1. For sources outside of the Legacy Survey, we calculated the photometric redshift using the Pan-STARRS grizy bands. We followed the method and used the code provided by Tarrío & Zarattini (2020), which estimates redshifts through local linear regression in a five-dimensional colour and magnitude space. The five-dimensional space consists of (r, g − r, r − i, i − z, z − y) where the letters indicate the extinction corrected Kron magnitudes (Kron 1980) of galaxies in the PanSTARRS grizy bands. The correction for interstellar extinction used the maps from Schlegel et al. (1998) and is described in detail in Sect. 2.3 of Tarrío & Zarattini (2020). To compute the photometric redshifts from the Pan-STARRS band, we found the 100 nearest neighbours in the five-dimensional space for each source, from a training set composed of 2 313 724 galaxies with spectroscopic redshifts, constructed by Tarrío & Zarattini (2020). The redshift was computed for all sources where at least four out of five colours were available. We note that missing features most often occur in very faint galaxies, which makes it likely that the sources are at a redshift z > 0.35, and are thus background sources. The quality of the Pan-STARRS photometric redshifts was checked by comparing to spectroscopic redshifts for sources where spectroscopic redshifts were available. Using standard literature metrics for the robust scatter σ_NMAD and outlier fraction OLF (cf. Dahlen et al. 2013; Duncan 2022) we find that the photometric redshifts have good quality, with σ_NMAD = 0.025 and OLF = 0.075.

The combination of all methods resulted in a redshift estimate for 77% (632/819) of the polarised sources and 67% (4544/6807) of the unpolarised sources. The distribution of redshifts estimates is given in Table 1. The final catalogues of polarised and unpolarised radio sources are provided with this paper, shown in Tables C.1 and D.1 respectively. These tables contain the polarised and unpolarised source properties, the best estimate for the redshift of the sources and the method used to get this estimate.

4. Magnetic field modelling

To compare observations with theoretical expectations, we modelled the magnetic field as a three-dimensional Gaussian random field, characterised by a single power-law spectrum. We followed the approach proposed by Tribble (1991), used in various works in the literature (e.g., Murgia et al. 2004; Guidetti et al. 2008; Bonafede et al. 2011, 2013; Vacca et al. 2012; Govoni et al. 2017; Stuardi et al. 2021). This approach starts with generating the vector potential of the magnetic field, A, in Fourier space, denoted by $\stackrel{\sim }{A}$. The amplitude and phase of the components of the vector potential were generated such that the phases are randomly distributed between 0 and 2π and the amplitudes follow a power-law given by

$$\begin{array}{c}\hfill |A_k{|}^2\propto k^{-\xi },\end{array}$$(17)

where k denotes the magnitude of the three-dimensional wave-vector k. The wave numbers k are related to the spatial scales Λ as $k=0.5·\frac{2\pi }{\mathrm{\Lambda }}$, where Λ refers to the reversal scale of the magnetic field, following the definition used by Murgia et al. (2004) Because the vector potential and magnetic field are real quantities, we made sure that the matrix $\stackrel{\sim }{A}$ is Hermitian (i.e. equal to its conjugate transpose). The components of the Fourier transform of the magnetic field are then given by the following cross product

$$\begin{array}{c}\hfill \stackrel{\sim }{B}(k)=ik\times \stackrel{\sim }{A}(k).\end{array}$$(18)

This results in the magnetic field B, which is simply calculated by (fast) Fourier transform, being divergence-free, isotropic and component-wise Gaussian random, with a power-law spectrum

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

where n = ξ − 2. A power-law spectral index of n = 3 implies that the magnetic field energy density is scale-invariant, for n < 3 the energy density is larger on smaller scales and for n > 3 the energy density is mostly in the larger scales (Murgia et al. 2004). The range of spatial scales Λ that can be explored is given by the size of the computational grid. The simulated maximum scale on which the magnetic field reverses is equal to Λ_max = π/k while the minimum scale Λ_min that can be probed is determined by the cell size.

The normalisation of the magnetic field was set after the Fourier transform such that the magnetic field strength approximately follows an assumed magnetic field profile. Like previous literature, we assumed that the magnetic field profile is proportional to the gas density profile, which is expected to happen during cluster formation from simulations (Dolag et al. 2008),

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

where B₀ is the average magnetic field strength at the cluster centre, n_e is the thermal electron gas density profile and η denotes the proportionality between the magnetic field strength and electron density. For η = 0.5, the magnetic field energy density is linearly proportional to the thermal gas density. The thermal electron density profile is available for every cluster in the Chandra-Planck sample, from the X-ray observations presented in Andrade-Santos et al. (2017), where the fitted profile was assumed to follow a modified double β model (see Vikhlinin et al. 2006 for more details):

$$\begin{array}{cc}& n_{\mathrm{e}}{n}_{\mathrm{p}}=n_0^2\frac{{(r/r_{\mathrm{c}})}^{-\alpha }}{{(1+r^2/r_{\mathrm{c}}^2)}^{3\beta -\alpha /2}}\frac{1}{{(1+r^{\gamma }/r_{\mathrm{s}}^{\gamma })}^{ϵ/\gamma }}\hfill \\ \hfill & +\frac{n_{02}^2}{{(1+r^2/r_{\mathrm{c}2}^2)}^{3{\beta }_2}}.\hfill \end{array}$$(21)

Given the modelled magnetic field and observed electron density profile, we calculated the expected Faraday rotation in the clusters by numerical integration of Eq. (4). Then, assuming an intrinsic polarisation p₀, and a single polarisation angle χ₀ for a polarised background screen, we computed the polarisation angle of the radio emission at 1.5 GHz at the cluster redshift (χ_obs) with Eq. (5). The predicted Stokes Q and U intensities were then obtained by inversion of $P=\sqrt{Q^2+U^2}$ and Eq. (3):

$$\begin{array}{c}\hfill Q=\pm \sqrt{\frac{p_0^2}{1+{tan}^2(2{\chi }_{\mathrm{obs}})}}U=\pm \sqrt{p_0^2-Q^2}.\end{array}$$(22)

Using the convention that Stokes Q is positive for −π/2 ≤ χ_obs ≤ π/2 and Stokes U is positive for 0 ≤ χ_obs < π. Finally, the images were convolved with a beam corresponding to a 6″ FWHM at the cluster redshift. From the convolved Stokes Q and U images, we calculated the expected depolarisation fraction at 1.5 GHz in the cluster rest-frame.

5. Results – observations

The intrinsic polarisation fraction of radio sources is often assumed to be the same for sources irrespective of their projected distance from the cluster centre. To test whether this is a good assumption, we plot in Fig. 6 the intrinsic polarisation fraction p₀ (see Eq. (12)) as a function of projected radius to the cluster centre. As the figure shows, there is a relatively large scatter in the intrinsic polarisation fraction of radio sources. This indicates that the assumption does not hold for this dataset and that the intrinsic polarisation fraction should be taken into account when estimating the amount of depolarisation.

Fig. 6. Best-fit intrinsic polarisation fraction against the projected distance to the cluster centre, normalised by the cluster R500. Polarised sources are coloured based on their position along the line-of-sight with respect to the nearest cluster, defined in Sect. 5.2.

To minimise the effect of the scatter introduced by source-dependent intrinsic polarisation, we calculated for every source the depolarisation ratio DP. We defined this as the ratio of the polarisation fraction at 1.5 GHz in the cluster rest-frame p_1.5 GHz to the intrinsic polarisation fraction p₀, using the best-fit model (Eq. (12)). In this way, we do not assume the same intrinsic polarisation fraction for all radio sources and take into account the cosmological redshift.

We combined the information from the upper limits (unpolarised sources) with the depolarisation ratio of polarised sources using the Kaplan-Meier (KM) estimator (Feigelson & Nelson 1985; see also Bonafede et al. 2011). The KM estimator is a non-parametric statistic used to estimate the complement of the cumulative distribution function, called the survival function. With x₁ < x₂ < ⋯ < x_r denoting distinct, ordered, observed values, the survival function is given by:

$$\begin{array}{c}\hfill S_{\mathrm{KM}}(x)=P(X\ge x)=1-F(x),\end{array}$$(23)

where F(x) denotes the cumulative distribution function of the random variable x, in our case the random variable is the depolarisation ratio measured in the centre of the band (i.e. DP). The KM estimator of the survival function is given by

$$\begin{array}{c}\hfill {\widehat{S}}_{\mathrm{KM}}(x)=\{\begin{array}{c}\begin{array}{cc}={\mathrm{\Pi }}_{i,x_i<x}{(1-\frac{d_i}{n_i})}^{{\delta }_i}\hfill & \mathrm{for}x>x_1\hfill \\ =1\hfill & \mathrm{for}x\le x_1,\hfill \end{array}\hfill \end{array}\end{array}$$(24)

with x_i the observed or censored depolarisation fraction of source i, d_i the number of sources with fractional polarisation equal to x_i, n_i the number of sources with (upper limits on) fractional polarisation ≥x_i and δ_i = (1, 0) if x_i is polarised or unpolarised, respectively. The KM estimator is here expressed in the case of a right-censored sample, and most algorithms indeed only support right-censored data. Thus, we transformed our left-censored data to right-censored data by subtracting the data from a constant, following Feigelson & Nelson (1985).

For unpolarised sources, we calculated upper limits on p_1.5 GHz as explained in Sect. 3.2, so an assumption on the intrinsic polarisation fraction must be made to translate this upper limit on the fractional polarisation to an upper limit on the depolarisation ratio. Thus, for these sources we calculated the depolarisation ratio assuming p₀ = 0.022, which is the median of the Kaplan-Meier estimator of all polarised radio sources detected at r > 1.5R₅₀₀. All KM estimates were calculated using the lifelines⁴ package (Davidson-Pilon 2019).

5.1. Full sample

Figure 7 shows in the left panel the depolarisation fraction for all sources in our sample, including upper limits that are below DP = 1. The right panel shows the median depolarisation fraction, calculated using the KM estimator by splitting the sample into bins of projected radius to the cluster centre. Each bin was chosen such that it contains an equal number of sources. The error bars reflect the 68% confidence interval of the KM estimator, added in quadrature with the uncertainty introduced by the fitting procedure. The uncertainty introduced by the fitting was estimated using a Monte Carlo method. For every source, we draw 1000 samples from a Gaussian distribution with a standard deviation equal to the one-sigma uncertainty on the depolarisation ratio given by the MCMC chain. We note that the error is dominated by the confidence interval of the KM estimator, thus that the effect of the uncertainty in the best-fit polarisation parameters is small. This means that we are limited by the number of polarised radio sources, and not by the quality of the data.

Fig. 7. Depolarisation against normalised projected radius. Left panel: full sample of sources and relevant upper limits. The error bars reflect the 68% confidence interval from the MCMC fitting procedure. Right panel: median of the Kaplan-Meier estimate of the depolarisation ratio survival function in different bins of projected radius to the cluster centre. The bin width is chosen such that each bin contains an equal number of polarised sources and is denoted by the horizontal lines (i.e. 0th and 100th percentile). The points are plotted at the median radius in each bin.

Figure 7 shows a clear trend of sources being more depolarised as they move towards the cluster centre, where the magnetic field strength and the line-of-sight column densities increase. The depolarisation ratio is around 0.92 beyond 2R₅₀₀, which is likely not an external, but an internal depolarisation effect because at these distances the column density and magnetic field strength of the intracluster medium would be too low to result in significant external depolarisation.

5.2. Background versus cluster members

To investigate whether there is a difference between depolarisation of cluster members and depolarisation of background radio sources, we classified each radio source according to the following definitions. We defined a source to be in front of a cluster if it lies at least 1σ_z away from a chosen boundary (Δz) around the cluster redshift

$$\begin{array}{c}\hfill z_{\mathrm{cluster}}-(z_{\mathrm{source}}+{\sigma }_{\mathrm{z}})>\mathrm{\Delta }z,\end{array}$$(25)

where z_cluster is the cluster redshift, z_source the source redshift and σ_zsource the one-sigma uncertainty on the source redshift. The values of z_source and σ_zsource are given in Tables C.1 and D.1 in the column ‘z_best’ for every source. Similarly, a source was defined to be behind a cluster if

$$\begin{array}{c}\hfill (z_{\mathrm{source}}-{\sigma }_{\mathrm{z}})-z_{\mathrm{cluster}}>\mathrm{\Delta }z.\end{array}$$(26)

All other sources were defined as inside the cluster. We have set Δz = 0.04(1 + z), following the definition of cluster membership used by Wen & Han (2015). Sources without an optical counterpart are likely faint sources at redshifts higher than z = 0.35, particularly because radio galaxies are often hosted by massive elliptical galaxies which should be easily detectable at z < 0.35 at the depth of Legacy and PanSTARRS. Therefore, sources without an optical counterpart were also defined as background sources. We verified, through a two-sample KS test, that the measured polarisation fraction of the sample of sources without an optical counterpart does not significantly differ from the sample of sources with an optical counterpart (p-value 0.14), implying that they pass similar Faraday screens.

We investigated the depolarisation effect as a function of radius and electron column density to partially split the degeneracy between magnetic field strength and electron column density, which are both a function of radius, and both increase the amount of depolarisation (Eq. (4)). To determine the electron column density for sources inside the clusters we integrated the best-fit electron density profile along the line-of-sight, from the centre of the cluster out to R₅₀₀, thus effectively integrating over half the sphere. For sources located behind the cluster, this column density was multiplied by two because we have assumed spherical symmetry in the electron profiles.

In Fig. 8 we show the depolarisation ratio calculated in different bins of normalised projected distance or electron column densities, for cluster members and background sources separately. Firstly, the figure shows the difficulty of detecting background sources close to the cluster centre, which means larger bins need to be used for background sources than for cluster members. For completeness, the full sample of data-points is plotted in Fig. B.1.

Fig. 8. Median of the Kaplan-Meier estimate of the depolarisation ratio survival function in different bins of radius and electron column density. Sources inside clusters are shown in green and sources behind clusters are shown in red. The bin width is chosen such that each bin contains an equal number of polarised sources and is denoted by the horizontal lines. The points are plotted at the median radius or column density in each bin.

Secondly, we detect a significant difference between cluster members and background sources in the highest column density bin (right panel). This could arise because at similar column densities background sources are projected further away from the cluster centre than cluster members and thus probe smaller magnetic field strengths, causing less depolarisation. Additionally, because cluster members are easier to detect near the cluster centre, the largest column density bin also samples preferentially higher column densities for cluster members.

Conversely, we expect that at similar radii, background sources probe higher column densities and are thus more depolarised. However, we do not significantly detect this difference given the uncertainties and the large bin size of background sources near the centre of the cluster.

Thirdly, at radii where we have similar sampling (i.e. r > 0.5R₅₀₀), we do not see a significant difference between cluster members and background sources. To statistically confirm this, we used the non-parametric log-rank test (Feigelson & Nelson 1985), used frequently with other astronomical works dealing with survival analysis (e.g., Bonafede et al. 2011; Kang et al. 2020; van Terwisga et al. 2022). The resulting survival curves of cluster members and background sources are shown in Fig. 9, and according to the log-rank test with p-value 0.89 we cannot reject the null hypothesis that there is no difference between background sources and cluster members.

Fig. 9. Survival function (i.e. 1-CDF) inferred from the Kaplan-Meier estimator for all sources located at r/R500 > 0.5, and r/R500 < 0.5. Background sources and cluster members are indicated by red and green, respectively. The grey dashed line shows the location of the 50th percentile, indicating the median for both populations.

Lastly, at the inner radii (r < 0.5R₅₀₀) where we do not have a similar sampling of background sources and cluster members, we see a hint of more depolarisation detected in the cluster member population. However, the log-rank test returns p = 0.13, indicating that with the current sampling, this result is not statistically significant.

5.3. Dynamical state

The magnetic field evolution of galaxy clusters remains poorly constrained. During the lifetime of clusters, mergers with other clusters or smaller substructures can alter the structure and strength of the magnetic field significantly. This section focuses on possible differences between merging and relaxed systems.

Generally, relaxed clusters show strongly peaked, symmetrical X-ray emission that has a radiative cooling time much shorter than the Hubble time (e.g., Fabian 1994). These clusters show the shortest cooling times in their cores and are therefore often referred to as cool-core clusters. Cluster mergers can destroy the cool core and significantly disturb the observed X-ray morphology (Burns et al. 2008). Thus, X-ray morphological parameters such as the concentration or cuspiness of the gas density profile can be used to determine whether a system has a cool core (e.g., Andrade-Santos et al. 2017).

We use the X-ray morphology parameters derived from the Chandra observations of 93 clusters in our sample in Andrade-Santos et al. (2017) to determine the presence or absence of a cool core. We note that this split does not perfectly correspond to a split in the dynamical state, as there are rare examples of merging clusters that still show a cool core (e.g., Somboonpanyakul et al. 2021). However, this split is sufficient to generally divide the sample into merging and relaxed systems. Using the concentration parameter calculated in the 0.15–1.0R₅₀₀ range (C_SB) by Andrade-Santos et al. (2017) to classify clusters as cool-core or non-cool-core, we found that 65% (60/93) of the clusters in our sample are non-cool-core (NCC) and 35% (33/93) are cool-core (CC) clusters.

Figure 10 shows the depolarisation effect separately for CC and NCC clusters in equal frequency bins, the full sample is plotted in Fig. B.2. We see a hint in the first radius bin of detecting more depolarisation in NCC clusters than in CC clusters.

Fig. 10. Median of the Kaplan-Meier estimate of the depolarisation ratio survival function in different bins of radius. Cool-core clusters are shown in blue and non-cool-core clusters are shown in orange. The left panel shows bin widths (denoted by horizontal lines) chosen such that each bin contains an equal number of sources detected in polarisation and the right panel shows manually selected bins. The points are plotted at the median radius in each bin.

To separate the effect of the central cooling core region in CC clusters, we have manually defined bins of projected radius in the right panel of Fig. 10. We have chosen an inner radius bin of 0.0−0.2R₅₀₀ because the effect of the cooling core is significant only in the inner ∼0.2R₅₀₀ of CC clusters (e.g., Vogt & Enßlin 2006; Eckert et al. 2011). The right panel shows that the larger depolarisation fraction in NCC is dominated by sources detected at r > 0.2R₅₀₀. In fact, sources detected at r < 0.2R₅₀₀ show a hint that there is more depolarisation in the central cooling core region of CC clusters, although the uncertainties are large due to the low number of sources detected near the centre of CC clusters. At r < 0.2R₂₀₀, we have detected only 9 sources and 16 upper limits in CC clusters, and 36 sources and 14 upper limits in NCC clusters. The significance of these results was determined by comparing the survival functions of sources detected in the 0.0−0.2R₅₀₀ and 0.2−1.0R₅₀₀ bins. The survival functions are shown in Fig. 11 and the log-rank test yields p-values of 0.22 and 0.001 for the 0.0−0.2R₅₀₀ and 0.2−1.0R₅₀₀ bins, respectively. This implies that the hint is statistically significant, with less depolarisation in CC clusters than in NCC clusters outside the core region. Conversely, inside the core region we do not have enough sources to significantly detect a difference between the two samples.

Fig. 11. Survival functions (i.e. 1-CDF) inferred from the Kaplan-Meier estimator for all sources in the 0.0−0.2R500 and 0.2−1.0R500 bins, separately for sources detected around CC and NCC clusters. The grey dashed line shows the location of the 50th percentile, indicating the median for both populations.

To examine to what extent the results of the CC/NCC split are affected by the position along the line-of-sight of the sources, we repeat this analysis separately for background sources and cluster members. Figure 12 shows that we do not detect a difference between NCC and CC clusters in either sub-sample. This is likely because of the low number of sources left in each sub-sample. We are mainly limited by the number of polarised sources detected near the centre of CC clusters. Comparing the survival curves of the NCC and CC sample for sources detected below 0.5R₅₀₀, the log-rank test returns p-values of 0.07 and 0.24, for the sources inside clusters and behind clusters, respectively. Thus, we cannot significantly detect differences between NCC and CC clusters when splitting the sample into background and cluster members due to the limited amount of data points. We do see that most of the depolarisation found at small radii is from cluster members, although the uncertainties become quite large due to the small sample sizes.

Fig. 12. Same as the left panel of Fig. 10, but separately for cluster members (left panel) and background sources (right panel).

5.4. Cluster mass and redshift

Although the sample of clusters is constrained to a relatively low redshift range (z < 0.35), we can attempt to trace the evolution of the magnetic field of clusters, by splitting the sample based on cluster redshift. We note that the redshift of the host cluster should correlate with the amount of beam depolarisation because the same telescope resolution corresponds to larger physical areas probed at higher redshifts. This means that we effectively average over larger magnetic field scales, and thus expect more beam depolarisation. Another effect that we have to take into account is the selection function of the Planck cluster sample. There is a strong correlation between cluster mass and redshift, with the most massive clusters preferentially being detected at high redshift (see Fig. 26 in Planck Collaboration XXVII 2016). This means that a cut in redshift effectively also corresponds to a mass-cut, as shown in Fig. 13. As the figure shows, it is not possible to separate the effects of cluster redshift and cluster mass because there is almost no overlap in the same mass range.

Fig. 13. Distribution of cluster mass in the low- and high-redshift samples.

Figure 14 shows the depolarisation trend for low- and high-redshift clusters separately, and the full sample is shown in Fig. B.3. The low-redshift sample contains clusters with z < 0.175 and the high-redshift sample contains clusters with 0.175 < z < 0.35. The first thing to note is that we detect significantly more sources at lower projected radii through the low-redshift clusters. Each projected radius bin has 78 detected sources through low-redshift clusters, while the high-redshift clusters have only 43 detected sources per bin. This is expected because the larger angular size of low-redshift clusters makes it easier to detect polarised sources, especially in the centre of the cluster. The low number of polarised sources detected at low projected radii in the high-redshift sample makes it difficult to compare the two populations. Therefore, we performed a bootstrap re-sampling to enforce that we have a similar sampling of radius in both low- and high-redshift clusters. This was repeated 1000 times, with one realisation of the re-sampled values shown in Fig. 15. Out of 1000 log-rank tests comparing the high-redshift sample to the sub-sampled low-redshift sample, 4% (43/1000) of the tests returned p < 0.05, indicating that we cannot distinguish a difference in depolarisation in the low- and high-redshift sample of clusters. However, this is likely due to the low number of sources in the high-redshift sample.

Fig. 14. Median of the Kaplan-Meier estimate of the depolarisation ratio survival function in bins of projected radii and column density, for the low- and high-redshift clusters. The bin width is chosen such that each bin contains an equal number of sources detected in polarisation and is denoted by the horizontal lines. The points are plotted at the median projected radius.

Fig. 15. Distribution of the projected radius of sources detected in the low and high-redshift cluster samples. The solid blue histogram shows the result of one realisation of randomly sub-sampling sources in the low-redshift cluster sample to match the distribution of sources in the high-redshift sample.

5.5. Presence of a radio halo

There is an apparent dichotomy in clusters regarding the presence of a radio halo, where clusters that show a radio halo are almost always found to be dynamically disturbed, while clusters without a radio halo are more relaxed (Cuciti et al. 2021). However, there are some cases of merging clusters without radio halos or with much fainter radio halos than usual (e.g., Cuciti et al. 2018; Russell et al. 2011). While these might be special cases, it is interesting to investigate whether there are differences between the magnetic fields in merging clusters that show a radio halo and those that do not.

We searched the literature for every cluster in our sample and found that out of the 60 clusters classified as merging, 26 have a radio halo detection, also incorporating the results of the second Data Release of the LOFAR Two-meter Sky Survey (Botteon et al. 2022). This thus splits the sample in about half, allowing for the same bins to be used. The resulting depolarisation curves are very similar, as shown in Fig. 16, and the log-rank test for similarity returned a p-value of 0.79. It is thus clear that with the current sample size we see no evidence of a difference in depolarisation between clusters with radio halos and clusters without radio halos.

Fig. 16. Comparison of the median observed depolarisation ratio against projected distance between the merging clusters with and without detected radio halos.

6. Results – modelling

This section details the results of the modelling, and the comparison of theory with observations. We simulated magnetic fields following the approach laid out in Sect. 4. We used the density profiles presented in Andrade-Santos et al. (2017), which were fitted to the modified double β model shown in Eq. (21). Profiles were only available for the 102 clusters from the ESZ catalogue, so we could not model the depolarisation in the 24 new clusters from the PSZ1 (Planck Collaboration XXXII 2015) and PSZ2 (Planck Collaboration XXVII 2016) catalogues.

6.1. Effect of density profiles

In previous works, when clusters were stacked often a mean profile was assumed (e.g., Murgia et al. 2004; Bonafede et al. 2011). We first investigated the effect of the different electron density profiles. We used a subset of six arbitrarily chosen clusters around the same redshift z = 0.1, such that we probe about the same physical scales. The modelled magnetic field parameters for this experiment are B₀ = 5.0 μG, n = 11/3 and η = 0.5 with a box-size of 1024³ pixels, where each pixel represents 2 kpc. The minimum magnetic field reversal scale Λ_min is thus 4 kpc, and the maximum reversal length scale, Λ_max = 1024 kpc. At the redshift of 0.1, the 6″ beam corresponds to a physical scale of 11 kpc. All models start from the same random initialisation of the magnetic field vector potential A, meaning that the only difference between the simulated clusters is the assumed electron density profile. The properties of the clusters are given in Table 2.

The resulting depolarisation ratio as a function of electron column density is shown in Fig. 17. From this figure it is clear that even at the same electron column densities, the amount of depolarisation can be quite different in different clusters, depending on the electron density profile. This is because at the same electron column densities, the local electron density profile along the line of sight can still differ quite a lot between clusters. This also influences the magnetic field strength along the line of sight because we assumed a relation between the magnetic field strength and the electron density profile. This means that it is important to take into account the different electron density profiles of the clusters, rather than define a mean electron density profile to stack the clusters.

Fig. 17. Modelled depolarisation ratio as a function electron column density for six different clusters around z = 0.1. The assumed parameters for the magnetic field are B0 = 5.0 μG, n = 11/3 and η = 0.5, with Λmin = 4 kpc and Λmax = 1024 kpc.

6.2. Effect of the magnetic field strength and fluctuation scales

We have chosen an arbitrary cluster, PLCKESZ G039.85−39.98, located at z = 0.176, to investigate the qualitative effect on the depolarisation profiles of changing the scales on which the magnetic field fluctuations and the central magnetic field strength B₀. The effect of increasing the magnetic field is easily understood to result in more depolarisation because the scatter in RM increases (e.g., Murgia et al. 2004). To understand the effects of changing the fluctuation scales, we must consider two different competing effects. First, as more power is put into larger scale fluctuations (i.e. increasing n), the scatter in RM over the entire cluster increases because one is integrating coherently over longer path lengths (cf. Eq. (4)). At the same time, because the fluctuations on smaller scales are reduced, the scatter in RM over the region probed by each individual observing beam decreases. Thus, depending on the size of the observing beam, this will either increase or decrease the amount of depolarisation as n changes.

We plot the modelled depolarisation profiles in Fig. 18 as a function of different parameters. As expected, the amount of depolarisation increases with increasing magnetic field strength. When the magnetic field energy density is mostly on large scales (i.e. n = 4), the depolarisation profile becomes quite flat as a function of projected radius because the magnetic field becomes correlated on scales larger than the beam. As the magnetic field becomes more correlated on smaller scales (i.e. from n = 4 to n = 2, green lines), the amount of depolarisation increases. However, as we put even more power on smaller scales (i.e. from n = 2 to n = 1) we reach the turn-over point where the effect of decreasing the RM scatter over the entire cluster dominates increasing the RM scatter on regions probed by the beam, resulting in less depolarisation. The exact turn-over point in the slope n depends on the size of the sampling region (i.e. the observing beam size) that the simulated radio images are smoothed with and can be different for different locations in the cluster, which have different magnetic field strengths.

Fig. 18. Effect of the magnetic field parameters on the depolarisation profile of the cluster G039.85−39.98, with R500 = 1.2 Mpc. The varied parameters are the central magnetic field strength B0 in units of μG, the magnetic field power-law index n and the maximum correlation scale Λmax in kpc. For all but one model, the maximum correlation scale was set to the maximum allowed by the image.

The degeneracy between a steep and strong magnetic field (e.g., B₀ = 10 μG, n = 4) and a shallower and weaker magnetic field (e.g., B₀ = 5 μG, n = 2) is also clear from this figure. This implies that using depolarisation alone makes it difficult to disentangle between a weaker magnetic field with a steep power-law index, or a shallower magnetic field with a flatter power-law index. The effect of setting a maximum fluctuation scale of Λ = 16 kpc does not strongly influence the depolarisation ratio except somewhat at the cluster outskirts. This can be explained by the fact that the observing resolution (FWHM of 18 kpc at the cluster redshift) is comparable to the maximum fluctuation scale. However, the amount of depolarisation does decrease slightly because the scatter in RM over the entire cluster will be smaller due to the magnetic field being less correlated along the line of sight.

6.3. Comparison with observations

To compare the models to the data, we modelled the depolarisation ratio as a function of projected radius for every cluster for which a Chandra observation and thus electron density profile was available. Given the computational intensity (which scales as N³ where N is the number of pixels) of generating many magnetic field cubes out to typical values of R₅₀₀, we decided to simulate different clusters with different resolution depending on the cluster redshift. To simulate the depolarisation effect, the model resolution should be at least a few times the physical resolution given by the synthesised beam of the radio observations. The resolution of the synthesised beam is given in Table 3 with the accompanying model resolution and model grid sizes used. For clusters below z = 0.05 it was not feasible to generate simulations, since the physical resolution (FWHM) of the 6″ synthesised beam corresponds to less than 6 kpc, and thus the resolution of the models should be higher than 1 kpc pixel⁻¹, which made the cube size unfeasibly big to generate. All modelled depolarisation profiles are shown in Fig. B.5 for completeness.

One final effect that we have to take into account when comparing the model to the data is internal depolarisation. Figure 8 showed that the depolarisation ratio is around 0.92 at the cluster outskirts. The fact that this is not 1.0 is likely caused by internal depolarisation effects. This internal depolarisation effect should not affect the trend, and therefore we multiply the simulated depolarisation ratio by the depolarisation ratio measured in the cluster outskirts to incorporate this effect.

6.3.1. Background versus cluster members

We can model the expected difference between the depolarisation of cluster members and background sources assuming that this difference can be fully attributed to the larger path length of background sources through the cluster. We assume that the radio emission from cluster members on average intersects about half the ICM column density and that emission from background sources travels through the full column. Theoretically, this is expected to give on average a factor of two larger Faraday depth for background sources (cf. Eq. (4)) and a factor $\sqrt{2}$ in σ_RM, which theoretically should not result in more than a factor two in depolarisation (cf. Eq. (12)) for the wavelength range that we are probing. Indeed, when modelling the depolarisation profiles occurring as a result of a Faraday screen halfway inside the cluster versus a Faraday screen behind the cluster, Fig. 19 shows that the location of the Faraday screen only has a marginal effect on the resulting depolarisation. This is in agreement with the results shown in Fig. 8, where no clear difference between background sources and cluster members was found.

Fig. 19. Average depolarisation ratio profile over all simulated clusters using a Faraday screen located behind (blue dashed line) or inside (orange solid line) the cluster. The uncertainty interval indicates the standard error on the mean of the simulated profiles.

6.3.2. Average magnetic field properties

As shown in Sect. 5.2, we did not detect a significant difference in the depolarisation of cluster members and sources located behind the clusters. This is consistent with a picture where only the difference in path length between background radio emission and radio emission from the cluster medium affects the depolarisation of radio sources. To estimate the average properties of magnetic fields in galaxy clusters we can therefore compare our models with the depolarisation calculated over all sources (background and cluster members) to maximise the signal-to-noise ratio.

Although η is an important parameter that can influence the magnetic field estimates and the dependence on the fluctuation power-law slope n (Johnson et al. 2020), we fixed η = 0.5 to reduce the number of free parameters. This value is chosen such that the magnetic field energy density follows the thermal gas density (as found in e.g., Coma and Abell 2382 Bonafede et al. 2010; Guidetti et al. 2008). We then varied B₀ = [1, 5, 10] μG and n = [1, 2, 3, 4]. The maximum and minimum correlation scales are also fixed to the minimum and maximum size allowed by the computational grid.

The observed depolarisation trend was re-calculated in five equal-width bins between 0−1R₅₀₀ using only sources detected in clusters that are part of the modelling, to make a fair comparison. The results of the comparison of the data with the modelled profiles are shown in Fig. 20. Because of the degeneracy between n and B₀ (shown in Sect. 6.1) and the fact that we are averaging over many individual clusters with different electron density profiles, there is a large overlap between the different models. Still, it is clear that B₀ = 1 μG does not fit the data for all values of n. For values of n between 1 and 4, the best fitting average central magnetic field strength is between 5−10, μG, but due to the degeneracy it is not possible to distinguish between these models.

Fig. 20. Comparison of the median observed depolarisation ratio in bins of projected radius, denoted by the horizontal lines, with the modelled depolarisation ratio profile for models with different magnetic field strengths (in μG) with η = 0.5. The model uncertainty interval indicates the standard error on the mean of the simulated profiles.

6.3.3. Dynamical state

Section 5.3 showed that NCC clusters appear to cause significantly more depolarisation than CC clusters outside 0.2R₅₀₀. To investigate in more detail to what degree this is caused by a difference in the magnetic field properties, we average the CC and NCC clusters separately. This allows us to quantify the effect of the different electron density profiles of the two cluster samples. If the thermal gas profiles are the only cause of the discrepancy in depolarisation between NCC and CC, then the same magnetic field parameters would fit both samples.

Figure 21 shows that we indeed expect more depolarisation outside the core region from NCC clusters than from CC clusters when they have the same magnetic field properties. This can be understood from the assumption that was made in Eq. (20), where the magnetic field energy density was assumed to follow the thermal gas density, normalised by the central electron density of the cluster. Because CC clusters generally have denser cores than NCC clusters, the magnetic field strength a few hundred kiloparsec away from the central cooling core declines faster than in NCC systems, where the denominator of Eq. (20) is smaller. Indeed the models also show that the amount of depolarisation increases more steeply towards the centre of cooling cores than in non-cool cores, which is in line with the observations shown in the right panel of Fig. 10.

Fig. 21. Comparison of the median observed depolarisation ratio with the modelled depolarisation ratio profile separately for the cool-core (CC) and non-cool-core (NCC) cluster sample. The modelled magnetic field parameters are B0 = 5 μG, η = 0.5 and n = 3 for both samples. The model uncertainty interval indicates the standard error on the mean of the simulated profiles.

6.3.4. Mass and redshift

Due to the low number of sources detected in the high-redshift sample, it was not possible to detect differences as a function of mass or redshift. Similar to the previous section, we can investigate to what extent we would expect a difference simply from the fact that we are probing a larger physical region at high redshift. However, in this case, the number of polarised sources detected in high-redshift clusters was already low due to the smaller angular size of the clusters and is even lower for the sample of clusters for which we also have density profiles available.

Within 1.0R₅₀₀, we have detected only 26 polarised sources in the high-redshift cluster sample, and 132 in the low-redshift sample with density profiles available. This causes large uncertainties for the high-redshift sample, particularly closest to the cluster centre. Figure 22 shows that, for similar magnetic field parameters, we would expect slightly more depolarisation from the high-redshift sample than the low-redshift sample, although again this effect is not strong enough to be observed in our data. We thus do not find evidence for a difference between the magnetic field properties of the high-redshift, high-mass and low-redshift, low-mass sample of clusters.

Fig. 22. Comparison of the median observed depolarisation ratio with the modelled depolarisation ratio profile separately for the low- and high-redshift cluster sample. The modelled magnetic field parameters are B0 = 5 μG, η = 0.5 and n = 3 for both samples. The model uncertainty interval indicates the standard error on the mean of the simulated profiles.

7. Discussion

We investigated the magnetic field properties in a sample of galaxy clusters through the effect of beam depolarisation. We confirm the hint of the depolarisation trend with cluster projected radius seen in previous studies (Bonafede et al. 2011; Stuardi et al. 2020) with a highly statistically significant result, as shown in Fig. 7. In this section, we discuss the implications of the results and the possible limitations of this study.

7.1. Cluster members versus background sources

One of the main questions that this paper addressed is whether there is a difference between using cluster radio galaxies and background radio sources to probe the magnetic fields in galaxy clusters. One could expect such a difference because cluster radio galaxies might locally reshape the magnetic field and density structure, causing a bias in the RM and amount of depolarisation. Interactions with surrounding gas have been suggested to affect observed RM distributions in various powerful radio sources (e.g., 3C75, 3C465, 3C270 and 3C353; Rudnick & Blundell 2003; Guidetti et al. 2011, 2012). However, this has not yet been shown in a statistical study.

Figure 8 and the log-rank comparison of the survival curves shown in Fig. 9 demonstrated that, although near the cluster centres (i.e. r < 0.5R₅₀₀) it is difficult to have similar sampling of cluster members and background sources, at radii where we have similar sampling (i.e. r > 0.5R₅₀₀) the depolarisation of cluster members and background radio sources is similar.

In Sect. 6.3.1 we modelled the difference between the amount of depolarisation expected for cluster members and background sources based on the different locations of the Faraday screens. We showed that this difference is minimal, and this difference is in line with the observed depolarisation trend for background sources and cluster members. This implies that there is no significant difference between using the depolarisation properties of cluster members or background sources as a probe of the cluster magnetic field.

Our results are in line with the findings by Bonafede et al. (2011) that used source angular size as a proxy of cluster membership and by Ensslin et al. (2003) that found that the biases from cluster members are not statistically significant. We note the caveat that the central region is still not well constrained with background sources.

When splitting the sample into NCC and CC clusters, there does seem to be a hint that there is a difference between background sources and cluster members near the cluster centres (Fig. 12), although only a few sources were detected near the central regions in these splits. A possible explanation for this is that there might be a pronounced effect on the cluster ICM from a select number of powerful cluster radio galaxies, which is averaged out when using a larger sample of sources. This means that when only a few cluster members are used to probe the magnetic field strength, the results may still be biased.

We thus did not find any strong differences between the depolarisation of cluster members and background sources in the full sample. However, larger samples might be able to pick up more subtle effects.

7.2. Magnetic field parameters

The average magnetic field properties of the cluster sample were explored by combining the depolarisation of all detected sources, irrespective of their redshift. The results in Fig. 20 showed that for all power spectrum indices, a central magnetic field strength higher than B₀ = 1 μG is needed to explain the observed depolarisation trend. For models with power-law indices between n = 1 and n = 4, an average central magnetic field strength between 5 and 10 μG proved to be the best fit, although it was not possible to distinguish between these models. Our results agree with previous radio observations that have shown that clusters have central magnetic field strengths between 1 and 10 μG with power spectrum indices between n = 2 and n = 4 (Murgia et al. 2004; Govoni et al. 2006, 2017; Guidetti et al. 2008, 2010; Laing et al. 2008; Bonafede et al. 2010; Vacca et al. 2012, 2018) and values from magneto-hydrodynamic simulations of clusters (Domínguez-Fernández et al. 2019).

With larger cluster samples or deeper cluster surveys with polarisation information, such as the MeerKAT Galaxy Cluster Legacy Survey (MGCLS; Knowles et al. 2022), it might be possible to group clusters with similar density profiles together. This would reduce the scatter in the modelled depolarisation trend and allow for a more accurate determination of magnetic field parameters.

7.3. Cluster properties

We investigated possible differences in observed depolarisation as a function of various cluster properties, such as whether a cluster is undergoing a merger. The magnetic field might be altered by cluster mergers, during which a massive amount of energy is released (up to 10⁶⁴ ergs on a few Gigayear timescales; Sarazin 2002). It is expected that this energy is injected on large spatial scales and released to smaller and smaller scales through turbulent cascades (Vacca et al. 2018; Domínguez-Fernández et al. 2019). Observations find central magnetic field strengths of around ∼1 μG and fluctuation scales up to a few hundreds of kpc in merging systems, while relaxed systems show higher central field strengths of around ∼10 μG and much smaller fluctuation scales (less than a few tens of kpc) (Taylor et al. 2002; Clarke et al. 2001; Vacca et al. 2018). This implies that, theoretically, we would expect a stronger depolarisation effect in CC clusters.

We investigated whether there are differences in the depolarisation found in CC and NCC clusters in Fig. 10. Surprisingly, we found that NCC clusters show more depolarisation than CC clusters outside the cooling-core region defined as r > 0.2R₅₀₀ (Fig. 11). When modelling (Fig. 21), it was found that the same central magnetic field strength in CC clusters results in less depolarisation outside the core than in NCC clusters because the magnetic field was assumed to scale with the electron density normalised by the central electron density, which is generally higher in CC clusters than in NCC clusters. Hence, the observed differences could be explained by the same magnetic field parameters in the CC and NCC sample.

To investigate to what extent this result is dependent on the cluster classification method, we also checked different morphological parameters, splitting the sample using the cuspiness and central gas density parameters from Andrade-Santos et al. (2017). In these splits, NCC clusters still showed more depolarisation than CC clusters outside the core region. We note that we could not use the entire sample of clusters in this analysis because only 93 out of 124 clusters are observed with Chandra in Andrade-Santos et al. (2017). A literature search resulted in the dynamical states for 9 more galaxy clusters, which also does not change the observed depolarisation trends significantly.

Thus, we found no strong evidence that CC clusters have significantly higher magnetic field strengths or smaller fluctuations scales than NCC clusters in the central regions, although the uncertainties were large as shown in Fig. 10. However, there is a hint that CC clusters indeed show more depolarisation inside the core region, as also found tentatively in Bonafede et al. (2011). Unfortunately, the typical size of the cooling cores in galaxy clusters is only 50–100 kpc, which is a region that is still poorly constrained in this study. The potential difference between the depolarisation in CC clusters and NCC clusters both inside and outside the core region should be investigated further because the sample size is still relatively small when splitting into multiple bins.

We also checked whether there is a correlation between magnetic field parameters and cluster mass or redshift. A positive correlation between magnetic field strength and cluster mass might be expected, as the observed radio power of giant radio halos is found to correlate with cluster mass, which can be reproduced by turbulent re-acceleration models with a positive scaling of the magnetic field strength with the cluster mass (Cassano et al. 2006a). However, the number of polarised sources detected in the high-redshift and high-mass cluster sample was too low to investigate a trend or differences between the low-mass and high-mass samples. A deeper survey such as MGCLS might be able to overcome this problem, although the sample of clusters should be large enough or carefully selected to break the redshift-mass selection bias discussed in Sect. 5.4.

Finally, we checked whether the presence of a radio halo in merging systems influences the observed depolarisation trend. Models based on the turbulent re-acceleration scenario usually define a radio halo as observable if the break frequency of the radio halo spectrum is above the observing frequency. In these models, the break frequency of the spectrum depends on the magnetic field strength, the cluster mass and the merging state (Cassano et al. 2006b). To investigate whether the magnetic field properties of clusters with a radio halo are different, we split the merging cluster sample based on the detection of a radio halo. No significant differences were observed between the depolarisation of clusters with radio halos and without radio halos, suggesting that they have similar magnetic field parameters. We checked that the cluster mass and redshift distributions are similar, (with a KS-test resulting in p-values of 0.2 and 0.8, respectively,) so the results are not biased by this. While these clusters are all classified as merging clusters according to the X-ray morphological parameters, a more in-depth study of their merging state might reveal that the clusters without radio halos are only minor merging systems where less turbulence is generated than in major mergers, which would be in line with the findings by Cuciti et al. (2021).

7.4. Possible caveats

In this section, we focus on the possible shortcomings of this work. Firstly, a single component external depolarisation model was used to fit the data. In reality, multiple interfering RM components can produce behaviour that is not proportional to λ² and even cause re-polarisation with decreasing frequencies (e.g., Pasetto et al. 2016). Observations of bright polarised sources observed at two different GHz frequencies have found that more than 25% of sources can show re-polarisation behaviour (Lamee et al. 2016). We have found in Sect. 3.2 that about 25% (193/819) polarised sources detected in this work are not well-fitted by the single component external depolarisation model. While fitting these sources with more complicated models (e.g., Brown et al. 2019) is beyond the scope of this work, we can briefly investigate how many sources show evidence of re-polarisation by allowing σ_RM to take negative values. This test resulted in 61 sources out of 819 that show a better fit with negative values of σ_RM. However, most of these sources do not show strong evidence for re-polarisation and could be fit almost equally well with a value of σ_RM that is around 0, and as such do not change the resulting depolarisation curve significantly. Additionally, the resulting median depolarisation as a function of radius is similar when incorporating the sources with bad fits, which reinforces the fact that we are not biasing the results by omitting these sources. The difference between the number of re-polarising sources found here and in the literature could be caused by the fact we measure the polarisation over the entire bandwidth, where almost always some depolarisation occurs, rather than at only two points in frequency points where multiple components might interfere and show re-polarisation.

Secondly, to derive upper limits an assumption on the intrinsic polarisation p₀ had to be made. For these sources we assumed p₀ = 0.022, which was the median intrinsic polarisation fraction of sources detected at r > 1.5R₅₀₀. If this assumption was too high, we are biasing the results through the inclusion of the upper limits by calculating too much depolarisation. However, not including upper limits would cause a bias in the opposite direction by omitting preferentially the most depolarised sources. While it is impossible to determine the intrinsic polarisation value for unpolarised sources, we can show that the depolarisation curve is not dominated by upper limits. Because the upper limits were computed conservatively over the entire extent of the total intensity sources, only 196 out of the 6807 unpolarised sources have an upper limit on the 1.5 GHz polarisation fraction that is below the assumed p₀ = 0.022. Thus, in terms of the number of sources, the upper limits are not dominating the results. The depolarisation trend with radius when omitting upper limits entirely is shown in Fig. 23. Cluster members now show less depolarisation in the centre of the cluster because the most constraining upper limits on the depolarisation fraction are found near the cluster centres, where the brightest sources are detected. However, it is clear that even without the inclusion of the upper limits, the depolarisation trend with radius is still clearly detected.

Fig. 23. Median depolarisation ratio in different bins of radius without incorporating upper limits from unpolarised sources. Sources inside clusters are shown in green and sources behind clusters are shown in red. Compared to Fig. 8, there is significantly less depolarisation in the lowest radius bin.

Thirdly, for the upper limits, we computed a polarised flux threshold as described in Sect. 3.2, which was dependent on the varying background noise level. This introduces a bias because all clusters are observed approximately in the pointing centre, so the upper limits are generally higher at larger projected radii. This means we are underestimating the amount of depolarisation more strongly at the edges of the field. Section A.2 investigates trends with angular distance from the pointing centre in detail and shows that this bias is very small compared to the observed depolarisation trend. Because the clusters are all observed near the pointing centre, other trends with projected distance from the cluster centre could also (partially) be due to instrumental or observational trends with angular distance from the pointing centre. These biases are also investigated in detail in Appendix A, where we present that there are indeed sources of bias, but through a Monte Carlo experiment we show that the effect of these biases is minimal compared to the observed depolarisation effect.

Fourthly, the effect of off-axis polarisation leakage can also mimic depolarisation because of the frequency dependence of the primary beam at a fixed angular distance from the pointing centre. This effect is expected to be at the order of the 1% level for VLA L-band observations (Jagannathan et al. 2017). To correct for this effect, full direction-dependent primary beam corrections need to be made (a-term corrections), which is possible for example with IDG (van der Tol et al. 2018), but is computationally expensive for large sample sizes and beyond the scope of this paper. However, the leakage effect is in the opposite direction from the observed trend because polarisation leakage effects are stronger near the periphery of the fields, while the observed depolarisation effect is stronger near the centre of the field. Additionally, we can see in Figs. 6 and A.1c that there is no significant increase in the measured intrinsic polarisation fraction of radio sources as a function of angular separation to the pointing centre. This implies that off-axis leakage effects are negligible for this study.

Lastly, electron density profiles were not available for all clusters studied in this work, with the 24 new clusters from the PSZ1 and PSZ2 catalogues not having Chandra observations. However, all clusters have been observed through the Sunyaev-Zel’dovich effect, which probes the integrated pressure along the line of sight. It has been shown that the pressure profile of galaxy clusters follows a relatively universal shape, called the universal pressure profile (UPP; Nagai et al. 2007; Arnaud et al. 2010). This profile scales in terms of the cluster properties M₅₀₀, R₅₀₀ and P₅₀₀, where P₅₀₀ is the characteristic pressure at an overdensity of 500. With an assumption on the cluster temperature, we can thus calculate a general electron density profile from the UPP to include these 24 clusters in our modelling. To derive the electron density profiles, we use the best-fit parameters for the UPP fit on Planck ESZ clusters from Planck Collaboration V (2013) combined with the best-fit mass-temperature relation on ESZ clusters from Lovisari et al. (2020). Including these 24 additional clusters in our modelling results in slightly more average depolarisation as a function of radius, but the final results do not change significantly, as shown in Fig. 24. Thus, assuming a universal pressure profile might be useful for future studies of larger, or higher redshift samples of clusters where X-ray observations are not available.

Fig. 24. Comparison of the median observed depolarisation ratio with the modelled depolarisation ratio profile using only clusters that have X-ray observations (blue) and all clusters (green) by calculating electron density profiles from a universal pressure profile (orange). The model uncertainty interval indicates the standard error on the mean of the simulated profiles.

8. Conclusion

In this work, we have utilised VLA L-band polarisation observations of a sample of 124 clusters from the Chandra-Planck Legacy Program for Massive Clusters of Galaxies to measure the depolarisation properties of radio sources inside and behind clusters. The main aims of the paper were to use the depolarisation ratio to (i) determine the average magnetic field properties in clusters, (ii) investigate whether there is a difference between using cluster members and background sources as probes and (iii) quantify the dependence of the magnetic field with cluster properties such as mass and dynamical state. We compared the data with modelled depolarisation trends by assuming the magnetic field is a Gaussian random field that follows the thermal electron density profile of the cluster. For the first time in a statistical polarisation study, we took into account the individual electron density profiles of different clusters when modelling the depolarisation ratio. We showed that the depolarisation ratio is a good probe of the magnetic fields in galaxy clusters. Our main results can be summarised as follows:

1. We clearly detect a trend of radio sources becoming more depolarised as they move (in projection) towards the cluster centre. This trend can be explained by models with a central magnetic field strength of 5−10 μG with power-law indices between n = 1 and n = 4 and cannot be easily attributed to observational or other systematic biases in the analysis.

2. The individual thermal electron density profiles of the clusters should be taken into account when modelling multiple clusters, as the theoretical depolarisation in separate clusters can be significantly different, even at similar electron column densities. This scatter might be overcome with larger samples when clusters with similar density profiles are grouped together.

3. The relation between the simulated beam depolarisation and fluctuation scale spectral slope is not monotonic. We found that simulated beam depolarisation can increase or decrease with increasing fluctuation scale spectral slope n depending on the size of the observing beam and the location in the cluster.

4. We found no statistically significant difference between the depolarisation properties of background and cluster sources, although background sources were rare to detect near the cluster centre, where cluster members were most often detected. The fact that we see no strong difference implies that the interaction between the radio sources in clusters and their local surrounding medium generally does not strongly influence their polarisation properties. Thus, in statistical studies, both in-cluster and background sources can be used as a probe of the magnetic fields.

5. Disturbed (non-cool-core) clusters showed more depolarisation in the 0.2−1.0R₅₀₀ region than cool-core clusters. After modelling, this effect was not strong enough to warrant different magnetic field parameters for disturbed or relaxed systems. While literature suggests that cool-core clusters have stronger magnetic fields inside the core region and should thus show more depolarisation, we did not significantly detect different polarisation fractions inside 0.2R₅₀₀ in cool-core and non-cool-core clusters. However, the uncertainties were large due to the low number of sources, and the most central (∼100 kpc) cooling core region is even more unconstrained in this study and should be investigated further.

6. The observed depolarisation in merging clusters that show a radio halo and merging galaxy clusters that do not show a radio halo is similar. This implies that the presence or absence of a radio halo in merging clusters is likely not dominated by the cluster magnetic field properties.

The biggest limitation of the study of magnetic fields in clusters through depolarisation is currently the number of polarised sources that are detected. With deeper cluster surveys and the advent of the SKA, depolarisation of radio sources will be a promising tool to study cluster magnetic fields.

Acknowledgments

The authors would like to thank the anonymous referee for the comments and suggestions that improved the quality of the 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). LR acknowledges partial support from U.S. National Science Foundation grant AST17-14205 to the University of Minnesota. AB acknowledges support from the ERC Starting Grant ‘DRANOEL’, number 714245 and from the from Italian MIUR grant FARE ‘SMS’. KJD acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 892117 (HIZRAD). Basic research in radio astronomy at the Naval Research Laboratory is supported by 6.1 Base funding. This work was performed using the compute resources from the Academic Leiden Interdisciplinary Cluster Environment (ALICE) provided by Leiden University. 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).

References

Appendix A: Possible biases

This appendix aims to investigate possible biases that might affect the analysis. The polarisation properties of radio sources might be influenced by other source properties such as size, total intensity, and whether the source is a single or multi-component source. If, for example, multi-component sources are easier to detect near the pointing centre and have different polarisation properties than single-component sources, the observed depolarisation effect with projected radius might be biased. This is because all clusters have been observed approximately in the pointing centre. As mentioned in Section 7.4, this also introduces a bias through the inclusion of the upper limits because the upper limits on the fractional polarisation are generally higher at larger projected radii due to the primary beam response. In the first section of this appendix, we check whether biases because of source properties are present, and in the second section we quantify the biases through a Monte Carlo experiment.

A.1. Source properties

To investigate whether there is a dependence between observed source properties and polarisation, taking into account the distance to the pointing centre, we plot running medians of various quantities versus the angular distance to the pointing centre. The uncertainty, σ_±, on the running median M is calculated as in Lamee et al. (2016),

$$\begin{array}{c}\hfill {\sigma }_{\pm }=|M-[p_{16},p_{84}]/\sqrt{N},\end{array}$$(A.1)

where p_x denotes the x-th percentile of the distribution, and N is the number of points in a bin. The amount of correlation is quantified by the Pearson (Pearson 1895) and Spearman (Spearman 1904) coefficients, shown in Table A.1. We define weak correlation for values of Pearson |r|≤0.3, moderate correlation for 0.3 < |r|≤0.7 and strong correlation for r > 0.7, using the common cutoff p-value of 0.05 for statistical significance. Because the Pearson coefficient only measures the linear relationship we also report the Spearman coefficient. The monotonicity of the correlation is given by the Spearman coefficient, where an absolute value of 1 indicates a perfectly monotonic relationship.

We plot in Figure A.1a the running median of the observed polarised source major axis versus the angular distance to the pointing centre. As the figure shows, there is some dependence of polarised source size on angular separation. Cluster members show an increase in median source size around 5 arcminutes from the pointing centre, while for background sources there is no significant correlation. The median source size is larger for cluster members than for background sources, which is expected simply because they are nearby sources.

Fig. A.1. Various quantities against angular separation separately for cluster members and background sources. The running median is shown with uncertainties in the shaded region. Histograms show the projected distributions along the axes.

The total flux density versus angular separation is shown in Figure A.1b. The background population clearly shows the effect of the primary beam response. Cluster members are less affected by this, possibly because they are all low redshift (z < 0.35) sources for which we are already sensitive enough to probe the majority of the cluster population. There is an excess of bright cluster members in the centre of the image, indicated by the peak at low angular separation. This means that we are detecting the brightest cluster members preferentially in the centre of the images. This is not unexpected because the centres of the clusters lie near the centres of the images, but it might bias the results if the total flux density is correlated with polarisation properties.

Although we are mainly interested in the depolarisation of radio sources, it is important to investigate the best-fit intrinsic polarisation fraction as a function of angular distance. Because of the chosen model and finite amount of bandwidth, there is a degeneracy between p₀ and σ_RM (see Fig. 4). While this degeneracy is simply a result of the fitting, real correlations between p₀ and σ_RM have been claimed before in the literature (e.g. Lamee et al. 2016). It is therefore important to inspect trends of p₀ with angular separation because that could create biases in the depolarisation trend. In Figure A.1c we plot the best-fit intrinsic polarisation fraction. Here we see that there is a weak correlation with angular distance for cluster members and no significant correlation for background sources. The fact that the median p₀ is lower near the centre of the image is likely (at least partly) caused by the fact that the brightest cluster members often lie in the centre of the cluster (Fig.A.1b), where it is thus possible to detect smaller polarisation fractions.

Lastly, there might be a difference in polarisation properties of single and multi-component sources, so we also plot this separation in Figure A.1d. We note that multi-component sources are often, but not always, cluster members. Both populations show no evidence for a strong correlation of the intrinsic polarisation with angular separation, and both populations have similar distributions of intrinsic polarisation, so it is unlikely that this is biasing the results significantly.

Now that we have established that the cluster population has a higher median flux density in the centre of the images and that the intrinsic polarisation fraction of cluster members is generally lower in the centre of the images, it is important to know whether these variables correlate with the depolarisation parameter σ_RM. Figures A.2 and A.3 show the trend of total flux density and intrinsic polarisation fraction with σ_RM. We see for background sources no clear correlation between total flux density and depolarisation. Cluster members do show that the brightest sources show more depolarisation. However, this is to be expected in the case of a magnetised depolarising intracluster medium (ICM) in the cluster centre if the brightest cluster members are also found preferentially in the cluster centre, which is indeed the case as shown by Figure A.1b. The question remains how much of this effect is a real effect and how much is caused by biases such as only picking up the most depolarised sources near the centre of the cluster. This is addressed in the next section. Figure A.3 shows the degeneracy between p₀ and σ_RM, particularly for large values of σ_RM. It is interesting that this trend implies that more depolarised sources have larger p₀, while Figure A.1c showed that smaller values of p₀ are generally found more towards the cluster centre. This trend would therefore cause a bias in the direction opposite to the trend expected from a depolarising ICM.

Fig. A.2. Same as Figure A.1, but for total Stokes I flux density against depolarisation parameter σRM.

Fig. A.3. Same as Figure A.1, but for intrinsic polarisation fraction against depolarisation parameter σRM.

A.2. Monte Carlo experiment

The previous section showed that there are no strong trends detected between source properties and polarisation properties or angular radius to the pointing centre, but there are weak trends in the data that possibly bias the results. Most notably, we have seen that the brightest cluster members are preferentially detected in the centre of the images and it is this class of sources that shows the most depolarisation. To quantify the bias introduced by selection effects and the fitting procedure, we took a Monte Carlo approach, simulating polarised radio sources with random properties that are taken from distributions that are representative of the data. If, through the effect of the choices made during the analysis or because of the radio source properties such as size and flux density a bias is introduced in the depolarisation curve, we should find that bias when employing the same methods on a sample of completely randomly (de)polarised sources. The distributions used to generate random polarised sources for this experiment are shown in Figure A.4, where we have fit gamma distributions to the strictly positive values p₀ and σ_RM and used a Cauchy distribution for RM to account for the large peak around RM = 0. The initial polarisation angle χ₀ was simply drawn from a uniform distribution between 0 and π. The steps taken in the Monte Carlo experiment were as follows. First, we calculated the total flux density per channel of all detected total intensity sources. Then, for each total intensity source, we randomly drew a value of p₀, RM, σ_RM and χ₀ from the representative probability distribution functions. Third, we computed simulated Stokes Q and U emission using the external depolarisation model given in Eq. 12 with the randomly drawn parameters. Then, we determined which sources are detected in polarisation at a 5σ level given the varying background noise due to the primary beam response. Finally, we fitted the detected sources with the MCMC IQU fitting code and calculated upper limits on the undetected sources as explained in Section 3.2. The resulting best-fit parameters and upper limits were again used to find the median depolarisation in bins of projected radius.

Fig. A.4. Probability distribution functions used to generate polarised radio sources during the Monte Carlo analysis shown in orange. The distribution of the real data is shown in blue.

This approach resulted in 1050 simulated sources detected in polarisation. The resulting median depolarisation trend with projected distance to the cluster centre is shown in Figure A.5 for simulated sources. Figure A.5 shows that there is indeed a very small bias from selection effects or fitting, with a minor trend showing slightly more depolarisation near the cluster centres than at the cluster outskirts. However, this trend is only a small fraction of the real detected trend in Figure 8. Thus, even though there are significant correlations as shown in Table A.1, they cause only a minimal bias because they have small correlation coefficients. This means that if the observed sources were a population with random polarisation parameters, there would not be a strongly detected depolarisation trend with radius. The striking difference between Figure A.5 and the depolarisation trend detected in the real data shown in Figure 8 makes a strong case that the observed depolarisation trend cannot be explained only by selection effects or biases.

Fig. A.5. Median of the Kaplan-Meier estimate of the depolarisation ratio survival function in different bins of radius for simulated data. The bin width is chosen such that each bin contains an equal number of polarised sources and is denoted by the horizontal lines. The points are plotted at the median radius in each bin

Appendix B: Full sample plots

For completeness, we plot in Figures B.1 to B.3 the full sample of data-points that have been summarised with the KM estimator in Figures 8 to 14. Figure B.4 shows an alternative visualisation of the full sample of data points, which also clearly shows the trend of sources becoming more depolarised as the projected radius and column densities increase. Finally, Figure B.5 shows the modelled depolarisation profiles for all clusters where X-ray data was available, for a single set of magnetic field parameters.

Fig. B.1. Depolarisation ratio and relevant upper limits for all sources as a function of projected radius and column density. Points are coloured by their position along the line-of-sight with respect to the nearest cluster.

Fig. B.2. Depolarisation ratio and relevant upper limits for all sources as a function of projected radius and column density. Points are coloured by the dynamic state of the cluster, indicated by non-cool-core (NCC) and cool-core (CC).

Fig. B.3. Depolarisation ratio and relevant upper limits for all sources as a function of projected radius and column density. Points are coloured according to the redshift of the cluster.

Fig. B.4. Full sample of measured values and upper limits on the depolarisation of radio sources as a function of both electron column density and projected radius. The points are coloured by the depolarisation value.

Fig. B.5. All simulated depolarisation ratio profiles for the full sample of simulated clusters. The assumed parameters for the magnetic field are B0 = 5.0 μG, n = 3 and η = 0.5.

Appendix C: Polarised source catalogue

Appendix D: Full source catalogue

Appendix E: Cluster catalogue

All Tables

All Figures

	Fig. 1. Central RMS noise in Stokes I in the 124 observed target fields.
In the text

Fig. 2. VLA 1–2 GHz primary beam corrected total intensity images of the 124 Planck clusters. The images are smoothed to the resolution of the lowest frequency channel (typically 6″) and the size is equal to the field-of-view at 2 GHz (0.35  ×  0.35 deg2). The colour scale is logarithmic with the scale range determined individually per cluster for visualisation purposes. The order of clusters follows the order in the Table E.1, in row-major order.

In the text

	Fig. 3. Same as Fig. 2, but for peak polarisation intensity. The colour scale used here is arcsinh. We note that the calibration artefacts visible in cluster position 107 (zero-based row-major index Beck & Krause 2005; Bonafede et al. 2011) were not used in the analysis, but the field was kept as two polarised sources at the edge of the primary beam were relatively unaffected by the artefacts.
In the text

Fig. 4. Results of the MCMC fitting of the external depolarisation model on the Stokes I, Q and U channels to a polarised source with good signal-to-noise ratio. Panel a shows the measured flux densities per channel in black and the best-fit model in blue. Panel b shows the posterior distribution of the model parameters visualised as one and two-dimensional projections in a corner plot. Panel c shows the inferred polarisation fraction from the data and the best-fit model together with the residuals.

In the text

	Fig. 5. Distribution of the sum of the squared residuals (χ2) of the model fits to all polarised components. The plot is truncated at χ2 = 400 for visibility.
In the text

	Fig. 6. Best-fit intrinsic polarisation fraction against the projected distance to the cluster centre, normalised by the cluster R500. Polarised sources are coloured based on their position along the line-of-sight with respect to the nearest cluster, defined in Sect. 5.2.
In the text

Fig. 7. Depolarisation against normalised projected radius. Left panel: full sample of sources and relevant upper limits. The error bars reflect the 68% confidence interval from the MCMC fitting procedure. Right panel: median of the Kaplan-Meier estimate of the depolarisation ratio survival function in different bins of projected radius to the cluster centre. The bin width is chosen such that each bin contains an equal number of polarised sources and is denoted by the horizontal lines (i.e. 0th and 100th percentile). The points are plotted at the median radius in each bin.

In the text

Fig. 8. Median of the Kaplan-Meier estimate of the depolarisation ratio survival function in different bins of radius and electron column density. Sources inside clusters are shown in green and sources behind clusters are shown in red. The bin width is chosen such that each bin contains an equal number of polarised sources and is denoted by the horizontal lines. The points are plotted at the median radius or column density in each bin.

In the text

	Fig. 9. Survival function (i.e. 1-CDF) inferred from the Kaplan-Meier estimator for all sources located at r/R500 > 0.5, and r/R500 < 0.5. Background sources and cluster members are indicated by red and green, respectively. The grey dashed line shows the location of the 50th percentile, indicating the median for both populations.
In the text

Fig. 10. Median of the Kaplan-Meier estimate of the depolarisation ratio survival function in different bins of radius. Cool-core clusters are shown in blue and non-cool-core clusters are shown in orange. The left panel shows bin widths (denoted by horizontal lines) chosen such that each bin contains an equal number of sources detected in polarisation and the right panel shows manually selected bins. The points are plotted at the median radius in each bin.

In the text

	Fig. 11. Survival functions (i.e. 1-CDF) inferred from the Kaplan-Meier estimator for all sources in the 0.0−0.2R500 and 0.2−1.0R500 bins, separately for sources detected around CC and NCC clusters. The grey dashed line shows the location of the 50th percentile, indicating the median for both populations.
In the text

	Fig. 12. Same as the left panel of Fig. 10, but separately for cluster members (left panel) and background sources (right panel).
In the text

	Fig. 13. Distribution of cluster mass in the low- and high-redshift samples.
In the text

	Fig. 14. Median of the Kaplan-Meier estimate of the depolarisation ratio survival function in bins of projected radii and column density, for the low- and high-redshift clusters. The bin width is chosen such that each bin contains an equal number of sources detected in polarisation and is denoted by the horizontal lines. The points are plotted at the median projected radius.
In the text

	Fig. 15. Distribution of the projected radius of sources detected in the low and high-redshift cluster samples. The solid blue histogram shows the result of one realisation of randomly sub-sampling sources in the low-redshift cluster sample to match the distribution of sources in the high-redshift sample.
In the text

	Fig. 16. Comparison of the median observed depolarisation ratio against projected distance between the merging clusters with and without detected radio halos.
In the text

	Fig. 17. Modelled depolarisation ratio as a function electron column density for six different clusters around z = 0.1. The assumed parameters for the magnetic field are B0 = 5.0 μG, n = 11/3 and η = 0.5, with Λmin = 4 kpc and Λmax = 1024 kpc.
In the text

	Fig. 18. Effect of the magnetic field parameters on the depolarisation profile of the cluster G039.85−39.98, with R500 = 1.2 Mpc. The varied parameters are the central magnetic field strength B0 in units of μG, the magnetic field power-law index n and the maximum correlation scale Λmax in kpc. For all but one model, the maximum correlation scale was set to the maximum allowed by the image.
In the text

	Fig. 19. Average depolarisation ratio profile over all simulated clusters using a Faraday screen located behind (blue dashed line) or inside (orange solid line) the cluster. The uncertainty interval indicates the standard error on the mean of the simulated profiles.
In the text

	Fig. 20. Comparison of the median observed depolarisation ratio in bins of projected radius, denoted by the horizontal lines, with the modelled depolarisation ratio profile for models with different magnetic field strengths (in μG) with η = 0.5. The model uncertainty interval indicates the standard error on the mean of the simulated profiles.
In the text

	Fig. 21. Comparison of the median observed depolarisation ratio with the modelled depolarisation ratio profile separately for the cool-core (CC) and non-cool-core (NCC) cluster sample. The modelled magnetic field parameters are B0 = 5 μG, η = 0.5 and n = 3 for both samples. The model uncertainty interval indicates the standard error on the mean of the simulated profiles.
In the text

	Fig. 22. Comparison of the median observed depolarisation ratio with the modelled depolarisation ratio profile separately for the low- and high-redshift cluster sample. The modelled magnetic field parameters are B0 = 5 μG, η = 0.5 and n = 3 for both samples. The model uncertainty interval indicates the standard error on the mean of the simulated profiles.
In the text

	Fig. 23. Median depolarisation ratio in different bins of radius without incorporating upper limits from unpolarised sources. Sources inside clusters are shown in green and sources behind clusters are shown in red. Compared to Fig. 8, there is significantly less depolarisation in the lowest radius bin.
In the text

	Fig. 24. Comparison of the median observed depolarisation ratio with the modelled depolarisation ratio profile using only clusters that have X-ray observations (blue) and all clusters (green) by calculating electron density profiles from a universal pressure profile (orange). The model uncertainty interval indicates the standard error on the mean of the simulated profiles.
In the text

	Fig. A.1. Various quantities against angular separation separately for cluster members and background sources. The running median is shown with uncertainties in the shaded region. Histograms show the projected distributions along the axes.
In the text

	Fig. A.2. Same as Figure A.1, but for total Stokes I flux density against depolarisation parameter σRM.
In the text

	Fig. A.3. Same as Figure A.1, but for intrinsic polarisation fraction against depolarisation parameter σRM.
In the text

	Fig. A.4. Probability distribution functions used to generate polarised radio sources during the Monte Carlo analysis shown in orange. The distribution of the real data is shown in blue.
In the text

	Fig. A.5. Median of the Kaplan-Meier estimate of the depolarisation ratio survival function in different bins of radius for simulated data. The bin width is chosen such that each bin contains an equal number of polarised sources and is denoted by the horizontal lines. The points are plotted at the median radius in each bin
In the text

	Fig. B.1. Depolarisation ratio and relevant upper limits for all sources as a function of projected radius and column density. Points are coloured by their position along the line-of-sight with respect to the nearest cluster.
In the text

	Fig. B.2. Depolarisation ratio and relevant upper limits for all sources as a function of projected radius and column density. Points are coloured by the dynamic state of the cluster, indicated by non-cool-core (NCC) and cool-core (CC).
In the text

	Fig. B.3. Depolarisation ratio and relevant upper limits for all sources as a function of projected radius and column density. Points are coloured according to the redshift of the cluster.
In the text

	Fig. B.4. Full sample of measured values and upper limits on the depolarisation of radio sources as a function of both electron column density and projected radius. The points are coloured by the depolarisation value.
In the text

	Fig. B.5. All simulated depolarisation ratio profiles for the full sample of simulated clusters. The assumed parameters for the magnetic field are B0 = 5.0 μG, n = 3 and η = 0.5.
In the text