The detection of cluster magnetic fields via radio source depolarisation

It has been well established that galaxy clusters have magnetic fields. The exact properties and origin of these magnetic fields are still uncertain even though these fields play a key role in many astrophysical processes. Various attempts have been made to derive the magnetic field strength and structure of nearby galaxy clusters using Faraday rotation of extended cluster radio sources. This approach needs to make various assumptions that could be circumvented when using background radio sources. However, because the number of polarised radio sources behind clusters is low, at the moment such a study can only be done statistically. In this paper, we investigate the depolarisation of radio sources inside and behind clusters in a sample of 124 massive clusters at $z<0.35$ observed with the Karl G. Jansky Very Large Array. We detect a clear depolarisation trend with the cluster impact parameter, with sources at smaller projected distances to the cluster centre showing more depolarisation. By combining the radio observations with ancillary X-ray data from Chandra, we compare the observed depolarisation with expectations from cluster magnetic field models using individual cluster density profiles. The best-fitting models have a central magnetic field strength of $5-10\,\mu$G with power-law indices between $n=1$ and $n=4$. We find no strong difference in the depolarisation trend between sources embedded in clusters and background sources located at similar projected radii, although the central region of clusters is still poorly probed by background sources. We also examine the depolarisation trend as a function of cluster properties such as the dynamical state, mass, and redshift. Our findings show that the statistical depolarisation of radio sources is a good probe of cluster magnetic field parameters. [abridged]


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 ha-Full Tables C.1, D.1 and E.1 are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/665/A71 los) place the magnetic field strengths of galaxy clusters around the µG level (Ferrari et al. 2008). Recently, observations of highredshift 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 , 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;Guidetti et al. 2008; Article number, page 1 of 32 arXiv:2207.09717v2 [astro-ph.CO] 20 Sep 2022 A&A proofs: manuscript no. main Laing et al. 2008;Bonafede et al. 2010;Guidetti et al. 2010;Vacca et al. 2012;Govoni et al. 2017). 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 doublelobed 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 et al. 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 Section 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 Section 3. Theoretical depolarisation expectations are derived through modelling of the magnetic fields as Gaussian random fields in Section 4 and results are shown in Sections 5 and 6. Finally, we conclude with a discussion and summary in Sections 7 and 8. Possible biases are discussed in Appendix A. Throughout this paper, we assume a flat ΛCDM model with H 0 = 70 kms −1 Mpc −1 , Ω m = 0.3 and Ω Λ = 0.7. We refer to the intensity of linearly polarised light simply as the polarised intensity.

Chandra-Planck ESZ sample
The Planck Early Sunyaev-Zel'dovich (ESZ) results presented 189 cluster candidates all-sky (Planck Collaboration et al. 2011), of which 163 clusters are at a redshift of z < 0.35. The Chandra-Planck Legacy Program for Massive Clusters of Galaxies 1 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 highquality X-ray data available. The Chandra observations of 147 clusters from the ESZ sample are presented in Andrade-Santos et al. (2017, 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).

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 et al. 2015) and PSZ2 (Planck Collaboration et al. 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 minutes 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 CASA TFCrop 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 offaxis leakage in more detail in Section 7.4. Finally, the antennabased complex gain solutions were calculated using the calibrator sources, and another round of automatic flagging was performed using the CASA TFcrop 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 ) 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-tonoise 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 -joinpolarizations 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 arcseconds. The distribution of central root-meansquare (RMS) noise in the full-band Stokes I images is given in Figure 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.
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, 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 Figure 2. The five fields with RMS noise higher than 60 µJy/beam in Figure 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.

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.

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 where λ indicates the observed wavelength, p 0 the polarisation fraction, I refers to the Stokes I intensity and 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 where n e is the electron density in parts per cm −3 , 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 Faraday rotation may cause polarised sources to be undetected in the wide-band Stokes Q and U images or in the linearly polarised intensity ( 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. RMsynthesis aims to approximate the Faraday dispersion function F(φ) by Fourier inversion of the following equation where P(λ 2 ) 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 λ 2 0 . 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: δλ 2 ||φ max || ≈ √ 3 δλ 2 ≈ 1200 rad m −2 .
The resolution in φ space is determined by our wavelength coverage, with the full-width half-maximum given by 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 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 2 rad m −2 , going up to 10 3 rad m −2 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 2 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 −2 in steps of 10 rad m −2 . 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 Figure 3.
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 arcminutes in steps of 1 arcminute 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.

Fractional polarisation measurement
To determine the polarisation properties such as the intrinsic polarisation fraction p 0 and the Faraday depth φ of polarised source candidates, we can model the polarised emission as 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) where Σ 2 RM 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  Figure 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 [10,7]) 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.
(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) where σ 2 RM 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 Equation 7.
Distinguishing between internal and external depolarisation effects can be done by measuring the spectral index of the po-larised 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 Sec. 5.) Therefore we decided to fit only the external depolarisation model given by Equation 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 Section 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 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 Section 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: The initial values for the parameters were found through a leastsquare 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 1 000 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 Figure 4.
To judge whether the model m, given by Equation 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, 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 Figure 5 we plot the distribution of the sum of the squared residuals (i.e. the χ 2 value) of the best fitting external depolarisation model to each polarised component. Most polarised sources have around 84 to 89 data points (i.e. channels) after masking the bad channels. This would give 80 to 85 degrees of freedom if the model was linear with 4 parameters. For comparison, we show also the theoretical χ 2 distribution with 80 degrees of freedom. The main peak of the sum of the squared residuals shows good agreement with the theoretical χ 2 distribution, indicating that most sources have acceptable fits. There is however a long tail of large χ 2 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 χ 2 value that is 5σ away from the We note that there are no apparent correlations between the χ 2 parameter and the derived best-fitting parameters p 0 , χ0, φ, σ 2 RM , 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 χ 2 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.
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 −1 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 −1 . The one-sigma upper limit on the fractional polarisation value is then calculated as 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).

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 5 806 out of 6 807 total intensity source candidates. For the remaining 1 001 sources, the best candidate optical coun-terpart 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 6 807 total intensity source candidates, 411 candidates were components of another source, leaving 6 396 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).

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 3 ) and the Sloan Digital Sky Survey (SDSS DR16; Ahumada et al. 2020) with a matching radius of 0.5 and 3 arcseconds 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 Section 2.3 of Tarrío & Zarattini (2020). To compute the photo- Notes. (a) The 'z best source' key is used in the catalogue presented in Table C.1 to indicate the origin of the redshift estimate.
metric 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% (4 544/6 807) 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.

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. 2011Bonafede et al. , 2013Vacca 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Ã. 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 where k denotes the magnitude of the three-dimensional wavevector k. The wave numbers k are related to the spatial scales Λ as k = 0.5 · 2π Λ , 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Ã 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 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 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), where B 0 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): Given the modelled magnetic field and observed electron density profile, we calculated the expected Faraday rotation in the clusters by numerical integration of Equation 4. Then, assuming an intrinsic polarisation p 0 , and a single polarisation angle χ 0 for a polarised background screen, we computed the polarisation angle of the radio emission at 1.5 GHz at the cluster redshift (χ obs ) with Equation 5. The predicted Stokes Q and U intensities were then obtained by inversion of P = Q 2 + U 2 and Equation 3: 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. 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 Figure 6 the intrinsic polarisation fraction p 0 (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.

Results -Observations
To minimise the effect of the scatter introduced by sourcedependent 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.5GHz to the intrinsic polarisation fraction p 0 , 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 nonparametric statistic used to estimate the complement of the cumulative distribution function, called the survival function. With x 1 < x 2 < .. < x r denoting distinct, ordered, observed values, the survival function is given by: 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 bŷ 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 leftcensored 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.5GHz as explained in Section 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 = 0.022, which is the median of the Kaplan-Meier estimator of all polarised radio sources detected at r > 1.5R 500 . All KM estimates were calculated using the lifelines 4 package (Davidson-Pilon 2019). 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 1 000 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. 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 500 , 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.

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 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 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 500 , 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 Figure 8 we show the depolarisation ratio calculated in different bins of normalised projected distance or electron col-A&A proofs: manuscript no. main umn 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 Figure B.1.
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 500 ), 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 Figure  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.
Lastly, at the inner radii (r < 0.5R 500 ) 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. P(dp > DP ) > 0.5R 500 -Behind cluster > 0.5R 500 -Inside cluster < 0.5R 500 -Behind cluster < 0.5R 500 -Inside cluster Fig. 9: Survival function (i.e. 1-CDF) inferred from the Kaplan-Meier estimator for all sources located at r/R 500 > 0.5, and r/R 500 < 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.

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.0 R 500 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 Figure B.2. We see a hint in the first radius bin of detecting more depolarisation in NCC clusters than in CC clusters.
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 Figure 10. We have chosen an inner radius bin of 0.0 − 0.2R 500 because the effect of the cooling core is significant only in the inner ∼ 0.2R 500 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 500 . In fact, sources detected at r < 0.2R 500 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 200 , 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 500 and 0.2 − 1.0R 500 bins. The survival functions are shown in Figure 11 and the log-rank test yields p-values of 0.22 and 0.001 for the 0.0 − 0.2R 500 and 0.2 − 1.0R 500 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.
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 500 , 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.

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 et al. 2016). This means that a cut in redshift effectively also corresponds to a mass-cut, as shown in Figure 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. Figure 14 shows the depolarisation trend for low-and highredshift clusters separately, and the full sample is shown in Figure 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 Figure 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-Article number, page 13 of 32 A&A proofs: manuscript no. main  and high-redshift sample of clusters. However, this is likely due to the low number of sources in the high-redshift sample.

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 Figure 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.

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 Section 4. We used the density profiles presented in Andrade-Santos et al.  Fig. 16: Comparison of the median observed depolarisation ratio against projected distance between the merging clusters with and without detected radio halos.

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 0 = 5.0 µG, n = 11/3 and η = 0.5 with a box-size of 1024 3 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 Figure 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.

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 0 . The effect of increasing the magnetic field is easily understood to result in more depolarisation because the scatter in RM increases Notes. (a) N denotes the amount of clusters with X-ray observations available such that models could be generated.
(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 Figure 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.
The degeneracy between a steep and strong magnetic field (e.g. B 0 = 10 µG, n = 4) and a shallower and weaker magnetic field (e.g. B 0 = 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.

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 3 where N is the number of pixels) of generating many magnetic field cubes out to typical values of R 500 , we decided to simulate different clusters with different resolution depending on the cluster redshift. To simulate the depolarisation effect, the G039.85-39.98 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 −1 , which made the cube size unfeasibly big to generate. All modelled depolarisation profiles are shown in Figure 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.

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 √ 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, Figure 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 Figure 8, where no clear difference between background sources and cluster members was found.

Average magnetic field properties
As shown in Section 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 0 = [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 500 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 Figure 20. Because of the degeneracy between n and B 0 (shown in Section 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 0 = 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.

Dynamical state
Section 5.3 showed that NCC clusters appear to cause significantly more depolarisation than CC clusters outside 0.2R 500 . 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 Figure 10.

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 500 , 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 lowredshift, low-mass sample of clusters.

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 Figure 7. In this section, we discuss the implications of the results and the possible limitations of this study.

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. 2011Guidetti et al. , 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 Figure 9 demonstrated that, although near the cluster centres (i.e. r < 0.5R 500 ) 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 500 ) the depolarisation of cluster members and background radio sources is similar.
In Section 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 (Figure 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.

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 Figure 20 showed that for all power spectrum indices, a central magnetic field strength higher than B 0 = 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;Guidetti et al. 2008;Laing et al. 2008;Bonafede et al. 2010;Guidetti et al. 2010;Vacca et al. 2012;Govoni et al. 2017;) and values from magnetohydrodynamic 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.

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 64 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 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) Clarke et al. 2001;. 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 Figure 10. Surprisingly, we found that NCC clusters show more depolarisation than CC clusters outside the cooling-core region defined as r > 0.2R 500 (Figure 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 Figure 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 lowmass 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 Section 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).

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 λ 2 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 Section 3.2 that about 25% (193/819) polarised sources detected in this work are not wellfitted by the single component external depolarisation model. While fitting these sources with more complicated models (e.g. Article number, page 19 of 32 A&A proofs: manuscript no. main Brown et al. 2019) is beyond the scope of this work, we can briefly investigate how many sources show evidence of repolarisation 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 0 had to be made. For these sources we assumed p 0 = 0.022, which was the median intrinsic polarisation fraction of sources detected at r > 1.5R 500 . 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 6 807 unpolarised sources have an upper limit on the 1.5 GHz polarisation fraction that is below the assumed p 0 = 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 Figure 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.
Thirdly, for the upper limits, we computed a polarised flux threshold as described in Section 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 Figure 6 and Figure 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 500 , R 500 and P 500 , where P 500 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 et al. (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 Figure 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. Combined (109) Data 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.

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 500 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.2R500 in coolcore 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.
Acknowledgements. 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

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.

Appendix 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), 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.
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 0 and σ RM (see Fig. 4). While this degeneracy is simply a result of the fitting, real correlations between p 0 and σ RM have been claimed before in the literature (e.g. Lamee et al. 2016). It is therefore important to inspect trends of p 0 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 0 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 0 and σ RM , particularly for large values of σ RM . It is interesting that this trend implies that more depolarised sources have larger p 0 , while Figure A.1c showed that smaller values of p 0 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.

Appendix A.2: A2. 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 ex-

Notes.
The columns 'Maj' and 'Min' indicate the major and minor axis of the polarised sources and 'PA' the position angle. The upper limit on the polarisation fraction p 1.5GHz was determined for unpolarised radio sources as described in section 3.2. The columns starting with z are the results of the redshift estimates as detailed in Section 3.4. The angular distance to the pointing centre and projected radius to the cluster centre are given by the θ p and r/R 500 columns respectively. The 'Visual counterpart id' column indicates whether the source has undergone visual inspection to mark the optical counterpart or whether it was determined automatically as described in Section 3.3.
Article number, page 31 of 32  Notes.
The cluster centre coordinates from the X-ray profiles are given by cRA and cDEC and pointing coordinates are given by pRA and pDEC.
The dynamical state was determined as described in Section 5.3 and the presence of a radio halo by a literature search.