© The Authors 2023

1. Introduction

It is well established that galaxy clusters grow over cosmic time via accretion of infalling matter (e.g., sub-clusters or galaxy groups) along the filaments that constitute the cosmic web (e.g., Press & Schechter 1974). Mergers of sub-clusters and galaxy groups build the most massive virialised structures in the Universe, but they also generate shock waves and trigger magneto-hydrodynamical turbulence in the intracluster medium (ICM; see Markevitch & Vikhlinin 2007). The presence of diffuse radio emission in clusters, which is not associated with black holes in galaxies (i.e. active galactic nuclei or radio galaxies), is usually associated with these merger events. This radio emission is linked to the presence of (re-)accelerated particles, that is, cosmic rays (CRs), with a Lorentz factor γ_L ≳ 10³, and magnetic fields, on average with levels of a few μ Gauss (Brunetti & Jones 2014, for a theoretical review). Diffuse radio emission in merging galaxy clusters can generally be divided into two categories¹: radio relics and radio halos (van Weeren et al. 2019, for an observational review).

Radio relics are elongated structures with sizes that can reach up to ∼2 Mpc. Recent high-resolution GHz-frequency observations have shown that they exhibit filamentary morphologies (Owen et al. 2014; Di Gennaro et al. 2018; Rajpurohit et al. 2020, 2021c; de Gasperin et al. 2022). They are usually detected in the cluster outskirts (e.g., Roettiger et al. 1999; van Weeren et al. 2010; Pearce et al. 2017; Di Gennaro et al. 2018; Vazza et al. 2018) in the proximity of shock discontinuities detected in X-rays (e.g., Finoguenov et al. 2010; Akamatsu et al. 2015; Urdampilleta et al. 2018; Di Gennaro et al. 2019). Therefore, these sources are associated with merger shocks. The most favourable formation scenario involves particle acceleration and magnetic field compression and amplification due to merger-induced shock propagation. At the shock location, first-order Fermi acceleration mechanisms provide energy to the thermal electrons in the ICM, which then become ultra-relativistic (i.e. CRs). This is also known as the diffusive shock acceleration mechanism (Blandford & Eichler 1987). These CRs then lose energy via synchrotron and inverse Compton (IC) radiation and fade away in the post-shock region. In support of this scenario, the spectral index (α, being S_ν ∝ ν^α) of these sources is flat (i.e. α ∼ −0.8) at the shock location and steepens towards the cluster centre (i.e. α ≲ −1.2). However, several issues arise in the acceleration of thermal electrons from the ICM. For instance, cluster-merger shocks show in X-ray observations too low Mach numbers (ℳ ∼ 2 − 3) to accelerate particles efficiently from the thermal pool (e.g., Botteon et al. 2020). Moreover, a recent study on a sample of relics at ∼150 MHz with the Low Frequency Array (LOFAR) strongly suggests that the radio emission in the post-shock region extends too far downstream (Di Gennaro et al. 2018; Rajpurohit et al. 2018) if we simply consider synchrotron and IC energy losses (Jones et al. 2023). Observational evidence (e.g., Markevitch et al. 2005; van Weeren et al. 2017) and numerical simulations (Kang et al. 2017) have suggested the presence of pre-existing electrons that are re-accelerated by the shock. This idea is supported by observations of radio galaxies that appear to be close to radio relics. These could provide the source of fossil plasma needed in the re-acceleration scenario (van Weeren et al. 2017; Di Gennaro et al. 2018). Recently, another class of shocks have been proposed, the so-called equatorial shocks. These kinds of events are more difficult to observe, as they are thought to be formed in the very first stages of the cluster merger (Ha et al. 2018). The only equatorial shock known to date is the one in the cluster pair 1E 2216.0-0401 and 1E 2215.7-0404 at redshift z = 0.09 (Gu et al. 2019).

Radio relics are also found to be highly polarised (up to 60% at GHz frequencies, Ensslin et al. 1998). This is consistent with the picture of shock acceleration formation, as a shock wave can compress and amplify cluster magnetic fields (Iapichino & Brüggen 2012; Donnert et al. 2018; Wittor et al. 2019; Domínguez-Fernández et al. 2021; Hoeft et al. 2022). Recent GHz-frequency high-resolution observations have shown that the filamentary morphology is also detected in polarised emission (Di Gennaro et al. 2021c; Rajpurohit et al. 2022a; de Gasperin et al. 2022).

Radio halos have roundish shapes with Mpc sizes and emission that follows the ICM distribution. Several recent studies have shown a correlation between the cluster X-ray luminosity and mass and the halo radio emission (e.g., Cassano et al. 2013; Cuciti et al. 2021), as well as a correlation between the presence of radio halos and the cluster disturbance (e.g., Cassano et al. 2010; Cuciti et al. 2021). The most favourable scenario involves second-order Fermi acceleration due to merger-induced magneto-hydrodynamical turbulence. Turbulence stochastically re-accelerates electrons and triggers a small-scale dynamo that causes magnetic field amplification (e.g., Brunetti et al. 2001; Petrosian 2001; Brunetti & Lazarian 2007, 2016; Brunetti & Blasi 2005; Donnert et al. 2013). Additionally, re-acceleration of secondary electrons emerging from proton-proton collisions might also support the synchrotron radiation from radio halos (Brunetti & Lazarian 2011; Pinzke et al. 2017; Brunetti et al. 2017) while still being consistent with γ-ray limits from Fermi-LAT observations (e.g., Adam et al. 2021). Radio halos are mostly characterised by a rather uniform spectral index distribution (α ∼ −1.3), although ultra-steep spectra radio halos have also been found (i.e. α ≲ −1.5, e.g., Brunetti et al. 2008; Dallacasa et al. 2009; Bonafede et al. 2012; Wilber et al. 2018; Bruno et al. 2021; Duchesne et al. 2021; Di Gennaro et al. 2021b). The existence of these ultra-steep spectra sources is one of the most stringent predictions of turbulent re-acceleration models (Cassano et al. 2006, 2012, 2023).

Unlike radio relics, radio halos are yet to be detected in polarisation. Magnetic field profiles and average magnetic field values have been obtained for a handful of systems (e.g., Bonafede et al. 2013; Stuardi et al. 2019; Rajpurohit et al. 2022b; Osinga et al. 2022) via Rotation Measure (RM) or depolarisation analysis of polarised background and/or embedded radio galaxies. These studies suggest average magnetic fields of a few μ Gauss with higher values in the cluster centre. Recently, μ Gauss-level magnetic fields in the cluster volume were also estimated in a sample of high-redshift (i.e. z = 0.6 − 0.9) clusters observed with the LOFAR at 144 MHz (Di Gennaro et al. 2021a). Although these observations do not provide any constraints on the magnetic seeds, they yield strong limits on the evolution of magnetic fields when the Universe was only half of its age and only a few Gyr after the first large-scale structure formed. Moreover, these observations provide additional evidence on the re-acceleration mechanisms (Cassano et al. 2019), as all the less massive clusters (i.e. M_500, SZ < 6 × 10¹⁴ M_⊙) exhibit steep spectral index halos (i.e. α < −1.5, Di Gennaro et al. 2021b).

2. PSZ2 G091.83+26.11

The cluster PSZ2 G091.83+26.11 (hereafter, PSZ2G091; see Fig. 1 and Table 1) was spectroscopically confirmed to be at z = 0.822 by Amodeo et al. (2018). With a mass of M_SZ, 500 = (7.4 ± 0.4)×10¹⁴ M_⊙ estimated from the total Sunyaev-Zel’dovich (SZ) intensity (i.e. Y_SZ, 500 = (2.3 ± 0.3)×10⁻⁴ Mpc²; see Planck Collaboration XXVII 2016), it is the third most massive cluster hosting diffuse radio emission at z > 0.8 (“el Gordo” at z = 0.870 with M_SZ, 500 = (1.17 ± 0.17)×10¹⁵ M_⊙, Lindner et al. 2014, and PSZ2 G160.83+81.66 at z = 0.888 with $M_{\mathrm{SZ},500}=5.7_{-0.7}^{+0.6}\times {10}^{14}{M}_{\odot }$, Di Gennaro et al. 2021a,b, are the other two). In particular, observations with LOFAR at 144 MHz and upgraded Giant Metrewave Radio Telescope (uGMRT) at 400 MHz and 650 MHz showed that the cluster hosts a giant radio halo and an additional elongated source in the east, south-east side of the cluster (Di Gennaro et al. 2021a,b). Due to its elongated shape and hints of spectral steepening towards the cluster centre, the latter was associated with a candidate radio relic. Preliminary inspection of the X-ray morphology with Chandra and XMM-Newton showed that the ICM emission of the cluster is highly disturbed, with two main peaks of emission, and is elongated in the north-east, south-west direction. More recently, Artis et al. (2022) found that the thermal SZ emission of the cluster with the New IRAM Kids Arrays (NIKA2) at 2 mm has a morphology very similar to the X-ray emission. By comparing the X-ray and the SZ surface brightness maps, they suggested that the two sub-clusters, identified by the two peaks of emission, could be in the first stages of a major merger.

Fig. 1. Composite image of PSZ2G091. Optical (white): PanSTARSS gri observations. Radio (red): 3.0 GHz VLA observations. X-ray (blue): Chandra 0.8–4.0 keV observations.

Interestingly, while the central radio emission follows the ICM emission well, resulting in the radio halo classification, the position of the candidate radio relic is parallel to the ICM elongation and not perpendicular, as is usually observed. In this paper, we therefore focus on the total intensity and polarisation emission of the candidate relic in PSZ2G091 in order to shed light on the origin of this radio source and on the cluster-merger state. To this end, we used Karl Jansky Very Large Array (VLA) observations in the L and S band (i.e. covering the 1–4 GHz frequency range). We complemented the radio data with Chandra X-ray observations. The paper is organised as follows: In Sect. 3, we describe our observations. In Sects. 4 and 5, we describe and discuss the results of our analysis. Finally, we conclude the paper with a summary in Sect. 6.

Throughout the paper, we assume a standard ΛCDM cosmology, with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7. This translates to a luminosity distance of D_L = 5187.7 Mpc and a scale of 7.576 kpc/″ at the cluster redshift, z = 0.822.

3. Observations and data reduction

In this section, we describe the VLA and Chandra observation and calibration. We referred to Di Gennaro et al. (2021a,b) for the LOFAR (120–168 MHz) and uGMRT (250–500 MHz, band 3; 550–900 MHz, band 4) data, for which we repeated the imaging to match them with the present radio images.

3.1. VLA data

The cluster PSZ2G091 was observed with the Karl Jansky VLA in the L (1–2 GHz) and S bands (2–4 GHz) in the D, C, and B configurations (Table 2). The observations were mostly carried out during 2019 and 2021 (project codes: 19B-080 and 21A-025, PI: Di Gennaro), for a total observing time of 13 and 15 h in the L and S bands, respectively. The L band B configuration observations were complemented with archival snapshot observations (30 min on target, project code: 15A-270, PI: van Weeren).

The data reduction was carried out with CASA v5.4 (McMullin et al. 2007) using 3C286 as the bandpass, flux, and polarisation angle calibrator and J1927+6117 as the phase calibrator. In the absence of a standard leakage calibrator (e.g., 3C147), we used J1407+2827 instead². For both project codes, we followed the same data reduction strategy as described in Di Gennaro et al. (2018). Briefly, we split the wide-band observations into single spectral window (spw) datasets, and after a first removal of radio frequency interference (RFI) using the tfcrop mode, we calibrated the antenna delays, bandpass, cross-hand delays, and polarisation leakage and angle. Finally, we merged the spw datasets together and applied all the solutions to the target, averaging in time and frequency (factor of two and four, respectively). A final round of AOFlagger (Offringa et al. 2010) on the cross-hand polarisation (LR and RL) was run to remove additional RFI. Bad spw data were discarded from the final imaging. Due to bad quality, one L-band D configuration observation was also discarded. We performed self-calibration on the single configurations to refine the amplitude and phase solutions using WSClean v2.10 (Offringa et al. 2014; Offringa & Smirnov 2017) to produce the model images. Finally, we combined the visibilities for all the configurations together in order to produce the final images of the cluster at different resolutions (see Table 3). For the L-band dataset, an additional self-calibration on a bright source located at the edge of the primary beam (RA = 18^h31^m24.6^s Dec = +62°30′34.32″) was needed to improve the quality of the image (i.e.peeling; see Appendix A). Additionally, we performed a bandpass calibration on the source using the model derived from the self-calibration.

In order to perform a polarisation analysis, we created images with a spectral resolution of Δν = 8 MHz for all the Stokes parameters (i.e. I, Q, and U). Following Di Gennaro et al. (2021c), Stokes Q and U images for each channel were produced with WSClean with the options -join-channels and -join-polarizations. We note that the option -squared-channel-joining, which is used to prevent the polarisation from averaging to zero, was not used because Stokes Q and U were imaged separately³ using the Stokes I image to search for clean components (i.e. -pol IQ and -pol IU). We removed all the noisy or low-quality images as well as those covering flagged spw and re-gridded them to the same pixel grid and convolved to the same angular resolutions. All the VLA images presented in the manuscript were made with WSClean v2.10 with Briggs weighting and robust=0, unless otherwise stated, and using multiscale deconvolution.

3.2. Chandra

We complemented our VLA radio data with Chandra observations. The PSZ2G091 cluster was observed on 19 June 2016 with the Advanced CCD Imaging Spectrometer (ACIS) using the ACIS-I CCD configuration for a total time of 23 ks (see Table 4; PI: Rossetti). We processed the data with CIAO v4.14 (Fruscione et al. 2006), which applies the CALDB 4.9.8 calibration files. We used chandra_repro to generate the level=2 event file, including the option check_vf_pha=yes that screens the particle background in very faint mode observations for ‘bad’ events, which are most likely associated with cosmic rays. We screened the light curve of the observation using the CIAO tasks deflare and lc_clean. No time bins were identified where the count rate deviated from the mean by more than 2σ; the observation was therefore not affected by flares, and the full 23 ks exposure could be used for further science analysis. Standard blank sky files scaled to match the hard-band count rate of the observation were used for background subtraction. The final exposure-corrected image was made in the 0.8–4.0 keV energy band, where the signal-to-detector background ratio is optimal. Where appropriate, compact sources unrelated to the diffuse ICM were removed from the spatial and spectral analysis.

4. Results

4.1. Total intensity and polarised images

The VLA 1.5 GHz and 3.0 GHz images of PSZ2G091 (Fig. 2) are remarkably similar to those at lower frequencies (i.e. at 144 MHz, 400 MHz and 650 MHz, Di Gennaro et al. 2021b). At both frequencies, the radio halo is also visible in the full-resolution images (i.e. Θ_1.5 = 4.3″ × 3.2″ and Θ_3.0 = 2.9″ × 2.2″). The radio halo is elongated in the north-east south-west direction, and it follows the emission from the X-ray Chandra image, nearly completely filling the 0.5R_SZ, 500 area (see Sect. 4.2). The brightest radio feature, however, is the candidate radio relic located about 500 kpc to the east of the cluster. It is a narrow filamentary structure that is sharply broken into two separate pieces, R1 and R2, forming an angle of about 130°. Among the two pieces of emission, R1 is much fainter than R2, and it is not completely detected in the full-resolution images at both frequencies. The length of the two sub-filaments is about 80″ for R1 and 40″ for R2 at low resolution, corresponding to about 640 kpc and 300 kpc at the cluster redshift, respectively. This remains unchanged in both the 1.5 GHz and 3.0 GHz images as well as in the LOFAR and uGMRT images (Di Gennaro et al. 2021b). The radio source remains unresolved in its width, even at the highest resolution, that is, 1.8″ × 1.1″ (3.0 GHz image, with robust=-1.25), meaning that the filamentary structure is about 15 kpc (deconvolved, i.e. full width half maximum). At this resolution, we noticed that the south-west side of R2 (i.e. the brightest part of the candidate relic) appears to be broken (Fig. 3). No particular feature in the X-ray image is visible at the same location, although we stress that the X-ray observation is quite shallow (i.e. 23 ks) and the radio source is located where the X-ray emission is faint due to the low density of the ICM (see last panel in Fig. 3). Moreover, no optical galaxy appears to be associated with the radio source (see fourth panel in Fig. 3).

Fig. 2. VLA images of PSZ2G091 at 1.5 GHz (top row) and 3.0 GHz (bottom row). From left to right: Full-resolution (i.e. no taper applied) taper of 25 kpc, 50 kpc, and 100 kpc (see Table 3 for the final resolutions and map noise). The beam size is shown in the bottom-left corner of each panel. Labels in the leftmost panel follow those shown in Di Gennaro et al. (2021b). The white dashed circle shows the R = 0.5RSZ, 500 region, and the cluster centre is marked with a white cross. Radio contours, when shown, are drawn with white solid lines at the 2.5σrms × [1, 2, 4, 8, 16, 32] levels; a −2.5σrms level is also drawn with a dashed line.

Fig. 3. Zoom on the radio relic. From left to right: highest-resolution total intensity VLA image (i.e. 3.0 GHz at 1.8″ × 1.1″; see Table 3); total averaged polarisation intensity (without correction for Ricean bias, see Di Gennaro et al. 2021c) in the 2–4 GHz band (effective frequency 3.1 GHz) at 5″ resolution; high resolution-total intensity image at 3.0 GHz with the intrinsic polarisation magnetic field vectors at 5″ resolution, corrected for Faraday rotation; optical PanSTARRS gri image with radio contours from the leftmost panel starting from 2σrms; Chandra 0.8–4.0 keV image with radio contours from the leftmost panel starting from 2σrms.

The location of this radio feature is quite interesting. Commonly, radio relics are located perpendicularly to the merger axis, as shock waves are generated after the collision of two or more clusters. The X-ray extension of the ICM emission of PSZ2G091 suggests that the merger event happened – or is happening – in the north-east, south-west direction, and the candidate radio relic (at least R2) appears to be parallel to it.

In order to investigate the degree of polarisation of the cluster, we ran the pyrmsynth tool⁴, which combines the single-channel Stokes Q and U images to generate Faraday cubes through the RM synthesis technique (Brentjens & de Bruyn 2005). We assumed a Faraday depth in the range of ± 4000 rad m⁻², with a sampling bin of 2 rad m⁻². We produced Faraday spectra at both high and low resolutions (i.e. 5″ and 12.5″), the former only using observations in the 2–4 GHz band. These frequency bandwidths and resolutions correspond to a maximum Faraday depth ($||{\varphi }_{\max }||=\sqrt{(3)}/\delta {\lambda }^2$, being δλ² the channel width in wavelength squared) of ∼1400 rad m⁻², a resolution in Faraday space ($\delta \varphi =2\sqrt{(3)}/\mathrm{\Delta }{\lambda }^2$, being $\mathrm{\Delta }{\lambda }^2={\lambda }_{\text{max}}^2-{\lambda }_{\text{min}}^2$ the total bandwidth in wavelength squared) of ∼78 rad m⁻² and a maximum recovered scale in Faraday space ($\mathrm{\Delta }{\varphi }_{\max }=\pi /{\lambda }_{\min }^2$) of ∼558 rad m⁻² (Brentjens & de Bruyn 2005, see Table 5). The only part of the cluster where we detected linear polarisation is R2, with the polarisation magnetic field vectors well aligned with the direction of the radio surface brightness edge (see second and third panels in Fig. 3; see also Sect. 5.3). However, we note that the southern part of R2 is fainter in polarisation than the northern part. This part of the candidate relic is coincident with the ICM emission (see last panel in Fig. 3), which may suggest a larger depolarisation effect on this piece of diffuse radio emission. No polarised emission was detected from compact sources in the cluster area.

4.2. X-ray morphology

As also shown in Di Gennaro et al. (2021a, see Fig. 2 in that manuscript) and in Artis et al. (2022), the ICM morphology of PSZ2G091 is strongly disturbed, with the X-ray emission elongating in the north-east, south-west direction. We identified two main peaks of emission, one in the north and the other in the south, at coordinates RA_S = 18^h31^m09.2^s Dec_S = +62°14′09.5″ and RA_N = 18^h31^m12.8^s Dec_N = +62°14′58.7″ (Fig. 4). The distance of the two peaks from each other is ∼420 kpc, if we assume the cluster is perfectly on the plane of the sky.

Fig. 4. Thermal and non-thermal comparison. X-ray 0.8–4.0 keV Chandra image with 3.0 GHz radio contours tapered at 50 kpc (contour levels are the same as for Fig. 2). The radio beam size is shown in the bottom-left corner. The white dashed circle shows the R = 0.5RSZ, 500 region, and the white cross denotes the cluster centre. The two yellow stars show the cluster X-ray peaks, which are associated with the two candidate sub-clusters (Artis et al. 2022).

In order to determine the global X-ray cluster properties (i.e. temperature and luminosity), we extracted the spectrum within a circle of R = 1.01 Mpc (i.e. ∼2.2′), which corresponds to R_SZ, 500, given the cluster mass M_SZ, 500 centred on the cluster coordinates (see Table 1). We modelled the spectrum with a single-temperature (phabs*APEC) model that included the absorption from the hydrogen column density (N_H) of our Galaxy. We used the mean total value⁵, which includes both the molecular (H₂) and atomic (HI) hydrogen (Willingale et al. 2013), that is, N_H = 4.86 × 10²⁰ atoms cm⁻².

We found a global cluster temperature of $k{T}_{500}=12.7_{-1.5}^{+2.1}$ keV and an unabsorbed luminosity in the 0.1–2.4 keV range of L_{[0.1 − 2.4 keV],500} = (3.21 ± 0.14)×10⁴⁵ erg s⁻¹. Emission from highly ionised Fe was marginally detected even in this shallow data, constraining the average gas metallicity at 0.77±0.29 Solar using the reference units of Asplund et al. (2009).

4.3. Flux densities and integrated spectral indices

We obtained the relic flux density in both total intensity and total average polarisation from the 50 kpc tapered image (see Fig. 2). In this way, we were able to retrieve all the diffuse radio emission while minimising the blending with the radio halo. For the total intensity, in order to have a wide frequency coverage, we also used the 144 MHz LOFAR and the 400 MHz and 650 MHz uGMRT observations tapered at the same value (see Appendix B). All these images were re-gridded at the same pixel size and convolved at the same common resolution (i.e. 12″, see Table 6)⁶. The flux densities were obtained from the same region at all frequencies, which also includes source C (see Fig. 2 and right panel in Fig. 5). The flux density of this source was then mathematically subtracted⁷. The uncertainty on the radio flux densities for the radio relic were calculated as

$$\begin{array}{c}\hfill \mathrm{\Delta }{S}_{\nu ,\mathrm{relic}}=\sqrt{\mathrm{\Delta }{S}_{\nu ,\mathrm{total}}+\mathrm{\Delta }{S}_{\nu ,\mathrm{source}C}},\end{array}$$(1)

Fig. 5. Integrated spectra of the diffuse radio emission in PSZ2G091. Left: the flux densities were measured from the solid-line regions in the right panel, and the spectral indices are reported in the bottom-left corner in the panel (black for the halo and shades of purple for the relic). Right: Radio image and contours are from the 3.0 GHz VLA image at 12″ resolution. The separation of R1 and R2 is indicated by the purple solid line in the right panel. The dashed dashed-line sectors refer to the profiles displayed in Fig. 8. White circles mark the position of the compact sources that have been masked for the flux density, spectral and curvature profiles.

where

$$\begin{array}{c}\hfill \mathrm{\Delta }{S}_{\nu }=\sqrt{{(f{S}_{\nu })}^2+N_{\mathrm{beam}}{\sigma }_{\mathrm{rms}}^2}.\end{array}$$(2)

Here, f is the systematic uncertainty due to residual amplitude errors (we assumed 15% for LOFAR, 5% for uGMRT, and 3% for VLA, see Shimwell et al. 2022; Chandra et al. 2004; Perley & Butler 2013, respectively), σ_rms is the map noise level, and N_beam is the number of beams covering the relic region and source C. We measured 248.8 ± 40.9 mJy, 66.5 ± 3.9 mJy, 34.6 ± 2.0 mJy, 15.9 ± 0.6 mJy, and 5.0 ± 0.2 mJy for the full relic (i.e. R1+R2) at 144 MHz, 400 MHz, 650 MHz, 1.5 GHz, and 3.0 GHz, respectively (see Table 6). We also calculated the flux density for R1 and R2 separately at the same frequencies and resolution (Table 6). As is clear from the images, R2 represents the main source of emission, while R1 is barely detected at the highest frequency. We fit the flux densities with a power-law of the type y = ax + b, where y = log S, x = log ν, and the slope $a={\alpha }_{144\mathrm{MHz}}^{3.0\mathrm{GHz}}$ is the spectral index between 144 MHz and 3.0 GHz. We obtained ${\alpha }_{144\mathrm{MHz}}^{3.0\mathrm{GHz}}=-1.24\pm 0.03$, ${\alpha }_{144\mathrm{MHz}}^{3.0\mathrm{GHz}}=-1.16\pm 0.03$, and ${\alpha }_{144\mathrm{MHz}}^{3.0\mathrm{GHz}}=-1.25\pm 0.03$, for the full relic, R1, and R2, respectively (see Fig. 5). Given these spectral indices and flux densities, we calculated the radio luminosity at each frequency, L_ν, according to:

$$\begin{array}{c}\hfill L_{\nu }=\frac{4\pi D_L^2}{{(1+z)}^{1+\alpha }}{S}_{\nu },\end{array}$$(3)

where D_L is the luminosity distance at the cluster redshift (i.e. 5188 Mpc). These radio luminosities are reported in the last column of Table 6. Uncertainties regarding radio luminosity take into account both the flux density and the spectral index uncertainties through Monte Carlo simulations.

In order to estimate reliable values for the averaged polarised flux, we needed to correct for the Ricean bias (Rice 1945). Indeed, the polarised flux is given by:

$$\begin{array}{c}\hfill P=\sqrt{Q^2+U^2}.\end{array}$$(4)

This meant that we would always measure a positive polarised emission even where no “real” polarised signal is present (i.e. where the Stokes Q and U noise dominates). Following Stuardi et al. (2021), we corrected for the Ricean bias as $P=\sqrt{|F({\varphi }_{\mathrm{peak}}){|}^2-2.3{\sigma }_{\mathit{QU}}^2}$ (George et al. 2012). Here, |F(ϕ_peak)| is the peak of the Faraday dispersion function F(ϕ) obtained from the RM synthesis technique, and σ_QU is the average noise in the Stokes Q and U maps (σ_QU = 3.3 μJy beam⁻¹ at 12.5″). After this correction, we obtained P_3.0 GHz = 0.5 mJy⁸ for R2 at 12.5 ″ resolution, while we set an upper limit for R1 of P_3.0 GHz < 0.06 mJy.

We also measured the flux density from the radio halo. In order to do so, we produced source-subtracted images, that is, we cut all the visibilities corresponding to scales larger than 500 kpc to create a model of “compact” sources that were then subtracted from the visibilities. We then tapered the uv data at lower resolutions. This approach was also performed on LOFAR and uGMRT data (Di Gennaro et al. 2021b), and it also included multiscale deconvolution to take into account the diffuse emission associated with radio galaxies. However, similarly to the LOFAR and uGMRT cases, this procedure leaves part of the emission of the candidate radio relic, which subsequently needs to be masked manually. We chose the 100 kpc tapered images to retrieve all the diffuse emission from the radio halo (see rightmost panel in Fig. 2) at the VLA, uGMRT, and LOFAR frequencies. Similarly to the radio relic, to retrieve the integrated spectral index, we also added the information from the 144 MHz LOFAR and 400 MHz and 650 MHz uGMRT observations. In this case, the uncertainties on the flux densities are given by $\mathrm{\Delta }{S}_{\nu }=\sqrt{{(f{S}_{\nu })}^2+N_{\mathrm{beam}}{\sigma }_{\mathrm{rms}}^2+{\sigma }_{\mathrm{sub}}^2}$, where σ_sub are the residuals on the subtraction (i.e. a few percent of the residual flux from compact sources; see Cassano et al. 2013). We measured 73.8 ± 11.1 mJy, 28.3 ± 1.7 mJy, 17.9 ± 1.0 mJy, 6.1 ± 0.2 mJy, and 3.1 ± 0.1 mJy at 144 MHz, 400 MHz, 650 MHz, 1.5 GHz, and 3.0 GHz, respectively. We note that for the LOFAR and uGMRT flux densities, we obtained values that are slightly higher than those reported in Di Gennaro et al. (2021b), although they are within the error bars. Assuming a power-law, the total integrated spectrum is ${\alpha }_{144\mathrm{MHz}}^{3.0\mathrm{GHz}}=-1.06\pm 0.03$ (in agreement with that reported in Di Gennaro et al. 2021b). The radio luminosities at each frequency are reported in Table 6.

4.4. Spectral index and curvature maps

We produced spectral index (α) and curvature (C) maps for the full cluster at different resolutions using the approach already described in Di Gennaro et al. (2021b). The images were re-gridded at the same pixel scale and convolved and at the same resolution. Moreover, we checked and fixed for possible offset in the astrometry at the different frequencies, selecting point-like sources and cross-matching their flux peak. Using the VLA 3.0 GHz image as the reference, as it has the most accurate astrometry⁹, we found a separation of Δ_RA = 1″ and Δ_Dec = 1.5″ for the LOFAR and the uGMRT observations and a separation of about Δ_RA = Δ_Dec = −0.5″ for the VLA 1.5 GHz.

To produce the low-resolution (12″) spectral index and curvature maps, we used all the frequencies available (i.e. the LOFAR, uGMRT, and VLA) with a 50 kpc taper. For each pixel, we fit a second-order polynomial (i.e. y = ax² + bx + c) if the curvature parameter was larger than 2σ, where σ is the uncertainty associated with the second-order term. In this case, the spectral index and the curvature were calculated at the median of the total band (i.e. ν_ref = 650 MHz, as 2a log ν_ref + b and C_650MHz = a, respectively); otherwise, we fit a first-order polynomial and the spectral index was simply given by the slope of the fit, while the curvature parameter was set to zero. According to this convention, in the second-order fit, negative values of curvature correspond to steeper spectral indices at high frequencies; in other words, the spectrum is convex (see also, e.g., Di Gennaro et al. 2018; Stuardi et al. 2019; Rajpurohit et al. 2020). We blanked all the pixels below the 2.5σ_rms threshold for each frequency. The spectral index uncertainty maps were obtained via 150 Monte Carlo simulations of the first- or second-order polynomial fit. We assumed that the uncertainty of each flux was given by the sum in quadrature of the noise map and the systematic flux uncertainties, that is, $\mathrm{\Delta }{S}_{\nu }=\sqrt{{(f{S}_{\nu })}^2+{\sigma }_{\mathrm{rms}}^2}$. We found hints of steepening from the outer edge of the relic towards the cluster centre (i.e. from α_650 MHz ∼ −1 to α_650 MHz ∼ −1.5). For R2, we also found a significant curvature C_650 MHz ∼ −0.4 (see Fig. 6), meaning that the spectral index steepens at higher frequencies. The flattening near source C is probably due to the blending of the radio galaxy spectral index and the outermost edge of the relic, as this is not visible in the spectral index maps at higher resolutions (i.e. at 4.5″ using the 144 MHz, 650 MHz, 1.5 GHz, and 3.0 GHz maps, see left panel in Fig. 6, and at 3″ using the 1.5 GHz and 3.0 GHz maps, see Fig. 7). The highest-resolution spectral index map (i.e. 3″) also showed steepening towards the cluster centre for the lower resolutions (Fig. 7).

Fig. 6. Spectral index and curvature maps (top) with corresponding uncertainty maps (bottom) of PSZ2G091. Left column: resolution of 4.5″ using the 144 MHz, 650 MHz, 1.5 GHz, and 3.0 GHz observations. Central and Right columns: Resolution of 12″ using the 144 MHz, 400 MHz, 650 MHz, 1.5 GHz, and 3.0 GHz observations. Radio contours are from the 3.0 GHz observations, drawn at the 2.5σrms × [1, 2, 4, 8, 16, 32] levels (with ${\sigma }_{\mathrm{rms},4.5″}=3.1\mathrm{\mu }\mathrm{Jy}{\mathrm{beam}}^{-1}$) and ${\sigma }_{\mathrm{rms},12″}=6.2\mathrm{\mu }{\mathrm{Jy}\mathrm{beam}}^{-1}$).

Fig. 7. Spectral index and uncertainty maps on the radio relic at 3″ resolution using the 1.5 GHz and 3.0 GHz observations. Radio contours are from the 3.0 GHz observations, drawn at the 2.5σrms × [1, 2, 4, 8, 16, 32] levels (with σrms = 2.7 μJy beam−1).

The radio halo showed hints of spectral index variation around the integrated value (i.e. ${\alpha }_{144\mathrm{MHz}}^{3.0\mathrm{GHz}}=-1.06\pm 0.03$; see also dashed line in the middle-right panel in Fig. 8). No significant curvature was found (see Fig. 6 and last bottom panel in Fig. 8).

Fig. 8. Spectral and curvature analysis of the relic in PSZ2G091. Left panel, from top to bottom: flux density, spectral index, curvature profiles across the relic and the halo at 12″ resolution. Purple circles refer to the relic and black squares to the halo. The box and sector numbers are shown in the right panel in Fig. 5; the black dashed line in the spectral index profile plot shows the integrated spectral index for the radio halo (i.e. ${\alpha }_{144\mathrm{MHz}}^{3.0\mathrm{GHz}}=-1.06\pm 0.03$). Right panel: Radio colour-colour plot for the radio relic at 4.5″ resolution. (See inset for the location of the boxes.) Arrows refer to the upper limit in the spectral index; the dotted line shows the power-law case across the spectrum; and the solid and dashed lines display the KGJP (with a particle injection time tinj = 100 Myr) and JP ageing models.

5. Discussion

At z = 0.822, PSZ2G091 is a remarkably peculiar cluster. Along with a central radio halo, this cluster hosts an elongated diffuse radio source oriented parallel to the merger axis, which in turn is inferred from the X-ray morphology. The radio morphology in total and polarised intensity, the spectral index, and the lack of a clear optical counterpart suggest that this source can be classified as a radio relic. The PSZ2G091 cluster is then the second-most massive galaxy cluster hosting a radio relic at z > 0.6¹⁰. In the following sections, we discuss the possible scenarios for the origin of such emission.

5.1. Radio spectral analysis

In the framework of the diffusive shock acceleration scenario, an initial population of electrons with a momentum (p) distribution of the type N(p)∝p^−δ, with δ = 1 − 2α, is (re-)energised by the passage of a merger-shock. As a consequence, the freshly (re-)accelerated particles will still follow a similar energy distribution¹¹ dN(E)/dE ∝ E^−δ_inj at the shock location, with a “flat” injected spectral index (α_inj), while the spectral index in the post-shock region becomes steeper due to a loss of energy caused by synchrotron and inverse Compton (IC) radiation.

Across the width of the relic, we observed a similar trend. We extracted the spectral index and curvature profiles using the 12″-resolution LOFAR, uGMRT (band 3 and 4), and VLA (L and S band) images from beam-size sectors covering the full relic (see light purple regions in the right panel in Fig. 5). Each sector corresponds to a physical size of ∼90 kpc at the cluster redshift. We found an increasing spectral steepening towards the cluster centre, as expected from the diffuse shock acceleration prediction, and a negative curvature in at least two of the outermost sectors (see left panel in Fig. 8). However, the curvature decreases downstream, contrary to what is expected in the post-shock region. This may be due to the mixing of different populations of electrons, such as those in the post-shock region, which lose energy due to synchrotron and IC radiation after (re-)energisation, and those in the halo region, which experience turbulent re-acceleration (Di Gennaro et al. 2021a). The presence of steepening and hints of curvature in the post-shock region was also observed in the higher-resolution¹² (i.e. 4.5″) colour-colour plots (Katz-Stone et al. 1993), that is, in the comparison of the spectral index in the “low” and “high” frequency bands (${\alpha }_{144\mathrm{MHz}}^{650\mathrm{MHz}}$ and ${\alpha }_{1.5\mathrm{GHz}}^{3.0\mathrm{GHz}}$, respectively; see right panel in Fig. 8). We obtained the “low-” and “high-frequency” spectral indices from beam-size boxes covering the full extension of the peripheral diffuse radio emission (see inset in the right panel in Fig. 8) and calculated the spectral indices as α = log(S₁/S₂)/log(ν₁/ν₂). In these kinds of plots, points laying on the line ${\alpha }_{144\mathrm{MHz}}^{650\mathrm{MHz}}={\alpha }_{1.5\mathrm{GHz}}^{3.0\mathrm{GHz}}$ follow a power-law spectrum, while points below it have convex spectra that are characterised by steeper spectral indices at high frequencies. We also overlayed two theoretical ageing models, the Jaffe-Perola (JP; Jaffe & Perola 1973) and the Komissarov-Gubanov (KGJP; Komissarov & Gubanov 1994) spectral ageing models, which take into account the injection spectral index, the magnetic fields, and, for the latter model, the time of injection, to describe the radiative losses. In this analysis, we assume α_inj = −0.9, B = 5 μG, and for the KGJP, t_inj = 100 Myr (which mimics projection effects). As for other relics (e.g., Di Gennaro et al. 2018; Rajpurohit et al. 2021c), the KGJP ageing model better represents the data.

If we assume that the radio relic traces a shock wave, we can relate the injected spectral index and the shock Mach number as

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathrm{radio},\mathrm{inj}}=\sqrt{\frac{2{\alpha }_{\mathrm{inj}}-3}{2{\alpha }_{\mathrm{inj}}+1}}.\end{array}$$(5)

From the spectral index profile, we measured a value in the outermost sector of α_inj = −0.89 ± 0.03, which results in a Mach number of ℳ_radio, inj = 2.48 ± 0.15. In the case that the particle cooling time is much shorter than the shock lifetime, we can relate the injected spectral index to the integrated spectral index (α_int) through the relation α_inj = α_int + 0.5 (Kardashev 1962). In this case, the shock Mach number is

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathrm{radio},\mathrm{int}}=\sqrt{\frac{{\alpha }_{\mathrm{int}}-1}{{\alpha }_{\mathrm{int}}+1}}.\end{array}$$(6)

For the relic in PSZ2G091, the injected spectral index obtained directly from the map is steeper than the one obtained with the integrated spectrum (∼ − 0.89 compared to −0.75), as also observed for other radio relics (e.g., “Sausage” relic, “Toothbrush” relic; Di Gennaro et al. 2018; Rajpurohit et al. 2018, respectively). The integrated spectral index of the relic, that is, α_int = −1.25 ± 0.3, corresponds to a Mach number of ℳ_radio, int = 3.0 ± 0.19.

5.2. X-ray and radio shock analysis

In order to reliably classify a patch of diffuse radio emission as a radio relic, it is necessary to detect a jump in the pressure (P) profile across the source (Markevitch & Vikhlinin 2007). In particular, to be defined as a shock discontinuity, we would need to measure p_post/p_pre > 1, where p = n_ek_BT (here, n_e is the electron density, k_B is the Boltzmann constant, and T is the ICM temperature). However, our X-ray observations are too shallow to provide clear temperature measurements in the pre- and post-shock regions, and only hints on the surface brightness profile could be derived.

We made use of the Chandra observations to investigate whether a surface brightness discontinuity can be found at the location of the relic. We modelled the underlying density distribution with a broken power-law (Markevitch & Vikhlinin 2007), assuming spherical symmetry, where the compression factor is defined as the ratio of the electron densities in the post- and pre-shock regions (i.e. 𝒞 = n_post/n_pre). We show in Fig. 9 the background-subtracted surface brightness profile across an elliptical sector centred at the coordinates of the northern sub-cluster and with an ellipticity matching the shape of the candidate radio relic (see inset in the figure). The observed profile was binned to a minimum S/N ratio of two. We estimated the residual background level caused by any imperfections in the blank sky approximation by using a “local background” region, namely, a 2′-radius circle located on the same CCD (the ACIS-I3 chip) but as far as possible (i.e. 4.7′) from the cluster. The local background count rate is statistically consistent with that obtained from the blank sky background; their difference is −3.9 ± 2.9 × 10⁻⁷ ph cm⁻² s⁻¹ arcmin⁻². This value was included as a fixed constant model in the surface brightness profile fitting. The fit was performed with proffit (Eckert 2016), and the model was evaluated using Cash statistics (Cash 1979). We detected a density jump $\mathcal{C}=2.{22}_{-0.30}^{+0.37}$ at the location of the candidate radio relic, that is, at r ∼ 1.2′ (i.e. r ∼ 550 kpc; stat/d.o.f. = 32.3/6). Assuming that the measured discontinuity can be associated with a shock wave, the X-ray Mach number can be calculated as:

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathrm{Xray}}=\sqrt{\frac{2\mathcal{C}}{\gamma +1-\mathcal{C}(\gamma -1)}},\end{array}$$(7)

Fig. 9. Surface brightness profiles across R2 with the best-fit model overlaid. The model corresponds to the emission from a projected broken power-law density distribution plus a constant surface brightness that accounts for the residual sky background. The data has been binned to a minimum S/N of two. On the bottom of the panel the residuals of the best-fit (i.e. (SX,obs − SX,mod)/ΔSX,obs) are displayed. In the inset, we show the location of the best-fit position of the X-ray discontinuity with respect to the radio relic (dashed line).

where γ = 5/3 is the adiabatic index of a monoatomic gas. Given the measured compression factor, we obtained an X-ray Mach number of ${\mathcal{M}}_{\mathrm{Xray}}=1.{93}_{-0.32}^{+0.42}$. We noticed these values are marginally smaller than those measured from the radio band. Domínguez-Fernández et al. (2021) showed that an initially uniform Mach number can result in a distribution when it encounters a turbulent medium such as the ICM (e.g., Zhuravleva et al. 2019) and that relics trace the high end of such a distribution. Moreover, Wittor et al. (2021) pointed out that viewing angles also affect such a distribution and that projection effects tend to impact X-ray observations more.

We checked that our results are robust to uncertainties in the background subtraction by fixing the residual background constant included in the model to (conservatively) ±20% of the total blank sky count rate. In these cases, the density jump varies in the range 𝒞 = 2.04 − 2.27 if we over- and under-subtract the background (stat/d.o.f. = 41.4/6 and stat/d.o.f. = 48.1/6, respectively). For a single power-law model, we obtained stat/d.o.f. = 88.9/8. Therefore, this model is disfavoured by the data.

5.3. Polarisation properties

In general, three polarisation models are used to describe the polarised emission (Burn 1966; Tribble 1991; Sokoloff et al. 1998): The “simple” case of a single Faraday screen, which only rotates the polarisation vectors (i.e. no depolarisation; ND); the case where an external Faraday screen (e.g., ICM) varies the magnetic fields on scales smaller than the restoring beam (i.e. external Faraday depolarisation; EFD); and the case where the Faraday screen is internal to the emitting source (i.e. internal Faraday depolarisation; IFD). These cases can be expressed as:

$$\begin{array}{c}\hfill P({\lambda }^2)=\{\begin{array}{c}{p}_{0}Iexp[2i({\chi }_0+\mathrm{RM}{\lambda }^2)](\mathrm{ND})\hfill \\ p_{0}Iexp(-2{\sigma }_{\mathrm{RM}}^2{\lambda }^4)exp[2i({\chi }_0+\mathrm{RM}{\lambda }^2)](\mathrm{EFD})\hfill \\ p_{0}I\left[\frac{1-exp(-2{\varsigma }_{\mathrm{RM}}^2{\lambda }^4)}{2{\varsigma }_{\mathrm{RM}}^2{\lambda }^4}\right]exp[2i({\chi }_0+\mathrm{RM}{\lambda }^2)](\mathrm{IFD}),\hfill \end{array}\end{array}$$(8)

with λ being the wavelength, p₀ the intrinsic polarisation fraction, χ₀ the intrinsic polarisation angle, RM, and σ_RM and ς_RM as the external and internal depolarisations, respectively.

In order to assess the cluster polarisation properties, we ran the QU-fitting code¹³ presented in Di Gennaro et al. (2021c; see also Appendix C) to obtain the intrinsic polarisation fraction (p₀), intrinsic polarisation angle (χ₀), RM, and external depolarisation (σ_RM). We created Stokes Q, U, and I cubes at 12.5″ resolution with 41 and 155 channels in the L and S band, respectively, and with δν = 8 MHz. The final frequency coverage was Δν = 1.34 − 4.0 GHz (i.e. between 0.075 m and 0.22 m), with an effective frequency of 3.0 GHz (λ = 0.1 m). These same single-channel images were also used to generate the total averaged polarised intensity and the Faraday spectrum through the pyrmsynth tool (see Table 5). We ran the fit for all the pixels above a given threshold, defined as starting from the ${\stackrel{\sim }{\sigma }}_{\mathrm{rms},P}$ value, that is, the root mean square noise of the total averaged polarised emission obtained from the RM-Synthesis technique, including the Ricean bias. We used $1.5{\stackrel{\sim }{\sigma }}_{\mathrm{rms},P}$ (with ${\stackrel{\sim }{\sigma }}_{\mathrm{rms},P}=23.4\mathrm{\mu }\mathrm{Jy}{\mathrm{beam}}^{-1}$) and $1.8{\stackrel{\sim }{\sigma }}_{\mathrm{rms},P}$ (with ${\stackrel{\sim }{\sigma }}_{\mathrm{rms},P}=9.9\mathrm{\mu }\mathrm{Jy}{\mathrm{beam}}^{-1}$) for the 1–4 GHz low-resolution (12.5″) and for the 2–4 GHz high-resolution (5″) images.

The maps of the polarisation parameters for the external depolarisation model at low resolution are displayed in Fig. 10. The intrinsic polarisation fraction is constant overall, around a value of 0.2, and the polarisation vectors appear to follow well the elongation of the relic. In the presence of both ordered and random magnetic fields, the theoretical maximum intrinsic polarisation fraction is set by the spectral index (Sokoloff et al. 1998; Govoni & Feretti 2004), namely,

$$\begin{array}{c}\hfill p_0=\frac{3\delta +3}{3\delta +7}\frac{1}{1+{\left(\frac{B_{\mathrm{rand}}}{B_{\mathrm{ord}}}\right)}^2},\end{array}$$(9)

Fig. 10. Polarisation parameter maps for PSZ2G091 at 12.5″ resolution. From left to right: intrinsic polarisation fraction (p0); intrinsic polarisation angle (χ0), with the polarisation vectors overlaid in white; RM, including the Galactic contribution; and external depolarisation (σRM). On top of each map, the grey histogram shows the distribution per pixel of the polarisation parameters.

with δ = 1 − 2α and B_ord and B_rand respectively as the “ordered” (i.e. aligned with the shock surface) and the “random” (i.e. isotropic) components of the magnetic fields. From our spectral index analysis, we obtained a maximum intrinsic polarisation fraction of 0.74. This is much higher than our measured value, implying that either the ratio B_rand/B_ord should be greater than one (Sokoloff et al. 1998; Govoni & Feretti 2004) or that the shock viewing angle should be smaller than 40° (Ensslin et al. 1998).

The RM shows a sharp transition from almost constant values around ∼50 − 60 rad m⁻² in R2_N to more mixed values in R2_S (RM ∼ 20 − 60 rad m⁻²; see Fig. 10 and Table 7). At the same location where we observed this RM change, the external depolarisation shifts from “low” (σ_RM ∼ 10 − 20 rad m⁻²) to “high” (σ_RM ∼ 60 − 70 rad m⁻²) values. Interestingly, we note that this transition happens where, in projection, the relic locates in a denser region of the ICM (see last panel in Figs. 3 and 4). This indicates that part of the denser ICM is in front of the relic and acting as a Faraday screen.

The RM of an extragalactic source at redshift z is given as:

$$\begin{array}{c}\hfill \mathrm{RM}=0.81{\displaystyle {\int }_{\mathrm{source}}^{\mathrm{observer}}}\frac{n_{\mathrm{e}}{B}_{\Vert }}{{(1+z)}^2}\mathrm{d}l[\mathrm{rad}{\mathrm{m}}^{-2}].\end{array}$$(10)

Here, n_e is the electron density (in cm⁻³), B_∥ is the magnetic field (in μGauss) along the line of sight, and l is the path length through the magneto-ionic medium (in pc). Since the observed polarisation signal is a combination of all that comes along the line of sight, the measured RM is actually the combination of different RM layers:

$$\begin{array}{c}\hfill \mathrm{RM}={\mathrm{RM}}_{\mathrm{Gal}}+{\mathrm{RM}}_{\mathrm{extragal}}+{\mathrm{RM}}_{\mathrm{cluster}}.\end{array}$$(11)

We used the latest Galactic Faraday Rotation catalogue (Hutschenreuter et al. 2022) and estimated the mean value of RM_Gal at the cluster location (l = 91.830, b = +26.116; see Table 1) within a radius of 0.1°¹⁴ (corresponding to ∼2.7 Mpc) to be RM_Gal = 60 ± 12 rad m⁻². This Galactic RM value is consistent with what we measured in R2_N, while the larger spread in R2_S suggests that there is also some contribution from the cluster. If we assume that RM_extragal = 0 (i.e. the ICM is the only Faraday screen), the mean cluster RM value¹⁵ is of a few rad m⁻². For a density column of n_e = 2 × 10⁻⁵ cm⁻³ (i.e. the value we obtain from the X-ray surface brightness profile; see Sect. 5.2) and for a path length of the magnetised plasma of L = 260 kpc (i.e. $L\approx 2\sqrt{2d_s{r}_s}$, where d_s = 15 kpc and r_s ∼ 500 kpc are the intrinsic width of the shock and its distance from the cluster centre, respectively; see Kierdorf et al. 2017), we obtained B_∥ ∼ 5 − 10 μG. These values are similar to low-z clusters (i.e. CIZA J2242.8+5301; see van Weeren et al. 2010; Di Gennaro et al. 2021c).

5.4. A potentially more complex merger event

The position of the radio relic in PSZ2G091 is unusual compared to other relics in the literature. An idealised binary merger in the plane of the sky would produce double radio relics perpendicular to the merger axis. The double radio relic in CIZA J2242.8+5301 is a textbook example of this (van Weeren et al. 2010; Hoang et al. 2017; Di Gennaro et al. 2018). Instead, the radio relic in PSZ2G091 is located eastward of the ICM distribution, parallel to the merger direction.

One possible explanation for the relic formation is that the two sub-clusters in PSZ2G091 are merging with a large impact parameter (i.e. b > 0). We compared our observations with simulations presented in van Weeren et al. (2011), which used the FLASH code (Fryxell et al. 2000) based on adaptive mesh refinement to follow the merger with high resolution but only in the places of interest, such as shock discontinuities (ZuHone 2010). According to their simulations, asymmetrical spiral-like shocks would develop in case of an impact factor of b = 4r_c, with r_c as the core radius, after about 1 Gyr from the collision (see their Fig. 2). This could represent the case of the radio relic in PSZ2G091, and the lack of the shock counterpart could be explained by a large mass ratio or by the presence of pre-existing plasma only at the R2 location being re-accelerated by the mild shock wave (van Weeren et al. 2017). This early merger scenario would also explain the flat spectral index of the radio halo and the small spectral fluctuations, as only the most turbulent regions would have had time to re-accelerate particles (t_acc < t_merger, see Donnert et al. 2013).

Another possibility is that the northern sub-cluster is itself undergoing a merger event. In other words, PSZ2G091 is a triple merger. In this scenario, the northern sub-cluster would have merged in the north-west, south-east direction, generating the radio relic, while the southern sub-cluster moves in the north, north-east direction towards the northern sub-cluster. Deep optical observations aiming to determine the spectroscopic redshifts of the cluster galaxy members could provide insights into this possible scenario.

Finally, we cannot completely exclude the possibility that the X-ray surface brightness discontinuity is associated with an equatorial shock and that the radio relic in PSZ2G091 is a consequence of that. However, the role of these types of shocks in accelerating particles is not clear. Moreover, the presence of the radio halo, which would require approximately a few Gyr to form after the merger event (e.g., Beresnyak & Miniati 2016), would disfavour this scenario.

5.5. Thermal and non-thermal correlation for the radio halo

Correlations between emission at different wavelengths from the same spatial regions can provide hints regarding common physical processes. Particularly, the emission from radio halos is expected to be related with that of the ICM (Brunetti & Jones 2014), according to the following relation:

$$\begin{array}{c}\hfill I_{\mathrm{rad}}\propto I_X^k,\end{array}$$(12)

with I_rad and I_X as the radio and the X-ray surface brightness, respectively. In order to investigate the existence of such a correlation at all the frequencies available, we performed point-to-point analyses (Govoni et al. 2001) using the Point-to-point TRend EXtractor¹⁶ (Ignesti 2022). At each frequency, we calculated I_rad and I_X from boxes of 12″ resolution with a radio surface brightness above 2.5σ_rms, ν at each frequency (see Table 6 for the map noises). This threshold also provided enough X-ray counts in each single box, which is the limiting factor due to the short exposure, and enough statistics to investigate the correlation. Other radio sources, such as the radio relic and the radio galaxies, were masked. We then fit the data using the BCES (bivariate correlated errors and intrinsic scatter) orthogonal method (Akritas & Bershady 1996).

The point-to-point analysis on radio halos commonly shows sub-linear trends (i.e. k < 1; Hoang et al. 2019; Xie et al. 2020; Rajpurohit et al. 2021b,a, 2023; Bonafede et al. 2022). In PSZ2G091, we find a trend consistent with a linear or super-linear correlation between the radio and the ICM distribution (i.e. k ≳ 1; see Fig. 11 and Table 8). Although the super-linear correlation would be in line with the early merger scenario, as the turbulence is not uniformly spread in the cluster volume, we point out that the combination of shallow X-ray observations, possible mixing of radio emission from electrons in the radio halo and in the radio relic, the large intrinsic scatter of the best-fit (i.e. σ_k = 0.14 − 0.20), and the large physical scale of the cells (i.e. ∼90 kpc), which hides possible small-scale fluctuations, does not allow for a clear investigation of the I_rad − I_X correlation. Finally, hints of a trend of k with the frequency were noted, similarly to what has been found by Rajpurohit et al. (2021c) and Hoang et al. (2021), in the radio halos in MACS J0717.5+3745 and CLG 0217+70, respectively.

Fig. 11. Point-to-point analysis of the radio halo in PSZ2G091 at 144 MHz, 400 MHz, 650 MHz, 1.5 GHz and 3.0 GHz at 12″ resolution. The inset of each panel shows the corresponding radio image with the beam-sized boxes used to compare the thermal and non-thermal emission.

6. Conclusions

In this paper, we present a multi-frequency analysis of the high-redshift galaxy cluster PSZ2 G091.83+26.11 (z = 0.822). We made use of LOFAR (120–168 MHz), uGMRT (250–500 MHz and 550–900 MHz), and VLA (1–2 GHz and 2–4 GHz) radio observations and data from the Chandra X-ray satellite to investigate the properties of the diffuse radio emission in the cluster, with a particular focus on the candidate radio relic. Below, we summarise the main results of our work:

– The diffuse radio emission in PSZ2 G091.83+26.11 is visible up to 3.0 GHz. Particularly, the radio halo extends in the north-east, south-west direction for about 1.2 Mpc, while at high resolution (∼2″) the bright extended relic-like source eastward of the cluster appears to be broken into two pieces, with lengths of about 640 kpc and 300 kpc and widths of about 15 kpc.

– The diffuse radio source eastward of the cluster (i.e. the candidate radio relic) has a flux density of 15.9 ± 0.6 mJy and 5.0 ± 0.2 mJy at 1.5 GHz and 3.0 GHz, respectively, at 12″ resolution. The resulting integrated spectral index, which also includes the LOFAR and uGMRT data, is α = −1.25 ± 0.03.

– The spectral index map and profile across the candidate radio relic steepens towards the cluster centre, from −0.89 ± 0.03 to −1.39 ± 0.03, with hints of non-negligible spectral curvature. Moreover, part of the source is also visible in polarisation, with polarisation vectors that follow the source surface. All these pieces of evidence suggest that this source can be classified as a radio relic. Using QU-fitting and assuming the only screens to rotate the polarisation vectors are the Galactic foreground and the ICM, we obtained an intrinsic polarisation fraction of about 20%. From the RM, we estimate a relic magnetic field of 5 − 10 μ Gauss.

– The X-ray Chandra analysis revealed a surface brightness jump at the location of the cluster elongated radio source (i.e. $\mathcal{C}=2.{22}_{-0.30}^{+0.37}$). Assuming this is due to shock compression, this corresponds to a Mach number $\mathcal{M}=1.{93}_{-0.32}^{+0.42}$.

– We compared our observational results with simulations based on adaptive mesh refinement, and we find that the location of the radio relic can be explained by an offset merger with impact factor b = 4r_c. However, we cannot entirely exclude a scenario also involving multiple merger events or equatorial shocks.

– The emission from the radio halo in PSZ2 G091.83+26.11 is in line with the previously published studies (Di Gennaro et al. 2021a,b). Its flux density is 6.1 ± 0.2 mJy and 3.1 ± 0.1 mJy at 1.5 GHz and 3.0 GHz, respectively, at 19″ resolution. Combining these values with the flux densities at LOFAR and uGMRT frequencies, we obtained an integrated spectral index of α = −1.06 ± 0.03. The I_rad − I_X point-to-point analysis suggests a trend consistent with a linear or super-linear correlation.

Given the above results – the offset position of the radio relic, its relatively short distance from the cluster centre (i.e. about 500 kpc) and low polarisation (i.e. ∼20%), the flat spectral index of the radio halo with no sign of curvature nor frequency break, and the cluster redshift – we speculate that PSZ2G091 could be in the early phase of its merger event.

Acknowledgments

We thank the referee for the suggestions which improved the quality of the manuscript. G.D.G. acknowledges support from the Alexander von Humboldt Foundation. M.B. acknowledges support from the Deutsche Forschungsgemeinschaft under Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306. R.J.v.W. acknowledges support from the ERC Starting Grant ClusterWeb 804208. A.S. is supported by the Women In Science Excel (WISE) programme of the Netherlands Organisation for Scientific Research (NWO), and acknowledges the Kavli IPMU for the continued hospitality. SRON Netherlands Institute for Space Research is supported financially by NWO. G.B. and R.C. acknowledge support from INAF through the mainstream project “Cluster science with LOFAR”. W.R.F. acknowledges support from the Smithsonian Institution, the Chandra High Resolution Camera Project through NASA contract NAS8-03060, and NASA Grants 80NSSC19K0116, GO1-22132X, and GO9-20109X. A.I. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 833824). H.J.A.R. acknowledges support from the ERC Advanced Investigator programme NewClusters 321271. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This paper is based on data obtained with the LOw Frequency Array (LOFAR). LOFAR (van Haarlem et al. 2013) is the Low Frequency Array designed and constructed by ASTRON. It has observing, data processing, and data storage facilities in several countries, which are owned by various parties (each with their own funding sources), and which are collectively operated by the ILT foundation under a joint scientific policy. The ILT resources have benefited from the following recent major funding sources: CNRS-INSU, Observatoire de Paris and Université d’Orléans, France; BMBF, MIWF-NRW, MPG, Germany; Science Foundation Ireland (SFI), Department of Business, Enterprise and Innovation (DBEI), Ireland; NWO, The Netherlands; The Science and Technology Facilities Council, UK; Ministry of Science and Higher Education, Poland; The Istituto Nazionale di Astrofisica (INAF), Italy. This research made use of the Dutch national e-infrastructure with support of the SURF Cooperative (e-infra 180169) and the LOFAR e-infra group. The Jülich LOFAR Long Term Archive and the German LOFAR network are both coordinated and operated by the Jülich Supercomputing Centre (JSC), and computing resources on the supercomputer JUWELS at JSC were provided by the Gauss Centre for Supercomputing e.V. (grant CHTB00) through the John von Neumann Institute for Computing (NIC). This research made use of the University of Hertfordshire high-performance computing facility and the LOFAR-UK computing facility located at the University of Hertfordshire and supported by STFC [ST/P000096/1], and of the Italian LOFAR IT computing infrastructure supported and operated by INAF, and by the Physics Department of Turin university (under an agreement with Consorzio Interuniversitario per la Fisica Spaziale) at the C3S Supercomputing Centre, Italy. This paper is based on data obtained with the Giant Metrewave Radio Telescope (GMRT). We thank the staff of the GMRT that made these observations possible. GMRT is run by the National Centre for Radio Astrophysics of the Tata Institute of Fundamental Research. This research has made use of data obtained from the Chandra Data Archive and the Chandra Source Catalog, and software provided by the Chandra X-ray Center (CXC) in the application packages CIAO and Sherpa. This research made use of APLpy, an open-source plotting package for Python (Robitaille & Bressert 2012). G.D.G. acknowledges Luca Di Mascolo for the help and support while dealing with APLpy.

References

Appendix A: Peeling on L-band data

In this section, we present the additional calibration performed on the point source at the location RA = 18^h31^m24.6^s DEC = +62°30′34.32″, whose side lobes were affecting the visibilities on the observing target. In Fig. A.1, we present the images of the target and peeled source after the “standard” self-cal calibration and after the on-source self-cal and bandpass calibrations. After the additional calibration, the source was subtracted from the visibilities.

Fig. A.1. Peeling results on the 1–2 GHz observations of PSZ2G091. Left panel: Cluster image after standard selfcal. Middle panel: Cluster image after standard peeling. Right panel: Cluster image after additional rounds of bandpass calibration on the troubling source. In the inset of each panel, we show the improvement on the calibration of the troubling source, with the limits on the colourmap set to the last step.

Appendix B: Additional frequency observations

In this section, we present the images for the additional frequency observations used for the analysis (Fig. B.1). Details on these images are reported in Table B.1. A comparison between all the observations used for this study at 12″ is shown in Fig. B.2.

Fig. B.1. Same as Fig. 2 but for the LOFAR 144 MHz and uGMRT 400 MHz and 650 MHz observations (top to bottom rows).

Fig. B.2. Comparison of the 12″ observations at each frequency. Radio contours start from 2.5σrms, ν, where σrms,144MHz = 137.1 μJy beam−1, σrms,400MHz = 98.7 μJy beam−1, σrms,650MHz = 41.6 μJy beam−1, σrms,1.5MHz = 17.4 μJy beam−1, and σrms,3.0GHz = 6.2 μJy beam−1.

Appendix C: Polarisation QU-fitting

In this section we briefly summarise the QU-fitting procedure to retrieve the polarisation parameters. We also present additional fitting plots for single high- and low-S/N pixels and for the full relic.

The Q(λ²) and U(λ²) emission was simultaneously fitted with cosine and sine models, while the I(λ²) emission was fitted with a log-parabolic model¹⁷ in order to take into account the curvature in the spectrum. We then used a Markov chain Monte Carlo (MCMC; Foreman-Mackey et al. 2013) approach in the fitting procedure to evaluate the posterior probability of the polarisation parameters (i.e. p₀, χ₀, RM, and eventually either σ_RM or ς_RM) and employed the following initial priors:

$$\begin{array}{c}\hfill \{\begin{array}{c}{I}_0\in [0,+\infty ]\hfill \\ A\in [-\infty ,+\infty ]\hfill \\ B\in [-\infty ,+\infty ]\hfill \\ p_0\in [0,1]\hfill \\ {\chi }_0\in [-\infty ,+\infty ]\hfill \\ \mathrm{RM}\in [-400,+400]\hfill \\ {\sigma }_{\mathrm{RM}}^2\in [0,+\infty ]\mathrm{or}{\varsigma }_{\mathrm{RM}}^2\in [0,+\infty ].\hfill \end{array}\end{array}$$(C.1)

We note that with this prior on χ₀, we needed to subsequently fold the polarisation angle to π to account for multiple angle rotations. Angle values of zero or π and ±π/2 reflect the north-south and east-west magnetic field directions, respectively.

In Fig. C.1, we show the corresponding negative (top) and positive (bottom) uncertainties maps of Fig. 10. In Fig. C.2, we show examples of fitting results on a single pixel on R2_N and R2_S, while in Fig. C.3, we show the integrated fitting results.

Fig. C.1. Corresponding negative (top) and positive (bottom) uncertainties maps of Fig. 10

Fig. C.2. QU-fitting results for a single pixel in R2N (top) and for a single pixel in R2S (bottom), assuming the external depolarisation model (EDF, Eq. 8). Left panel: Fits on Stokes I, Q, and U fluxes. Middle panel: Resulting fractional polarisation, p(λ2), and polarisation angle, χ(λ2). Right panel: Corner plot for the distribution of the uncertainties in the fitted polarisation parameters (i.e. p0, χ0, RM, and ${\sigma }_{\text{RM}}^2$). The contour levels are drawn at [0.5, 1.0, 1.5, 2.0]σ, with σ as the 68% statistical uncertainty (see dashed lines in the 1D histogram).

Fig. C.3. Same as Fig. C.2 but for R2N (top) and R2S integrated (bottom).

All Tables

All Figures

	Fig. 1. Composite image of PSZ2G091. Optical (white): PanSTARSS gri observations. Radio (red): 3.0 GHz VLA observations. X-ray (blue): Chandra 0.8–4.0 keV observations.
In the text

Fig. 2. VLA images of PSZ2G091 at 1.5 GHz (top row) and 3.0 GHz (bottom row). From left to right: Full-resolution (i.e. no taper applied) taper of 25 kpc, 50 kpc, and 100 kpc (see Table 3 for the final resolutions and map noise). The beam size is shown in the bottom-left corner of each panel. Labels in the leftmost panel follow those shown in Di Gennaro et al. (2021b). The white dashed circle shows the R = 0.5RSZ, 500 region, and the cluster centre is marked with a white cross. Radio contours, when shown, are drawn with white solid lines at the 2.5σrms × [1, 2, 4, 8, 16, 32] levels; a −2.5σrms level is also drawn with a dashed line.

In the text

Fig. 3. Zoom on the radio relic. From left to right: highest-resolution total intensity VLA image (i.e. 3.0 GHz at 1.8″ × 1.1″; see Table 3); total averaged polarisation intensity (without correction for Ricean bias, see Di Gennaro et al. 2021c) in the 2–4 GHz band (effective frequency 3.1 GHz) at 5″ resolution; high resolution-total intensity image at 3.0 GHz with the intrinsic polarisation magnetic field vectors at 5″ resolution, corrected for Faraday rotation; optical PanSTARRS gri image with radio contours from the leftmost panel starting from 2σrms; Chandra 0.8–4.0 keV image with radio contours from the leftmost panel starting from 2σrms.

In the text

Fig. 4. Thermal and non-thermal comparison. X-ray 0.8–4.0 keV Chandra image with 3.0 GHz radio contours tapered at 50 kpc (contour levels are the same as for Fig. 2). The radio beam size is shown in the bottom-left corner. The white dashed circle shows the R = 0.5RSZ, 500 region, and the white cross denotes the cluster centre. The two yellow stars show the cluster X-ray peaks, which are associated with the two candidate sub-clusters (Artis et al. 2022).

In the text

Fig. 5. Integrated spectra of the diffuse radio emission in PSZ2G091. Left: the flux densities were measured from the solid-line regions in the right panel, and the spectral indices are reported in the bottom-left corner in the panel (black for the halo and shades of purple for the relic). Right: Radio image and contours are from the 3.0 GHz VLA image at 12″ resolution. The separation of R1 and R2 is indicated by the purple solid line in the right panel. The dashed dashed-line sectors refer to the profiles displayed in Fig. 8. White circles mark the position of the compact sources that have been masked for the flux density, spectral and curvature profiles.

In the text

Fig. 6. Spectral index and curvature maps (top) with corresponding uncertainty maps (bottom) of PSZ2G091. Left column: resolution of 4.5″ using the 144 MHz, 650 MHz, 1.5 GHz, and 3.0 GHz observations. Central and Right columns: Resolution of 12″ using the 144 MHz, 400 MHz, 650 MHz, 1.5 GHz, and 3.0 GHz observations. Radio contours are from the 3.0 GHz observations, drawn at the 2.5σrms × [1, 2, 4, 8, 16, 32] levels (with ${\sigma }_{\mathrm{rms},4.5″}=3.1\mathrm{\mu }\mathrm{Jy}{\mathrm{beam}}^{-1}$) and ${\sigma }_{\mathrm{rms},12″}=6.2\mathrm{\mu }{\mathrm{Jy}\mathrm{beam}}^{-1}$).

In the text

	Fig. 7. Spectral index and uncertainty maps on the radio relic at 3″ resolution using the 1.5 GHz and 3.0 GHz observations. Radio contours are from the 3.0 GHz observations, drawn at the 2.5σrms × [1, 2, 4, 8, 16, 32] levels (with σrms = 2.7 μJy beam−1).
In the text

Fig. 8. Spectral and curvature analysis of the relic in PSZ2G091. Left panel, from top to bottom: flux density, spectral index, curvature profiles across the relic and the halo at 12″ resolution. Purple circles refer to the relic and black squares to the halo. The box and sector numbers are shown in the right panel in Fig. 5; the black dashed line in the spectral index profile plot shows the integrated spectral index for the radio halo (i.e. ${\alpha }_{144\mathrm{MHz}}^{3.0\mathrm{GHz}}=-1.06\pm 0.03$). Right panel: Radio colour-colour plot for the radio relic at 4.5″ resolution. (See inset for the location of the boxes.) Arrows refer to the upper limit in the spectral index; the dotted line shows the power-law case across the spectrum; and the solid and dashed lines display the KGJP (with a particle injection time tinj = 100 Myr) and JP ageing models.

In the text

Fig. 9. Surface brightness profiles across R2 with the best-fit model overlaid. The model corresponds to the emission from a projected broken power-law density distribution plus a constant surface brightness that accounts for the residual sky background. The data has been binned to a minimum S/N of two. On the bottom of the panel the residuals of the best-fit (i.e. (SX,obs − SX,mod)/ΔSX,obs) are displayed. In the inset, we show the location of the best-fit position of the X-ray discontinuity with respect to the radio relic (dashed line).

In the text

	Fig. 10. Polarisation parameter maps for PSZ2G091 at 12.5″ resolution. From left to right: intrinsic polarisation fraction (p0); intrinsic polarisation angle (χ0), with the polarisation vectors overlaid in white; RM, including the Galactic contribution; and external depolarisation (σRM). On top of each map, the grey histogram shows the distribution per pixel of the polarisation parameters.
In the text

	Fig. 11. Point-to-point analysis of the radio halo in PSZ2G091 at 144 MHz, 400 MHz, 650 MHz, 1.5 GHz and 3.0 GHz at 12″ resolution. The inset of each panel shows the corresponding radio image with the beam-sized boxes used to compare the thermal and non-thermal emission.
In the text

Fig. A.1. Peeling results on the 1–2 GHz observations of PSZ2G091. Left panel: Cluster image after standard selfcal. Middle panel: Cluster image after standard peeling. Right panel: Cluster image after additional rounds of bandpass calibration on the troubling source. In the inset of each panel, we show the improvement on the calibration of the troubling source, with the limits on the colourmap set to the last step.

In the text

	Fig. B.1. Same as Fig. 2 but for the LOFAR 144 MHz and uGMRT 400 MHz and 650 MHz observations (top to bottom rows).
In the text

	Fig. B.2. Comparison of the 12″ observations at each frequency. Radio contours start from 2.5σrms, ν, where σrms,144MHz = 137.1 μJy beam−1, σrms,400MHz = 98.7 μJy beam−1, σrms,650MHz = 41.6 μJy beam−1, σrms,1.5MHz = 17.4 μJy beam−1, and σrms,3.0GHz = 6.2 μJy beam−1.
In the text

	Fig. C.1. Corresponding negative (top) and positive (bottom) uncertainties maps of Fig. 10
In the text

Fig. C.2. QU-fitting results for a single pixel in R2N (top) and for a single pixel in R2S (bottom), assuming the external depolarisation model (EDF, Eq. 8). Left panel: Fits on Stokes I, Q, and U fluxes. Middle panel: Resulting fractional polarisation, p(λ2), and polarisation angle, χ(λ2). Right panel: Corner plot for the distribution of the uncertainties in the fitted polarisation parameters (i.e. p0, χ0, RM, and ${\sigma }_{\text{RM}}^2$). The contour levels are drawn at [0.5, 1.0, 1.5, 2.0]σ, with σ as the 68% statistical uncertainty (see dashed lines in the 1D histogram).

In the text

	Fig. C.3. Same as Fig. C.2 but for R2N (top) and R2S integrated (bottom).
In the text