© H. W. Edler et al. 2022

1. Introduction

Galaxy clusters are the building blocks of the large-scale structures of the Universe. When observed at radio wavelengths, we can detect synchrotron emission that traces the nonthermal plasma in the cluster volume (for a review, we refer to van Weeren et al. 2019). Different phenomena are associated with this emission. There is consensus that diffuse cluster-scale sources such as radio halos develop in the intracluster medium (ICM) as a consequence of mergers between clusters (Brunetti & Jones 2014). Radio halos are megaparsec-scale sources with a steep spectral index of α ≲ −1.2¹ that are located in the central region of clusters, approximately following the thermal ICM brightness distribution as observed in the X-rays. The turbulent reacceleration model (Brunetti et al. 2001; Petrosian 2001; Brunetti & Lazarian 2007, 2016; Miniati 2015) is the favored scenario to explain the mechanism that generates radio halos. In this model, turbulence injected in the ICM by a merger event reaccelerates cosmic-ray electrons (CRe) from a population of at least mildly relativistic primary (e.g., accelerated by active galactic nuclei) or secondary (generated in proton-proton collisions) seed electrons (Brunetti & Lazarian 2011; Pinzke et al. 2017). This is supported by the observed connection between the dynamic state of galaxy clusters and radio halos (Cassano et al. 2010a, 2013; Cuciti et al. 2015, 2021). Another important expectation of this model is a cutoff in the electron spectrum that depends on the energetics of the merger event and translates into a gradual steepening in the observed synchrotron spectrum. This implies that at a frequency of ∼1 GHz, we would in general be able to detect radio halos that were generated during very energetic merger events, while at lower frequency, we might be sensitive to a population of radio halos that is characterized by very steep radio spectra (α < −1.5), which are in general referred to as ultra-steep spectrum radio halos (USSRH; Cassano et al. 2006; Brunetti et al. 2008). Such sources have been observed in only a small number of clusters so far (Brunetti et al. 2008; Macario et al. 2011; van Weeren et al. 2012; Wilber et al. 2018; Duchesne et al. 2021; di Gennaro et al. 2021).

Other extended radio sources in galaxy clusters, such as gently reenergized tails (GReETs; de Gasperin et al. 2017) and radio phoenixes (Kale & Dwarakanath 2011; Cohen & Clarke 2011; Duchesne et al. 2021), are thought to trace ancient tails and lobes of radio galaxies that have been reaccelerated. They are smaller sources (∼100 kpc) and are directly connected to individual radio galaxies. The observational difference between radio phoenixes and GReETs lies in the morphology, the spectral properties, and in the presence of a shock in the ICM. These differences can be explained by a difference in the underlying reacceleration mechanism. Phoenixes are mostly irregular sources with a steep spectrum without clear spatial trends, and they are thought to be accelerated by the adiabatic compression of old radio lobes through a shock wave (Enßlin & Gopal-Krishna 2001; Enßlin & Brüggen 2002). In contrast, GReETs are elongated tails of radio galaxies with highly unusual properties. Starting from the host radio galaxy, spectral aging causes an increasing steepening of the observed synchrotron emission, but at some point, the tail shows an unexpected rebrightening that coincides with a constant or even flattening spectral index. This is thought to be explained by a gentle energetization mechanism that is probably connected to microturbulence in the tail that is induced by interactions with the ICM (de Gasperin et al. 2017; van Weeren et al. 2021). The exact nature of the mechanism remains unclear because only a small number of such sources or candidates have been reported in the literature so far, with very limited spectral information (de Gasperin et al. 2017; Cuciti et al. 2018; Wilber et al. 2018; Botteon et al. 2021a; Ignesti et al. 2022).

Observations at low radio frequencies (< 1 GHz) are essential for studying ultra-steep sources such as USSRHs and GReETs because these sources remain undetectable at higher frequencies, which is due to their spectral properties. The LOw-Frequency Array (LOFAR; Van Haarlem et al. 2013) is the largest and most sensitive radio-interferometer operating at low frequencies. It is therefore perfectly suited for studying diffuse radio emission in clusters of galaxies.

Abell 1033 (PSZ2 G189.31+59.24, hereafter A1033) is a moderately massive ($M_{500}=3.{24}_{-0.32}^{+0.30}\times {10}^{14}{M}_{\odot }$) galaxy cluster at a redshift of z = 0.126 (Planck Collaboration XXVII 2016) that shows recent merger activity and hosts a steep spectrum source previously classified as a radio phoenix (de Gasperin et al. 2015). In addition, the cluster contains an elongated ultra-steep (α ≈ −4.0) source connected to a wide-angle tail radio galaxy (WAT). Gentle reacceleration has been suggested to explain the peculiar spectral properties (de Gasperin et al. 2017). This source is the prototype of the GReET category, and it remains the most extreme of the reported cases of radio galaxy tails that show signs of reacceleration. Observations in the ultra-low frequency regime (< 100 MHz) are necessary for a more detailed study of this source. Furthermore, a radio halo in the cluster has recently been reported in Botteon et al. (2022) based on 144 MHz LOFAR observations.

We analyze a LOFAR low-band antenna (LBA) observation of A1033 and data of the LOFAR Two-Meter Sky-Survey (LoTSS, Shimwell et al. 2017, 2019, 2022) to study the properties and spectral behavior of the peculiar radio sources of this cluster. Throughout this work, we assume a fiducial flat ΛCDM cosmology with Ω_m = 0.3 and H₀ = 70 km s⁻¹ Mpc⁻¹. All spatial distances are given with respect to a redshift of z = 0.126, at which one arcsecond corresponds to 2.53 kpc.

2. Observations and data reduction

For the ultra-low frequency study of A1033, we reduced and analyzed a dedicated LOFAR LBA observation (project LC10-008, PI: de Gasperin) between 30 and 78 MHz. Furthermore, we used data of two 8 h LOFAR HBA pointings from the LOFAR Two-metre Sky Survey (Shimwell et al. 2017, 2019, 2022) centered at 144 MHz, as well as archival data of a 5.5 h GMRT observation at 323 MHz and of a 64 ks Chandra X-ray observation in the energy band 0.5−7 keV. The data reduction and analysis of the Chandra observation is described in de Gasperin et al. (2015), and for the reduction of the GMRT data we refer to de Gasperin et al. (2017). We summarize information about the radio observations in Table 1.

2.1. LOFAR LBA data

The LOFAR LBA observation was conducted in observation mode LBA_OUTER, where the outer half of the station dipoles are used to minimize electromagnetic crosstalk. During the preprocessing, the data were flagged at a resolution of 1 s and 3 kHz to mitigate radio frequency interference. The data were subsequently averaged in time and frequency to 4 s and 49 kHz. The multibeam capability of LOFAR LBA allowed it to continuously point one beam toward the calibrator source 3C196 during the observation. This source is a well-behaved standard calibrator for the low band (Heald et al. 2015). The calibrator data were reduced with the LOFAR LBA calibrator pipeline described in de Gasperin et al. (2019). This pipeline is used to find solutions for station-based direction-independent effects as well as an initial estimate for the ionospheric phase errors. The calibrator solutions revealed a poor data quality for the station CS031LBA, and therefore, we excluded this station from further processing. The calibrator solutions were applied to the target field dataset. Additionally, we corrected the data for the primary beam toward the phase center. Next, we performed direction-independent self-calibration of the target field to find solutions for the average ionospheric effects. We employed the self-calibration pipeline presented in de Gasperin et al. (2020). Starting from an initial model based on a collection of radio surveys (Intema et al. 2017, TGSS; Condon et al. 1998, NVSS; Rengelink et al. 1997, WENNS; Lane et al. 2014, VLSSr), solutions for the direction-independent ionospheric total electron content (TEC), Faraday rotation, and second-order primary beam effect were derived in two rounds of self-calibration. The root mean square (rms) background noise of the self-calibrated image is 2.5 mJy beam⁻¹ at a resolution of 35″, the image quality is limited by significant direction-dependent ionospheric errors still present in the data. These errors were addressed in direction-dependent calibration, where we employed the calibration pipeline that is being developed for the LOFAR LBA Sky Survey (LoLSS; de Gasperin et al. 2021) and based on the facet calibration strategy outlined in van Weeren et al. (2016). For the next processing steps, we reduced the data volume by averaging the data to a resolution of 8 s in time. We used the direction-independently calibrated image to isolate suitable direction-dependent calibrators. We employed the Python blob detector and source finder (PyBDSF; Mohan & Rafferty 2015) as source finder, and selected bright and compact sources by thresholding based on flux density and source area. Compact sources in a proximity of 6′ at most were merged by a grouping algorithm. We selected all sources and groups of sources with a flux density above S_min > 1.0 Jy × (ν/60 MHz)^−0.8 at ν = 30 MHz or ν = 54 MHz. The reason for this is that spectral variation, mainly of the primary beam, can make it hard to identify all bright sources at a single frequency. This yielded 29 calibrator directions for the first of two major iterations of our strategy. To prepare the data for further calibration, we subtracted the model obtained after direction-independent calibration to create a dataset that was empty up to model inaccuracies. Then, we looped through the calibrators ordered by their flux density at 30 MHz and repeated the following steps. First, the model of the calibrator was readded to the empty data. Second, the data were phase shifted to the calibrator direction and subsequently averaged further to a resolution of 16 s (if S₅₄ > 10 Jy) or 32 s (otherwise) in time and 0.39 MHz in frequency. This averaging was possible because we only imaged the single calibrator and not the whole field of view (FoV), so that we were not limited by time or frequency smearing. Third, we corrected this averaged dataset for the primary beam in the new phase center. Then, we self-calibrated on this dataset: By imaging a small field around the calibrator, we obtained a model of the calibrator source(s). The data were smoothed in time and frequency with a baseline-dependent kernel size. Here, shorter baselines were smoothed with a larger kernel to increase the signal-to-noise ratio. We solved on the smoothed data for scalar phases against the model. The frequency interval of this solve step was 0.39 MHz, the time interval varied from 64/128 s in the first iteration down to 16/32 s in the following minor cycles. These solutions were then applied to the data, and the calibrator direction was imaged again to obtain an improved model for the next self-calibration iteration. For sources brighter than 10 Jy, we performed two additional diagonal amplitude calibration steps: one step that enforced constant solutions between stations to capture beam model errors, and a subsequent step solved to correct for remaining errors. If the rms noise of the image did not improve by more than 1% during a self-calibration iteration and the noise level was lower than the initial value, we considered this direction to be converged. Otherwise, it was not further used as a calibrator direction. In this case, the solutions were discarded and the direction was divided between the neighboring calibrators. We resubtracted the converged calibrator model that was corrupted with the final solutions in that direction from the data to obtain a cleaner empty dataset. Artifacts in the image around the calibrator related to direction-dependent errors were now strongly reduced. This procedure was repeated for all 29 calibrator directions.

We then started again with the averaged data after direction-independent calibration. From these data, we peeled all calibrator sources outside of a 2.3° region around the pointing center, meaning that we subtracted the model after corrupting it with the corresponding direction-dependent solutions without readding it to the data at any point. This strongly reduced the impact of bright sources outside of the primary beam. These data were then imaged using DDFacet (Tasse et al. 2018). This imager simultaneously applied calibration solutions for the remaining 22 calibrator directions based on Voronoi-tesselated facets around the calibrators. The resulting image and model are the products of the first major iteration of our direction-dependent calibration pipeline. We performed a second major iteration using the improved model as input. This time, we searched for calibrators on the DDE-calibrated image and went 20% lower (to 0.8 Jy) in terms of minimum calibrator flux. We found 25 suitable directions for the second iteration, which were now all inside the primary beam main lobe. Again we looped over these sources ordered according to their brightness and repeated the individual self-calibration cycles as described above. Since A1033 is a somewhat special source because it is both very bright (> 10 Jy at 54 MHz) and extended (around 5′), our default algorithm had problems in correctly detecting all emission associated with A1033 as a calibrator direction. For this target, we therefore employd an LBA-specific implementation of the extraction-pipeline described in van Weeren et al. (2021). This pipeline subtracted all sources outside of a 0.2° radius around A1033 from the data using the model and calibration solutions from the direction-dependent pipeline. Then, it performed multiple rounds of self-calibration solving for TEC and a constant phase in each time interval. This tecandphase-constraint reduces the number of free parameters compared to the phase-solve used in direction-dependent calibration. We imaged the final target-extracted dataset using a robust weighting of −0.3 and −1.0. The images have a background noise of 1.4 and 1.7 mJy beam⁻¹ at a resolution of 21″ × 13″ and 12″ × 9″. They are displayed in Fig. 1.

Fig. 1. Four images (a)–(d) showing A1033 at different frequencies and resolutions. Contours start from 3σ and increase in powers of two. Dashed contours correspond to −3σ. The white ellipse in the bottom left corner of each image corresponds to the beam size as reported in Table 2. Image (d) was produced using archival GMRT data from de Gasperin et al. (2017). (a) LBA 54 MHZ, 21″ × 14″ resolution. (b) LBA 54 MHZ, 12″ × 9″ resolution. (c) HBA 144 MHZ, 12″ × 6″ resolution. (d) GMRT 330 MHZ, 10″ × 6″ resolution.

2.2. LOFAR HBA data

A1033 was covered by the pointings P155+27 and P157+35 observed in the context of LoTSS (Shimwell et al. 2017, 2019, 2022). The analysis of P157+35 has been presented in de Gasperin et al. (2017). Here, we exploit both datasets, taking advantage of the latest developments in the processing pipelines, to increase the sensitivity of the images toward the cluster. LoTSS is observing the entire northern sky with 8 h observations per pointing, covering the frequency range 120−168 MHz. LoTSS pointings are processed by the LOFAR Surveys Key Science Project using the automated direction-independent and -dependent calibration pipelines described in van Weeren et al. (2016), Williams et al. (2016), de Gasperin et al. (2019), and Tasse et al. (2018, 2021). A1033 is located within the newly released area of LoTSS Data Release 2 (DR2; Shimwell et al. 2022), and because a cluster in the PSZ2 catalog (Planck Collaboration XXVII 2016), it is also included in the sample of PSZ2/LoTSS-DR2 clusters recently analyzed by Botteon et al. (2022). Here, we used the reprocessed dataset presented in Botteon et al. (2022), in which the data quality has been improved using the extraction and self-calibration procedure described in van Weeren et al. (2021). We used WSCLEAN (Offringa et al. 2014; Offringa & Smirnov 2016) to reimage the data and study the properties of the diffuse emission in the cluster. The final HBA image is shown in Fig. 1c. It has a resolution of 12 × 6, and with a noise level of 110 μJy beam⁻¹, it is about 40% deeper than the image previously presented in de Gasperin et al. (2017).

3. Results

In Table 2 we summarize the medium- and high-resolution LBA as well as the HBA and GMRT images of A1033 displayed in Fig. 1. In Fig. 2, the radio sources in the cluster are labeled. A connection between the WAT radio galaxy east of the cluster and the GReET is visible at 54 and 144 MHz. A second radio galaxy with a head-tail morphology is marked HT in the figure.

Fig. 2. Color-composite of the high-resolution LOFAR HBA (orange), the Chandra X-ray (blue) and the SDSS g (white, Alam et al. 2015) images. The white line corresponds to the path that is used to analyze the spectral properties of the GReET, and the dots mark the distance from the starting point in increments of 100 kpc. Paths A (green) and B (yellow) highlight the location of possible interrupted structures. The dashed section shows the discontinuity.

For our analysis, we assumed a systematic flux scale uncertainty of 10% for the LOFAR LBA, LOFAR HBA, and the GMRT (de Gasperin et al. 2021; Shimwell et al. 2022; Chandra et al. 2004). For all radio images, a common lower uv-cut of 100 λ was used during imaging, and all imaging was carried out using WSCLEAN. When information from multiple frequencies is combined for spectral index mapping or spectral shape-modeling, a number of preprocessing steps are required to ensure that the images are comparable. For each image, a point-source catalog was created using PyBDSF, and the images were aligned by cross-matching point sources to account for potential astrometric errors at a level below 5″. Next, all images were convolved to the smallest common circular beam. Last, they were regridded on the same pixel layout.

3.1. Spectral index maps

We created spectral index maps in the ranges 54−144 MHz and 144−323 MHz. For this, we discarded all pixels with a surface brightness significance below 3σ_rms in at least one of the two relevant images. We used the linear least-squares method in log-space to determine the spectral index values. The uncertainties were estimated according to Gaussian propagation of uncertainty. To investigate the spectral shape, we created a curvature map from the spectral index maps. The curvature is defined as ${\mathrm{\Delta }}_{\alpha }={\alpha }_{54}^{144}-{\alpha }_{144}^{325}$. The high-resolution spectral index and curvature maps are displayed in Fig. 3, and the corresponding uncertainty maps are provided in Fig. A.1. At the position of the WAT host galaxy, we measure a spectral index of α = −0.65 in both frequency intervals. Along the tail, the spectral index gradually steepens. The steepening is stronger at higher frequency. At the location of the bend in the tail, the spectral index does not steepen further and remains approximately constant along the GReET for both frequency ranges. For this section of the tail, we measure an extremely steep spectral index of −3.8 ± 0.18 between 144 and 325 MHz, which is fairly uniform for a projected distance of 300 kpc. Toward ultra-low frequencies, the spectrum of this region shows extreme curvature (Δ_α = 1.8 ± 0.23).

Fig. 3. High-resolution spectral index map of A1033. Left image: spectral index between 54 and 144 MHz, and right image: same between 144 and 323 MHz (center). The images are on the same color scale and have an angular resolution of 12″ × 9″. Upper limits below −1.6 are outside the black contours and marked by arrows. The saturation scales inversely with the uncertainty on the spectral index. Right image: corresponding curvature map, defined as the difference ${\mathrm{\Delta }}_{\alpha }={\alpha }_{54}^{144}-{\alpha }_{144}^{325}$.

3.2. Spectral aging analysis

In absence of acceleration mechanisms, the electron population in a radio galaxy tail is subject to aging due to synchrotron and inverse Compton losses. This causes a spatial evolution of the spectral shape along the tail depending on the magnetic field, injection history, and injection spectral index of the active galactic nucleus (AGN) as well as the projected velocity of the radio galaxy. In the following, we consider this pure-aging scenario for the WAT and examine whether it agrees with our multifrequency observations. This is the expected scenario when no acceleration mechanisms, usually related to ICM shocks, are present in the tail. For our spectral aging model, we probed a path along the brightest part of the tail where we measured the flux densities at 54, 144, and 323 MHz and the spectral indices between 54−144 and 144−323 MHz in 23 different 12.5″ beam-sized regions that were separated by the beam FWHM. This path is shown in Fig. 2, and the flux density and spectral index evolution as a function of distance is shown in Fig. 4. We considered the Jaffe-Perola (JP) plasma aging model (Jaffe & Perola 1973), which is the most commonly used model for the aging of radio galaxies, to describe the spectral distribution for the first 300 kpc of the tail, before the brightness of the tail increases again (de Gasperin et al. 2017). We performed a nonlinear least-squares fit of the JP model to the spectral index points, assuming that the projected velocity v_⊥ of the galaxy along the path is constant.

Fig. 4. Flux density (top) and spectral index (bottom) measured in the 12.5″ beam-sized regions shown in Fig. 2 as a function of distance and estimated age. The continuous lines show the best-fitting model (χ2/d.o.f. = 29.5/17) assuming minimum aging (B = 2.3 μG) with a projected velocity of 780 km s−1. The dotted lines show an extrapolation of this model.

From the spectral index maps, we estimate the injection index to be α_inj = −0.65. The aging history of our model then depends on the magnetic field B and the velocity v_⊥. Unfortunately, the magnetic field and the projected velocity are strongly correlated, implying that it is not possible to distinguish between a scenario with a magnetic field close to the minimum aging field strength and a lower velocity and a scenario with more rapid losses and a higher velocity of the WAT galaxy (as illustrated in Fig. B.1). Therefore, we considered a scenario that maximizes the lifetime of the electrons, which provides an upper limit on the radiative age of the plasma. The lifetime of the CRe in the relevant energy regime depends on the magnetic field B as ${\tau }_{\mathrm{age}}\propto \sqrt{B}/(B^2+B_{\mathrm{CMB}}^2)$ (Stroe et al. 2014), where B_CMB = 3.2 × (1 + z)² μG (Longair 2011) is the magnetic field strength equivalent to the cosmic microwave background (CMB) energy density at redshift z. The expression for τ_age takes its maximum value at the minimum loss magnetic field B = B_min, where the loss rates due to synchrotron radiation and inverse Compton scattering with the CMB are equal. This is given by $B_{\min }=B_{\mathrm{CMB}}/\sqrt{3}$ at the redshift of A1033, B_min = 2.3 μG. We used the nonlinear least-squares method to fit a JP model to the spectral index data points (see Appendix B). By fitting the model to the spectral index points and not the flux densities, we can eliminate the normalization factors and thereby reduce dependences on the injection history of the AGN. The best-fitting model has a projected velocity of v_⊥ = 780 km s⁻¹ with a reduced χ² statistic of χ²/d.o.f. = 29.5/17, which agrees reasonably well with the ≈730 km s⁻¹ found in de Gasperin et al. (2017). This corresponds to a maximum age of the GReET of 800 Myr. The spectral index and flux density data points as well as the best-fitting model are shown in Fig. 4. For the flux density models, the flux normalization was fit in log-space independently as a free parameter for each beam-sized region.

To complement the spectral analysis, we created a color-color map. To do this, we calculated the 54−144 and 144−323 MHz spectral indices in 5″ square boxes in a region covering the WAT and GReET and a second region covering the radio phoenix. We considered all boxes within the two regions with a significance above 3σ_rms at all three frequencies. For comparison, we also show the spectral trajectory of a power-law injection and the JP aging model discussed above. For points along the WAT and GReET, we additionally color-code the distance from the injection point along the tail. The resulting diagram is shown in Fig. 5. The boxes associated with the radio phoenix occupy a rather confined region in the spectral index space without clear trends. In contrast, along the WAT, there is a clear steepening with increased distance from the injection point. The spectral index points departure from the PL-injection spectrum along a trajectory that agrees with a JP model with increased distance. After 300 kpc, the steepening between 54−144 and 144−323 MHz stops and the boxes corresponding to the following part of the GReET scatter around ${\alpha }_{144}^{323}=-4$ and ${\alpha }_{54}^{144}=-2$.

Fig. 5. Color-color plot of A1033. Each square represents a 5″ × 5″ region of the WAT/GReET (orange edges) or phoenix (blue). For the WAT/GReeT area, the square faces are color-coded according to their approximate distance from the injection point. Additionally, two lines show the spectral properties for a power-law injection (dotted gray) and for a JP aging model with B = 2.3 μG (dashed green).

3.3. Source-subtracted low-resolution images

The medium-resolution images Figs. 1a and c show an extended diffuse radio emission superimposed on the bright sources associated with individual radio galaxies, the GReET and the radio phoenix. In the LOFAR images, emission with an extension of several 100 kpc is visible at 54 MHz as well as at 144 MHz. This emission is associated with the radio halo detected in A1033 (Botteon et al. 2022). To isolate the emission of the radio halo from the foreground sources, we obtained a model for the compact sources based on high-resolution (Briggs = −1.0) images at 54 MHz, 144 MHz, and 323 MHz. This model is restricted to the clean components in a region that only includes the sources in the high-resolution images (e.g., the radio galaxies, GReET, and phoenix), and subtracted from the data, identical regions are used at all three frequencies. We reimaged the subtracted data using baselines in the range 100−20 000 λ and tapering the imaging weights to a resolution of 25″ to enhance the extended emission. The resulting images are displayed in Fig. 6 together with the regions we used to restrict our subtraction- model. In the LOFAR maps at 54 and 144 MHz in Fig. 6, the extended emission is significantly detected at a surface brightness above 5σ_rms for a region with a diameter of ≈700 kpc except for the subtracted regions. The shape is approximately circular in both frequencies, but some morphological differences exist. At 54 MHz, the halo appears to extend further to the southwest, while a hole in the emission is visible in the northeast. Given the complex procedure of the calibration of a bright, extended source with complex morphology superimposed on the faint, even more extended halo emission together with the source subtraction, we cannot exclude the possibility that these differences in morphology in the lower-significance regime are not physical, but rather artifacts in the data. At 323 MHz, the detected emission appears more patchy, probably because of artifacts caused by calibration errors or radio frequency interference in the GMRT map. In the low-resolution images, we also display the point-source-subtracted Chandra X-ray contours from de Gasperin et al. (2015). The shape of the radio emission is approximately circular and mostly confined within the X-ray contours.

Fig. 6. Source-subtracted low-resolution images. The image weights are tapered to a resolution of 25″. Red contours indicate sources that were subtracted from the data, and cyan contours are smoothed Chandra X-ray contours in the energy band 0.5−7 keV (de Gasperin et al. 2015), linearly spaced between 10−100σ, where 1σ = 6.6 × 10−6 photons s−1 cm−2 arcmin−2. The purple rectangle labeled P marks the brightest part of the visible halo. (a) 54 MHz, (b) 144 MHZ, (c) 325 MHZ.

3.4. Calculating the radio halo flux

To estimate the total flux density of the halo, we employed the halo flux density calculator Halo-FDCA (Boxelaar et al. 2021), an MCMC code for estimating halo flux densities with reduced human bias. The code can handle masks to exclude image regions from the analysis. This allowed us to exclude the subtracted regions that do not contain reliable emission. We used the code to fit a spherically symmetric exponential surface-brightness profile described by

$$\begin{array}{c}\hfill I(r)=I_0{e}^{-\frac{r}{r_{\mathrm{e}}}}\end{array}$$(1)

to the data. Here, I(r) is the surface brightness at a distance r from the profile center, I₀ is the peak surface brightness, and r_e is the e-folding radius (Murgia et al. 2009). This empirical model has frequently been found to describe radio halos accurately (e.g., Vacca et al. 2014; Osinga et al. 2021; Botteon et al. 2021a,b) and it comes with a minimum number of free parameters, which is beneficial for the halo in A1033, which is partly obscured by other sources. A similar analysis for A1033 in LOFAR HBA is presented in Botteon et al. (2022). To allow a consistent analysis of the LBA and HBA data specifically tailored to A1033, we repeated the source subtraction, imaging, and fitting procedures also for the HBA data.

The result of the halo flux calculation is a 1.46 ± 0.15 Jy at 54 MHz and 0.29 ± 0.03 Jy at 144 MHz. The detailed parameters are listed in Table 3, and plots of the halo fit result and fit residual are shown in Fig. C.1. The location and extension of the best-fitting models at 54 MHz and 144 MHz agree well with e-folding radii of 154 ± 2 and 165 ± 3 kpc, respectively. The LBA model is slightly farther south and the radius is smaller by 7%. We also manually measured the flux density within the respective 3σ_rms contours F_3σ of the halo in Fig. C.1, obtaining $F_{54\mathrm{MHz}}^{3\sigma }=0.87\pm 0.09$ and $F_{144\mathrm{MHz}}^{3\sigma }=0.15\pm 0.02$. Due to the missing extrapolation of the emission in the masked regions, the manually measured flux is a lower limit for the total flux of the halo. The ratio of the manually measured and the fitted flux density is different between LBA and HBA (0.6 and 0.52). We attribute this to a small difference in the location and size of the fitted surface brightness profile. We estimate the total radio power at 150 MHz of the halo based on the Halo-FDCA results. The resulting value is P_150 MHz = 1.22 ± 0.13 × 10²⁵ W Hz⁻¹, where a spectral index of α = −1.65 was used for extrapolation and k-correction, as we justify in the following Sect. 3.5. The independent estimate of the radio power in Botteon et al. (2022) (P₁₅₀ = 1.09 ± 0.28 × 10²⁵ W Hz⁻¹) is consistent with ours (at about 1σ), but their reconstructed morphology differs significantly (r_e = 165 ± 3 to r_e = 103 ± 1). We interpret this to mean that the halo morphology is more susceptible to differences in masking and fitting than the flux density of the halo. In general, the fitted parameters of the halo model are subject to a systematic uncertainty connected to the masking and subtraction procedure, which is hard to assess and was therefore not accounted for in our uncertainties.

3.5. Spectral index of the radio halo

We used the source-subtracted images presented in Fig. 6 to create maps of the low-resolution spectral index and to study the spatial variation of the spectral index within the radio halo. The spectral index maps were created as described in Sect. 3.1. We excluded the regions that were used to subtract the compact sources. These maps are displayed in Fig. 7, and the uncertainty maps are shown in Fig. A.2. In the range of 54−144 MHz, the spectral index of the halo could be determined significantly across a large area. The northern part of the halo is less steep with a spectral index of about α = −1.5, and the southwestern region is steeper than α = −2. The same trend extends to the ultra-steep upper limits in the southwest and to the flat-spectrum features expressed by lower limits in the northeast. As stated in Sect. 3.3, we cannot exclude the possibility that these features are at least partly a consequence of systematic errors in the calibration. Because only a part of the radio halo is detected above a significance of 3σ in the 323 MHz GMRT map, we can constrain the spectral index in the range 144−323 MHz only in a comparatively small area. In this region northwest of the GReET, the spectral index is ${\alpha }_{144}^{323}=-2.15\pm 0.19$. For the majority of the remaining halo area, we can constrain the spectrum to be steeper than α = −2. This suggests that the spectrum of the halo is curved.

Fig. 7. Source-subtracted spectral index maps at a spatial resolution of 30″. Top panel: spectral index between 54 MHz and 144 MHz, and bottom panel: same between 144 MHz and 325 MHz. The 3σ upper and lower limits are indicated by arrows. Dashed red lines highlight the regions from which sources were subtracted.

To estimate the integrated halo spectral index between 54 and 144 MHz, we compared different methods. Using the halo fits, we derive a spectral index of α_FDCA = −1.65 ± 0.17 based on the integrated surface brightness. If the flux density is measured manually in the region where we have > 5σ significance at both frequencies, the corresponding integrated spectral index is ${\alpha }_{\mathrm{man}}^{5\sigma }=-1.74\pm 0.14$. A more conservative estimation using only the area above 10σ_rms in both images yields ${\alpha }_{\mathrm{man}}^{10\sigma }=-1.69\pm 0.15$. We adopted the value based on the Halo-FDCA fluxes as reference value for the halo because it does not depend on the sensitivity of the observations.

4. Discussion

4.1. Spectral study of the GReET

In Fig. 4 we show the flux density and spectral index evolution along the tail together with the best-fitting JP model as described in Sect. 3.2. After about 200 kpc from the injection point, the spectral index appears to be systematically flatter than the best-fitting model. To some degree, this deviation might be caused by a change in velocity of the WAT. Presumably, the galaxy reached a maximum velocity at the minimum of its trajectory in the cluster potential well and slowed down afterward. This is not accounted for in our constant-velocity assumption. However, at about ≈350 kpc form the injection point, the spectrum does not steepen further and a plateau in spectral index associated with the GReET is reached. This constant and even flattening spectral index cannot be explained with the extrapolation of a pure aging model, even when a change in the WAT velocity is considered. The same stalling in the spectral index trend can also be observed in the color-color diagram in Fig. 5, where after initial JP-like aging, the spectral evolution stops and appears frozen for several 100 kpc. Consequently, an additional energetization mechanism is required. This mechanism needs to be just efficient enough to compensate for the radiation losses, but not so strong that the spectral shape is significantly flattened.

This was previously concluded in de Gasperin et al. (2017), where turbulent reacceleration of the electrons in the tail by magnetic pumping was proposed as a possible mechanism. The suggested scenario is that turbulence is forced into the tail by interactions with the surrounding ICM, which can provide mild acceleration that acts for sufficiently long timescales to explain the homogeneous nature of the GReET. The extreme curvature we detect fully supports the scenario outlined in the previous work: it is a consequence of a cutoff in the electron energy spectrum caused by the balance between the very gentle acceleration and the cooling of CRe, with an acceleration time of the radio-emitting electrons comparable to their cooling time (e.g., Brunetti & Jones 2014).

4.2. Turbulent velocity

In the high-resolution map shown in Fig. 2, the tail of the WAT shows disrupted features. We suggest that these can be attributed to large-scale shear flows and turbulent motions in the ICM. The tail of the WAT appears to be split into three filaments after about 250 kpc from the injection point. The brightest filament directly connects to the GReET, and the fainter filaments, labeled in Fig. 2 as A and B, are discontinuous. Because the GReET is clearly bifurcated, where one of the two parts has no clear connection to the WAT, we consider the possibility that a shear motion in the ICM, at a scale close to the injection scale of the turbulence, caused either A or B to disconnect from the southern segment of the GReET. This then allows us to constrain the large-scale turbulent velocity in the ICM: In the diffusion regime, CRs are displaced a distance d after a time t according to $d\sim \sqrt{\mathit{tD}}$, where D is the spatial diffusion coefficient. For super-Alfvénic turbulence, we can approximate D ∼ v_LL, where v_L is the turbulent velocity at scale L. If the displacement d ∼ L, the velocity of large-scale motions can be estimated as v_L ∼ d/t, where t can be constrained from the age of the plasma from our limit on the velocity of the WAT. We estimate the age at the discontinuities of path A and B to be t_A < 380 Myr and t_B < 310 Myr using the distance of the discontinuities A and B from the injection point (300 kpc and 250 kpc, respectively) as well as our limit on the WAT velocity of ≈780 km s⁻¹. The corresponding projected length scales of the discontinuity are d_A = 65 kpc and d_B = 144 kpc. This means that we can constrain the turbulent velocity to v_A ≥ 167 km s⁻¹ and v_B ≥ 464 km s⁻¹. These values fall within the typical range found for large-scale turbulent motions in simulations of cluster mergers, which are several 100 km s⁻¹ (Vazza et al. 2018). This velocity can also be compared to the sound speed c_s derived from the ICM temperature $k{T}_{\mathrm{ICM}}=6.{15}_{-0.88}^{+1.14}$ keV found in the Chandra X-ray analysis of de Gasperin et al. (2015), which yields

$$\begin{array}{c}\hfill c_{\mathrm{s}}=1480{\left(\frac{T_{\mathrm{ICM}}}{{10}^8\mathrm{K}}\right)}^{1/2}{\mathrm{km}\mathrm{s}}^{-1}={1250}_{-92}^{+111}{\mathrm{km}\mathrm{s}}^{-1}.\end{array}$$(2)

This implies that the turbulent motions estimated above are likely subsonic, in agreement with expectations in the ICM (e.g., Brunetti & Jones 2014).

This consideration raises the question of how the tail can resist turbulent diffusion for a timescale of several 100 Myr. When we assume a Kolmogorov cascade, the turbulent velocity at a scale L depends on the size of the scale as v_L ∝ L^1/3. Thus, at a scale equal to the thickness of the filaments, the velocity derived from our estimate should be v_A(10 kpc)≥89 km s⁻¹ or v_B(10 kpc)≥190 km s⁻¹. If this value is below the Alfvén velocity v_A, which is typically in the range ∼100 km s⁻¹ (Brunetti & Jones 2014), Reynolds and Maxwell stresses of magnetic fields may stabilize the narrow tails. Similar considerations have been put forward to explain the stability of radio filaments in AGN bubbles (Brienza et al. 2021). From the Navier-Stokes equation, the term corresponding to matter motions on a scale k: ρdv²k competes against magnetic field line tension on the same scale: B²k/(4π). Thus, for the line tension to dominate, we require

$$\begin{array}{c}\hfill B\ge v\sqrt{4\pi \rho }.\end{array}$$(3)

For a reference density 100 kpc form the cluster center (at the projected distance to the GReET) of ρ ≈ 10⁻²³ kg m⁻³, which is a typical value for clusters similar to A1033, and a velocity in the range derived above from v_A(10 kpc)≥89 km s⁻¹ to v_B(10 kpc)≥190 km s⁻¹, a magnetic field of at least 3.1−6.7 μG is required. For the upper half of this range, the magnetic field in the GReET would be greater than what is usually considered for magnetic fields in the ICM even for significantly more massive clusters (e.g., Bonafede et al. 2010), which would imply that it must be seeded from the WAT and sustained for the life time of the GReET. Lower values for the magnetic field in the tail are possible if the distance between the cluster center and the GReET is significantly larger than the projected distance, such that the density ρ entering Eq. (3) is lower.

Furthermore, we tried to quantify the conditions under which the electrons can be confined in the GReET for 800 Myr: In general, diffusion perpendicular to magnetic filaments is small, which allows an efficient trapping of the electrons. However, if the magnetic field in the GReET is turbulent, that is, sub-Alfvénic as we require to explain the stability, perpendicular diffusion of electrons may be driven by stochastic diffusion of field lines. The perpendicular diffusion length of the CRe $l_{\perp ,\text{cr}}^2$ after a time t can be estimated as (Lazarian & Yan 2014)

$$\begin{array}{c}\hfill l_{\perp ,\mathrm{cr}}^2\sim \frac{{(Dt)}^{3/2}}{27L}{M}_{\mathrm{A}}^4,\end{array}$$(4)

where in this case, the Alfvénic Mach-number M_A ∼ 1 and the injection scale L ∼ L_⊥ = 10 kpc is the traverse scale of the GReET. We define the lifetime of the electrons in the GReET as the time when the perpendicular diffusion is comparable to the traverse size l_⊥, cr ∼ L_⊥,

$$\begin{array}{c}\hfill t\sim 9L_{\perp }^2/D.\end{array}$$(5)

To be stable for 800 Myr, a diffusion coefficient below D = 3.4 × 10²⁹ cm² s⁻¹ is required, which agrees well with the previous estimate of de Gasperin et al. (2017) and observational results for other clusters based on studies of the metallicity profiles (Rebusco et al. 2005, 2006).

4.3. Radio halo

The spectral index of α = −1.65 ± 0.17 between 54 and 144 MHz places the radio halo in A1033 in the USSRH category. Together with the further steepening of the spectrum toward higher frequencies, this agrees with a cutoff in the electron spectrum as predicted by the turbulent acceleration model (Brunetti et al. 2001, 2008; Petrosian 2001; Kuo et al. 2003; Cassano et al. 2006). It is well known that radio halos in more massive clusters are more powerful (Cassano et al. 2013; van Weeren et al. 2021; Cuciti et al. 2021). Because the mass of A1033 is considerably lower than that of the vast majority of clusters that are known to host a radio halo, the cluster is expected to have a lower radio power. In Fig. 8 we compare the 150 MHz radio power of the A1033 halo to the radio halos and to the correlations of halo power to cluster mass presented in van Weeren et al. (2021). Furthermore, we show the two halos detected in the LOFAR deep fields (Osinga et al. 2021) and the low-power halos discovered in Hoang et al. (2021) and Botteon et al. (2021b). All these radio halos are detected close to 150 MHz and are hosted by clusters in the PSZ2 catalog. This comparison shows that while A1033 is one of the lowest-mass clusters hosting a radio halo, the halo is quite powerful. It lies considerably above the low-frequency correlations. Compared to the 150 MHz correlations from van Weeren et al. (2021), it is overluminous by a factor in the range of 7−55. This is much larger than the intrinsic scatter of the correlations, which is ∼3.

Fig. 8. Radio power at 150 MHz of low-frequency radio halos as a function of the PSZ2 mass. The lines indicate the extreme cases of the correlations of radio power to cluster mass reported in van Weeren et al. (2021), which were derived using different fitting methods and samples. The dark blue line shows the steepest (orthogonal method) and the light blue line shows the flattest (including candidate halos, Y|X method) correlation. The shaded regions show the 95% confidence interval. The halos reported as ‘van Weeren+21 & lit.’ include a literature radio halo sample (de Gasperin et al. 2015; Cordey 1985; Wilber et al. 2018; Botteon et al. 2019; Bîrzan et al. 2019; Savini et al. 2018, 2019; Bonafede et al. 2018; Hoang et al. 2019a,b; Macario et al. 2013; George et al. 2017; Duchesne et al. 2021).

The relatively large radio power is notable given the ultra-steep spectrum nature of the A1033 halo because in general, USSRH are expected to be less powerful (Cassano et al. 2010b; Cuciti et al. 2021). We investigated whether the power estimate of the halo might be contaminated by other sources. We found that while the region we used to subtract the GReET indeed contains the significant emission in the HBA high-resolution map, images at lower resolution indicate a less confined feature with a ten times lower surface brightness in extension of the GReET. This patch of emission is coincident with the brightest region of the unobscured part of the radio halo, labeled P in Fig. 6. To ensure that patch P does not contaminate our radio power estimate significantly in case it is not related to the halo, we repeated the flux density calculation, this time with a more conservative mask. The resulting radio power estimate is lower by only 1.3%. Consequently, contamination from embedded sources cannot explain the large observed power.

A radio halo with a power above the correlation should correspond to an energetic merger event and/or a merger in an evolutionary phase close to the peak power. In both cases, the X-ray surface brightness distribution should show strong signs of disruption. The X-ray morphological parameters for the cluster were determined in the Chandra analysis in de Gasperin et al. (2015). While the X-ray concentration parameter c (defined, e.g., in Cassano et al. 2013) does not unambiguously qualify the cluster as a merging one (c = 0.200 ± 0.004), the centroid shift w = 0.086 ± 0.006, the high differential velocity of the brightest cluster galaxy (BCG), and a distinct bimodality in the redshift distribution indicate a clear and possibly strong merger scenario (de Gasperin et al. 2015). The strong merger might not fully show in the X-ray morphology parameters because it is likely occurring with a significant line-of-sight component, as also concluded in de Gasperin et al. (2015) based on the BCG dynamics and redshift bimodality. Cuciti et al. (2021) found a trend of the distance from the correlation of radio power to cluster mass and the X-ray morphological disturbance. Following their definition of the disturbance d_X-ray, we have d_X-ray = 0.86 ± 0.03, which places A1033 in the region of disturbed clusters. While a more energetic merger should lead to a flatter spectral index, the relatively low cluster mass might be the dominant factor determining the spectral properties. An increasing number of spectral studies of halos in low- and intermediate-mass clusters will help to constrain the spectral properties of halos in this mass regime.

A further reason for the comparatively high power and the steep spectrum of the halo in A1033 could lie in a greater abundance of seed electrons in the cluster, possibly accumulated by the various other bright radio sources. The synchrotron luminosity in reacceleration models is given by (e.g., Brunetti & Vazza 2020; di Gennaro et al. 2021)

$$\begin{array}{c}\hfill P\propto F\eta \frac{B^2}{B^2+B_{\mathrm{cmb}}^2},\end{array}$$(6)

where η is the (re-)acceleration efficiency, and $F=\rho v_{L_{\text{inj}}}^3/L_{\text{inj}}$ is the turbulent energy flux in the emitting volume, which depends on the ICM density ρ as well as the injection scale L_inj and turbulent velocity v_Linj. The efficiency η is the fraction of the turbulent energy flux contributing to CRe acceleration and proportional to the CRe energy density u_CR: η ∝ u_CR (Brunetti & Lazarian 2007; Bonafede et al. 2022). Following from Eq. (6), for the same CRe energy spectrum, a higher electron density would accordingly increase the halo power while not causing a spectral flattening. We speculate that at the low surface brightness patch in extension of the GReET (marked P in Fig. 6) with spectral properties in between those of the GReET and the halo, we might witness electrons that were conserved by the GReET entering the CRe population associated with the halo. A similar connection was proposed in Wilber et al. (2018), for example.

Because only a very small number of halos are detected in clusters with masses similar to or lower than A1033, the shape and scatter of the correlation in this mass range is not yet well determined. The detailed analysis of the largest sample of radio halos so far, based on the second data release of LoTSS (Botteon et al. 2022), will provide strongly improved constraints in the near future (Cuciti et al., in prep.).

4.4. Radio phoenix

The bright extended source 150 kpc south of the cluster center was previously studied in de Gasperin et al. (2015) and was categorized as a radio phoenix. A connection to a bright elliptical galaxy (marked S in Fig. 2) that coincides with a radio point source at 1.4 GHz was proposed. The radio phoenix has an irregular and elongated shape with a largest linear size of ∼380 kpc. Toward the south, it shows a steep gradient, while toward the cluster center, the source is fading more slowly, and at lower frequencies, it extends toward the head-tail galaxy and the BCG. In Fig. 9, its flux density is compared to the head-tail radio galaxy and the GReET between 54 and 1400 MHz. The source shows spectral curvature. The spectral index flattens from ${\alpha }_{609}^{1425}=-1.58\pm 0.17$ (de Gasperin et al. 2015) to ${\alpha }_{54}^{144}=-0.99\pm 0.14$. The spatial variation of the spectral properties (Fig. 3) does not reveal a clear trend. The spectral curvature we observe is a characteristic feature of radio phoenixes (Enßlin & Gopal-Krishna 2001). Therefore, the ultra-low frequency picture of the source is in line with the phoenix classification. The source differs from the GReET in the spectral shape, its spectrum is significantly less curved between 54 and 323 MHz (Δ_α ≈ 0.3 to Δ_α ≈ 2.0), and furthermore, its morphology is more irregular and not clearly connected to a radio galaxy.

Fig. 9. Flux densities as a function of frequency. Additional GMRT and VLA data at 609 MHz and 1.4 GHz are taken from de Gasperin et al. (2015, 2017).

5. Conclusion

We presented an analysis of the galaxy cluster Abell 1033 based on new ultra-low frequency data taken with LOFAR LBA at 54 MHz and on existing data at higher frequencies (LOFAR HBA data at 144 MHz and GMRT data at 323 MHz). The cluster is an especially interesting target because despite its comparatively low mass, it hosts a variety of strong and extended radio sources such as a GReET, a radio phoenix, interacting AGN, and a radio halo, which show very heterogeneous properties. Our findings are listed below.

• For the GReET, we detected extreme spectral curvature. The spectrum steepens form α ≈ −4 between 144 and 323 MHz to α ≈ −2 between 54 and 144 MHz. This finding is in line with the emission being generated by strongly aged electrons that are reaccelerated by a mechanism with a very low acceleration efficiency.

• Assuming a maximum lifetime scenario, we found a lower limit for the projected velocity of the WAT of 780 km s⁻¹; this agrees with earlier findings. This corresponds to an maximum age of the tail of ∼800 Myr.

• We investigated the properties of the radio halo in A1033. We found a radio power of P_150 MHz ≈ 1.22 × 10²⁵ W Hz⁻¹. This power is higher than the correlations of radio power to cluster mass by a factor > 7. At the same time, we found the spectrum of the halo to be ultra-steep between 54 and 144 MHz (α ≈ −1.69), with a further steepening at higher frequencies.

• We presented possible reasons for both the high luminosity and the ultra-steep spectrum of the halo: The halo power might be explained by an energetic merger that does not fully reveal itself in the X-ray morphology parameters due to a significant line-of-sight component of the merger axis. This scenario is supported by the redshift bimodality and the high radial velocity of the BCG. This requires that the spectral properties are dominantly controlled by the cluster mass. Alternatively, a particularly high density of CRes, possibly accumulated by the various bright radio sources in the cluster, can explain an increased radio power without spectral flattening. Simulations could be used to investigate whether a significant per-cluster variation of seed electrons can be caused by sources such as the GReET and radio phoenix.

• We detected two candidate disrupted filaments in the GReET. We speculated that large-scale turbulent motions are responsible for the disruption. With the lower limit on the WAT velocity we derived, we constrained the turbulent velocity to be greater than 167 km s⁻¹ and 464 km s⁻¹, which is within the typical range found in simulations. The thin tail itself could be stabilized by magnetic fields if the turbulence at these small scales becomes sub-Alfvénic.

• We found that the spectrum of the radio phoenix flattens from α ≈ −1.6 to α ≈ −1.0 between 1.4 GHz and 54 MHz. This which is a characteristic property of radio phoenixes.

With the further advance of low-frequency radio surveys such as LoTSS and LoLSS, it will be possible to perform spectral analyses for a rapidly growing number of cluster radio sources. These will show whether the GReET in A1033 is representative for a greater category of sources or if it is truly special in terms of its low acceleration efficiency. Additionally, systematical analyses of radio halos in clusters of low and moderate mass will lead to a better understanding of the halo power and the spectral index distributions in this mass range.

Acknowledgments

We thank the referee for the useful comments which helped to improve the manuscript. LOFAR 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 NWO (grant 2019.056). 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. The data are published via the SURF Data Repository service which is supported by the EU funded DICE project (H2020-INFRAEOSC-2018-2020 under Grant Agreement no. 101017207). This project is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under project number 427771150. F.d.G. acknowledges support from the Deutsche Forschungsgemeinschaft under Germany’s excellence strategy – EXEC2121 “Quantum Universe” – 390833306. V.C. acknowledges support from the Alexander von Humboldt Foundation. R.J.v.W. acknowledges support from the ERC Starting Grant ClusterWeb 804208. A.B. acknowledges support from the VIDI research program with project number 639.042.729, which is financed by the Netherlands Organisation for Scientific Research (NWO), and from the ERC-StG DRANOEL n. 714245. A.D. acknowledges support by the BMBF Verbundforschung under the grant 05A20STA.

References

Appendix A: Uncertainty maps

Fig. A.1. Uncertainty maps to Figure 3. The left and center figure show the spectral index uncertainties between 54-144 MHz (left) and 144-323 MHz (center). The right image shows the uncertainty of the curvature map.

Fig. A.2. Source-subtracted spectral index uncertainty maps at a spatial resolution of 30″. The two figures are for the frequency ranges 54-144 MHz (left) and 144-323 MHz (right). Dashed red lines highlight the regions from which sources are subtracted.

Appendix B: Spectral aging model

We employed the standard JP spectral aging model to describe the aging in the first section of the WAT radio galaxy. The flux density of this model at a frequency ν and after a time t is given by (Harwood et al. 2013)

$$\begin{array}{c}\hfill S_{\mathrm{JP}}(N_0,\nu ,B,t,z,{\alpha }_{\mathrm{inj}})=N_0\frac{\sqrt{3}ln(10)e^{3}B}{16\pi {ϵ}_{0}c{m}_{\mathrm{e}}{(z+1)}^2}{\displaystyle {\int }_0^{\pi }\mathrm{d}\delta sin{(\delta )}^2\int \mathrm{d}log(E)EF(x)n_{\mathrm{e}}(E,B,t,z,{\alpha }_{\mathrm{inj}}),}\end{array}$$(B.1)

where e and m_e are the electron charge and mass, ϵ₀ is the vacuum permittivity, δ is the pitch angle, and F(x) is the following function:

$$\begin{array}{c}\hfill F(x)=x{\displaystyle {\int }_x^{\infty }{K}_{5/3}(z)\mathrm{d}z.}\end{array}$$(B.2)

The variable x = ν/ν_c is the ratio of the frequency ν and the critical frequency ν_c, which is given by

$$\begin{array}{c}\hfill {\nu }_{\mathrm{c}}=\frac{3E^{2}eBsin(\delta )}{4\pi m_{\mathrm{e}}^3{c}^4}·\end{array}$$(B.3)

In Equation B.1, the energy-integration is performed in log-space for numerical efficiency. We first fit a spectral index model to the observed spectral index values α to remove the dependence on the normalizations N₀. The spectral index model is calculated as

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{JP},i}({\nu }_1,{\nu }_2,B,t,z,{\alpha }_{\mathrm{inj}})=\frac{log\frac{S_{\mathrm{JP}}({\nu }_1,B,t,z,{\alpha }_{\mathrm{inj}})}{S_{\mathrm{JP}}({\nu }_2,B,t,z,{\alpha }_{\mathrm{inj}})}}{log\frac{{\nu }_1}{{\nu }_2}}·\end{array}$$(B.4)

The uncertainty for the observed spectral index values σ_α is obtained from the uncertainties of the flux densities σ_S according to Gaussian propagation of uncertainty,

$$\begin{array}{c}\hfill {\sigma }_{\alpha }=\frac{1}{log({\nu }_1/{\nu }_2)}\sqrt{\frac{{\sigma }_{S_1,\mathrm{stat}.}^2+{(0.1\times S_1)}^2}{S_1^2}+\frac{{\sigma }_{S_2,\mathrm{stat}.}^2+{(0.1\times S_2)}^2}{S_2^2}}.\end{array}$$(B.5)

Here we take into account the 10% systematic uncertainty on the flux scales. We fit the WAT projected velocity v_⊥ by minimizing the standard χ² statistics assuming a minimum aging magnetic field B_min and a linear motion t = d_⊥/v_⊥,

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{i=0}^9\sum _{j=0}^1{\left(\frac{{\alpha }_i({\nu }_j,{\nu }_{j+1})-{\alpha }_{\mathrm{JP},i}({\nu }_j,{\nu }_{j+1},B_{\min },d_{\perp }/v_{\perp },z,{\alpha }_{\mathrm{inj}})}{{\sigma }_{{\alpha }_i}({\nu }_j,{\nu }_{j+1})}\right)}^2.}\end{array}$$(B.6)

Subsequently, a normalization factor N_0, i is fitted in log-space for each beam-sized region along the WAT/GReET.

Fig. B.1. Aging model merit function against magnetic field (x-axis) and projected velocity (y-axis) in arbitrary units.

Appendix C: Halo fit results

Fig. C.1. Results of the halo fitting in LBA (a) and HBA (b). The circles in the left panels indicate the position of the best-fit halo on top of the flux density maps. The center panels show the model, and the right panels present the residual. The green and red regions are excluded from the fitting.

All Tables

All Figures

Fig. 1. Four images (a)–(d) showing A1033 at different frequencies and resolutions. Contours start from 3σ and increase in powers of two. Dashed contours correspond to −3σ. The white ellipse in the bottom left corner of each image corresponds to the beam size as reported in Table 2. Image (d) was produced using archival GMRT data from de Gasperin et al. (2017). (a) LBA 54 MHZ, 21″ × 14″ resolution. (b) LBA 54 MHZ, 12″ × 9″ resolution. (c) HBA 144 MHZ, 12″ × 6″ resolution. (d) GMRT 330 MHZ, 10″ × 6″ resolution.

In the text

Fig. 2. Color-composite of the high-resolution LOFAR HBA (orange), the Chandra X-ray (blue) and the SDSS g (white, Alam et al. 2015) images. The white line corresponds to the path that is used to analyze the spectral properties of the GReET, and the dots mark the distance from the starting point in increments of 100 kpc. Paths A (green) and B (yellow) highlight the location of possible interrupted structures. The dashed section shows the discontinuity.

In the text

Fig. 3. High-resolution spectral index map of A1033. Left image: spectral index between 54 and 144 MHz, and right image: same between 144 and 323 MHz (center). The images are on the same color scale and have an angular resolution of 12″ × 9″. Upper limits below −1.6 are outside the black contours and marked by arrows. The saturation scales inversely with the uncertainty on the spectral index. Right image: corresponding curvature map, defined as the difference ${\mathrm{\Delta }}_{\alpha }={\alpha }_{54}^{144}-{\alpha }_{144}^{325}$.

In the text

	Fig. 4. Flux density (top) and spectral index (bottom) measured in the 12.5″ beam-sized regions shown in Fig. 2 as a function of distance and estimated age. The continuous lines show the best-fitting model (χ2/d.o.f. = 29.5/17) assuming minimum aging (B = 2.3 μG) with a projected velocity of 780 km s−1. The dotted lines show an extrapolation of this model.
In the text

	Fig. 5. Color-color plot of A1033. Each square represents a 5″ × 5″ region of the WAT/GReET (orange edges) or phoenix (blue). For the WAT/GReeT area, the square faces are color-coded according to their approximate distance from the injection point. Additionally, two lines show the spectral properties for a power-law injection (dotted gray) and for a JP aging model with B = 2.3 μG (dashed green).
In the text

Fig. 6. Source-subtracted low-resolution images. The image weights are tapered to a resolution of 25″. Red contours indicate sources that were subtracted from the data, and cyan contours are smoothed Chandra X-ray contours in the energy band 0.5−7 keV (de Gasperin et al. 2015), linearly spaced between 10−100σ, where 1σ = 6.6 × 10−6 photons s−1 cm−2 arcmin−2. The purple rectangle labeled P marks the brightest part of the visible halo. (a) 54 MHz, (b) 144 MHZ, (c) 325 MHZ.

In the text

	Fig. 7. Source-subtracted spectral index maps at a spatial resolution of 30″. Top panel: spectral index between 54 MHz and 144 MHz, and bottom panel: same between 144 MHz and 325 MHz. The 3σ upper and lower limits are indicated by arrows. Dashed red lines highlight the regions from which sources were subtracted.
In the text

Fig. 8. Radio power at 150 MHz of low-frequency radio halos as a function of the PSZ2 mass. The lines indicate the extreme cases of the correlations of radio power to cluster mass reported in van Weeren et al. (2021), which were derived using different fitting methods and samples. The dark blue line shows the steepest (orthogonal method) and the light blue line shows the flattest (including candidate halos, Y|X method) correlation. The shaded regions show the 95% confidence interval. The halos reported as ‘van Weeren+21 & lit.’ include a literature radio halo sample (de Gasperin et al. 2015; Cordey 1985; Wilber et al. 2018; Botteon et al. 2019; Bîrzan et al. 2019; Savini et al. 2018, 2019; Bonafede et al. 2018; Hoang et al. 2019a,b; Macario et al. 2013; George et al. 2017; Duchesne et al. 2021).

In the text

	Fig. 9. Flux densities as a function of frequency. Additional GMRT and VLA data at 609 MHz and 1.4 GHz are taken from de Gasperin et al. (2015, 2017).
In the text

	Fig. A.1. Uncertainty maps to Figure 3. The left and center figure show the spectral index uncertainties between 54-144 MHz (left) and 144-323 MHz (center). The right image shows the uncertainty of the curvature map.
In the text

	Fig. A.2. Source-subtracted spectral index uncertainty maps at a spatial resolution of 30″. The two figures are for the frequency ranges 54-144 MHz (left) and 144-323 MHz (right). Dashed red lines highlight the regions from which sources are subtracted.
In the text

	Fig. B.1. Aging model merit function against magnetic field (x-axis) and projected velocity (y-axis) in arbitrary units.
In the text

	Fig. C.1. Results of the halo fitting in LBA (a) and HBA (b). The circles in the left panels indicate the position of the best-fit halo on top of the flux density maps. The center panels show the model, and the right panels present the residual. The green and red regions are excluded from the fitting.
In the text