Galaxy clusters enveloped by vast volumes of relativistic electrons

The central regions of galaxy clusters are permeated by magnetic fields and filled with relativistic electrons1. When clusters merge, the magnetic fields are amplified and relativistic electrons are re-accelerated by turbulence in the intracluster medium2,3. These electrons reach energies of 1–10 GeV and, in the presence of magnetic fields, produce diffuse radio halos4 that typically cover an area of  around 1 Mpc2. Here we report observations of four clusters whose radio halos are embedded in much more extended, diffuse radio emission, filling a volume 30 times larger than that of radio halos. The emissivity in these larger features is about 20 times lower than the emissivity in radio halos. We conclude that relativistic electrons and magnetic fields extend far beyond radio halos, and that the physical conditions in the outer regions of the clusters are quite different from those in the radio halos.

The central regions of galaxy clusters are permeated by magnetic fields and filled with relativistic electrons 1 .When clusters merge, the magnetic fields are amplified and relativistic electrons are re-accelerated by turbulence in the intra cluster medium 2;3 .These electrons reach energies of 1 -10 GeV and, in the presence of magnetic fields, produce diffuse radio halos 4 that typically cover an area of ∼ 1 Mpc 2 .Here we report observations of four clusters whose radio halos are embedded in much more extended, diffuse radio emission, filling a volume 30 times larger than that of radio halos.The emissivity in these larger features is about 20 times lower than the emissivity in radio halos.We conclude that relativistic electrons and magnetic fields extend far beyond radio halos, and that the physical conditions in the outer regions of the clusters are quite different from those in the radio halos.
We used the LOw Frequency ARray (LOFAR 5 ) to search for signatures of non-thermal radiation in the outer regions of galaxy clusters.We examined the radio emission at 144 MHz from 310 massive clusters from the Planck Sunyaev-Zel'dovich catalogue 6 in the LOFAR Two Meter Sky Survey (LoTSS 7;8 ).After a careful removal of contaminating emission from other astronomical sources (see Methods section) we found that in four cases the radio halo emission is embedded in a much larger emission that extends to cluster-centric distances of up to 2-3 Mpc and fills the volume of the cluster, at least up to R 500 (Fig. 1) that is the radius within which the mean mass over-density of the cluster is 500 times the cosmic critical density at the cluster redshift (z).The mass enclosed in a sphere with radius R 500 is M 500 .These clusters are ZwCl 0634.1+4750(z = 0.17, M 500 = 6.65 × 10 14 M ), Abell 665 (z = 0.18, M 500 = 8.86 × 10 14 M ), Abell 697 (z = 0.28, M 500 = 10.99 × 10 14 M ) and Abell 2218 (z = 0.17, M 500 = 5.58 × 10 14 M ).All these clusters are dynamically disturbed 9 and were known to host radio halos 10;11;12 , yet the larger-scale emission discovered with LOFAR, that we name "megahalo", allows us to probe the magnetised plasma in a volume that is almost 30 times larger than the volume occupied by the radio halos.
The radial surface brightness profile of ZwCl 0634.1+4750(see Fig. 2) clearly demonstrates the difference to the radio halo emission.The profiles of the other three clusters are shown in Extended Data Fig. 3.All the profiles show two components: a bright region dominated by the radio halo whose brightness decreases relatively fast with cluster-centric distance and an extended, low-surface brightness component.The region of the profile dominated by the radio halo can be fitted with an exponential function as commonly found for these type of sources 13 .The emission beyond the radio halo shows a shallower profile implying that at 600-800 kpc from the centre a transition occurs.The surface brightness of the large-scale emission is a factor of at least 10 lower than the surface brightness of the radio halo with an average emissivity ∼ 20−25 times lower than the emissivity of radio halos (see Methods section).The low surface brightness, combined with its large size, is the reason why this emission has eluded all previous searches but could be detected by LOFAR.
For two clusters, ZwCl 0634.1+4750 and Abell 665, we present deep observations at even lower frequencies (53 MHz and 44 MHz, respectively) where low-energy relativistic electrons shine brightly.In the observation of ZwCl 0634.1+4750only the brightest part of the largescale emission beyond the radio halo is detected, whereas in Abell 665 almost all of it is visible (Extended Data Fig. 1).The combination of the low (144 MHz) and ultra-low (∼50 MHz) frequency observations allows us to constrain the energetics of the particles responsible for the synchrotron emission via the radio spectral index α (defined as S(ν) ∝ ν α , with S being the flux density and ν the frequency).We obtained α = −1.62 ± 0.22 and α = −1.64 ± 0.13 for ZwCl 0634.1+4750 and Abell 665 respectively.Although the uncertianties on the spectral index measurements are relatively large, in both cases we found evidence that the spectrum is steeper than the spectrum of the central radio halos, which is α = −1.25 ± 0.15 for ZwCl 0634.1+4750 and α = −1.39 ± 0.12 for Abell 665.This adds to the evidence that the newly discovered megahalos are a new phenomenon, distinct from radio halos.
Our discovery confirms that magnetic fields and relativistic electrons fill a much larger volume than ever observed before, therefore requiring ubiquitous mechanisms for the energisation of particles on large scales.The existence of megahalos demonstrates that beyond the edge of radio halos mechanisms operate that maintain a sea of relativistic electrons at energies high enough to emit at frequencies of ∼ 100 MHz.
The surface brightness of this new feature stays fairly constant over more than 500 kpc while the underlying intra cluster medium (ICM) density decreases by a factor of ∼ 5 [14].This can be used to infer the relative contribution to the cluster energy content from thermal and non-thermal components, which has implications, for example, on the Sunyaev-Zeldovoch effect in cosmology 15 .If the magnetic field scales with the ICM density as found in Bonafede et al. (2010) 16 (B 2 ∝ n ICM , where n ICM is the density of the ICM), the ratio between the energy density of non-thermal electrons and the thermal gas energy must increase going towards the outer regions of the cluster.For example, for the cluster Abell 2218, whose ICM density profile has been studied in Pratt et al. ( 2005) 17 , we found that a constant energy density of relativistic electrons with radius can reproduce the observed surface brightness profile.In this case, the ratio between the energy density of non-thermal electrons and the thermal gas energy must increase by a factor of ∼ 3 going from ∼ 0.5 × R 500 to ∼ R 500 .Alternatively, the magnetic field strength must approximately be increased by a factor √ 3 over those distances in order to produce the same radio emission.Reproducing the observed trend of cosmic rays and magnetic fields in such peripheral regions of galaxy clusters, where a mixture of accretion modalities and (re)acceleration mechanisms is present, represents an important challenge for future theoretical models of galaxy clusters.
The steep spectrum that we observe in ZwCl 0634.1+4750 and Abell 665 suggests that turbulence might be responsible for maintaining the relativistic electrons inside a volume of the order of 10 Mpc 3 [18; 19; 20] and that we may be probing a turbulent component different from that powering radio halos.Numerical simulations seem to support this scenario and show that, in addition to the central merger-driven turbulence responsible for radio halos, there is a broader turbulent component, likely related to the accretion of matter onto the cluster, that can accelerate particles 21;22;23 (see Methods section).The observed characteristics of megahalos suggest a change either in the macrophysical or in the microphysical properties of the plasma when moving from the radio halo to the outer region.In the former case, the properties of turbulence may change going towards the periphery, as suggested by simulations.In the latter case microphysical properties such as the acceleration efficiency, the mean free path or transport properties, all of which are related, may change in the outer regions.
Although the mechanisms responsible for the formation of the large-scale emission are still unknown, it is reasonable to assume that the mass of clusters plays an important role in determining the energy budget available for particle acceleration, similar to what happens for radio halos.In fact, more powerful radio halos are hosted in more massive clusters 24 .To understand why we have detected this emission only and exactly in these four clusters, in Fig. 3 we show the mass-redshift distribution of the clusters inspected for this work.The three solid lines show the expected mass-redshift relation taking into account the cosmological surface brightness dimming (SB ∝ (1 + z) −4 ) and assuming a power-law dependence of the large-scale emission surface brightness with the mass of the clusters (SB ∝ M β ; where we assumed three possibilities for β).Since the large-scale emission in ZwCl 0634.1+4750 is detected at ∼ 3σ we impose that the lines go through that point in the diagram to set the normalisation.This means that with current LOFAR observations we expect to be able to detect at more than 3σ the large-scale emission that embeds radio halos in clusters that lie above (i.e. higher mass) and to the left (i.e.lower z) of the assumed dependencies.Interestingly, all the sources that we detected lie in this region.The other clusters with similar masses and redshift either have low-quality data (due to observations taken during disturbed ionospheric conditions that distort the radio signal at low frequencies) or they are particularly complex fields where an accurate subtraction of contaminating sources is not reliable 25;26 .The fact that the clusters where we discovered the new large scale emission lie close or above our estimated detection limit in mass-redshift space, suggests that we may be seeing just the tip of the iceberg of a phenomenon common to a large number of clusters that can come to light in deeper low-frequency radio observations.
In this paper, we have highlighted the differences between the diffuse emission in the central and outer regions of galaxy clusters.The results of this work may shed light on recent findings on a few clusters where it has been shown that radio halos become larger at lower frequencies, reaching largest linear sizes of the order of 2 Mpc 26;27;28 .Interestingly, all these clusters are massive and would lie on the top left region of the diagram in Fig. 3, supporting the idea of a new, emerging population of sources.Finally, evidence of a two component nature of diffuse sources has been claimed in a few non-merging clusters 29;30;31 .However, they do not probe regions of galaxy clusters different from those of radio halos.
Currently, we can only observe megahalos in clusters that meet a certain combination of mass and redshift.However, our study suggests that deeper observations, such as those that will be made with the upgraded LOFAR 2.0 and SKA 32 , will unveil many more clusters showing diffuse emission on such large scales (Fig. 3), thus opening the possibility to systematically explore the peripheral regions of galaxy clusters.Whether megahalos constitute a new class of sources sitting below or embedding radio halos, remains to be seen following deeper searches for this emission on larger samples.Still, the brightness and spectral index profiles suggest that a new type of phenomenology is at play when going to larger distances from the cluster centre.
This discovery shows that relativistic electrons and magnetic fields fill larges swathes of the cosmos.It helps us understand how energy is dissipated during the formation of large-scale structure as well as how particles are accelerated in low-density plasmas.Stars mark clusters where we found the megahalos.Solid lines define the region of the plot where we expect to be able to detect emission beyond radio halos with current LOFAR observations.Possible sources lying in the grey shaded regions would be resolved out (left) or unresolved (right).The dashed line marks the regions where we expect to be able to detect megahalos with LOFAR 2.0 ultra-low frequency observations (the 1σ rms noise for sources with spectral index -1.6 is expected to decrease by a factor of two with respect to current observations).Error bars represent the 1σ uncertainty on the mass estimate 6 .

Methods
Observations and data reduction The properties of the final images presented in this paper are listed in Extended Data Table 1, together with the performed source subtraction (see below) and the figure where they are presented.LOFAR HBA data reduction The LOFAR 144 MHz data presented in this paper are part of the LOFAR Two Meter Sky Survey (LoTSS 33;7;8 ) that is an on-going 120 − 168 MHz survey of the entire northern hemisphere performed with LOFAR high-band antennas (HBA).Data have been processed with the Surveys Key Science Project reduction pipeline v2.2 7;34 , which includes both corrections for direction-independent and direction-dependent effects (prefactor 35;36;37 ; killMS 38;39 ; DDFacet 40 ).We subtract all the sources outside a region of ∼ 0.5 deg 2 containing the galaxy clusters from the uv-data by using the model produced by the pipeline.
We then phase shift the dataset to the centre of the extracted region, we correct for the LOFAR station beam in that direction and we perform additional loops of phase and amplitude calibrations on the target field in order to improve the quality of the final image.This extraction procedure is presented in van Weeren et al. (2021) 41 and is routinely used in LOFAR 144 MHz observations.The data reduction of the clusters of the Planck sample is discussed in detail in Botteon et al. (2022) 42 .Here we performed a more accurate procedure for the source subtrac-tion aimed at removing the contribution of the extended radio galaxies embedded in the diffuse emission.The final images are produced with the multi-frequency deconvolution scheme in WSClean 43;44 .Different resolutions are obtained using different weighting schemes.We did not attempt an in-band spectral analysis because of the narrow frequency range of LoTTS, combined with the uncertainties in the flux density scale and the high rms noise of the in-band images, which especially affects the resolved low signal-to-noise sources, such as megahalos 8 .LOFAR LBA data reduction ZwCl 0634.1+4750 and Abell 665 have been observed with the LOFAR Low-Band Antennas (LBA) for 8 hours in the frequency range 30 -77 MHz and 30 -60 MHz respectively.We used the LBA_OUTER antenna configuration, where only the outermost antennas are selected, because this simplifies the calibration by reducing the primary beam size as well as the electromagnetic crosstalk between the dipoles.Observations were performed in multi-beam mode, with one beam continuously pointing at the calibrator and one beam continuously pointing at the target.The data reduction of the calibrator follows the procedure described in de Gasperin et al. (2019) 37 and it is used to isolate direction-independent systematic effects such as the bandpass, the stations' clocks drifts, and the polarisation misalignment caused by time delays between the X and Y polarisation signals.The solutions are then applied to the target field along with the primary beam correction in the direction of the target.
The self-calibration steps of the target field are described in de Gasperin et al. (2020) 45 .The self-calibration starts with a model obtained from the combination of existing surveys at higher frequencies: TGSS 46 , NVSS 47 , WENSS 48 , and VLSS 49 .We estimate the spectral indices of the sources present in these surveys and we extrapolate their flux densities to LBA frequencies.Then, we estimate direction-independent (field-averaged) differential total electron content (dTEC) solutions by calibrating against the predicted model.Next we solve and correct for the average differential Faraday rotation and second-order beam effects.Sources outside the main lobe of the primary beam are imaged and subtracted from the uv-data before proceeding with a second self-calibration cycle.The main errors that still affect the data at this point are the direction-dependent errors caused by the ionosphere.The procedure to correct for direction-dependent errors is discussed in de Gasperin et al. (2020) 45 and Edler et al. ( 2021) 50 .As a first step, sources in the direction-independent calibrated image are grouped by proximity and the brightest groups are identified to be used as calibrators in the direction-dependent calibration.All the sources in the field are subtracted using the direction-independent model image.We then iterate on the calibrators, starting with the brightest one.The calibrator's visibilities are added back to the data and the dataset is phase-shifted to the calibrator direction.We run several rounds of self calibration and the improved model of the calibrator is re-subtracted to produce a cleaner empty dataset, before repeating the procedure for the next brightest calibrator.
After the wide-field direction-dependent calibration, further improvement of the image quality can be achieved by extraction and self-calibration of a small region around the target of interest.We follow the idea of van Weeren et al. ( 2021) 41 and employ an LBA-specific implementation of the extraction strategy.We use the final model and the solutions of the directiondependent calibration to accurately subtract all sources in the field except for those in a circular region around the targets.The radius of these regions is 23 and 14 for ZwCl 0634.1+4750 and Abell 665 respectively and it is chosen to include sufficient compact source flux density to obtain robust calibration solutions.Then, the data is phased-shifted to the centre of this region and further averaged in time and frequency.After correcting for the primary beam in this direction, we perform several rounds of self-calibration on the target, starting from the model obtained using the solutions of the closest direction-dependent calibrator.The images at high and low resolution are shown in Extended Data Fig. 1.
High-resolution 144 MHz images of the clusters and source subtraction procedure Diffuse emission on scales larger than the radio halos and not associated with radio galaxies is already clearly visible for the four clusters at intermediate resolution (20 − 30 ).In order to highlight this emission, we remove the contribution of the discrete sources embedded in the diffuse emission subtracting them out from the uv-data.We proceeded in two steps: as a first step, we produced an image at high resolution (6 − 10 , shown in Extended Data Fig. 2).The model image that we obtained contains the compact sources in the field and part of the extended sources, such as tails of radio galaxies and the central radio halo.However, it does not contain the diffuse emission from regions beyond the radio halo.In order to produce a source-subtracted image where both the radio halo and the new emission are clearly visible, we removed the clean components in the region of the radio halo from the model image before we subtract the rest from the visibilities.These images, where only the high-resolution source subtraction has been performed, are marked in the Extended Data Table 1 with "HR".Then, we imaged these new visibilities at lower resolution (30 and 60 ).The 60 resolution source subtracted images are shown in Extended Data Fig. 3 (b, d, f and h panels).We note that the large-scale emission does not contain clear residuals or artefacts from the subtracted sources.We inspected the model images obtained at 30 resolution and we retained only the clean components associated with the radio halo and what was left of the extended radio galaxies.Then we subtracted them from the uv-data.Finally, we produced an image at very low resolution (2 , Fig. 1) which only contains the residual large-scale emission.The images where also the radio halo has been subtracted are marked in the Extended Data Table 1 with "HR+RH".
Surface brightness radial profiles Following the approach described in Murgia et al. (2009) 13 and Cuciti et al. (2021) 9 , we derived the azimuthally averaged surface brightness radial profile of the radio emission of the four clusters (Extended Data Fig. 3).We used low-resolution (60 ) images after subtracting only the sources visible in the high-resolution image to have a good compromise between sensitivity to the extended emission and resolution to characterise the profile and possibly distinguish between the radio halo and the more extended emission.We note that the whole extension of the diffuse emission is best detected in the 2 resolution images (Fig. 1) that we used to measure the size of the sources.If the image contains some residuals from bright diffuse sources, we masked them and we excluded the masked pixels when calculating the surface brightness.This is the case in Abell 665, where we masked the residuals of a diffuse patch of emission that does not appear to be associated with the megahalo.Also in Abell 697 we masked the residuals of a bright compact source in the north and a radio galaxy in the south (these sources are visible in Extended Data Fig. 2).For Abell 665 we considered only the southern part of the cluster because the northern part has been crossed by a shock front 51 , which may have altered the properties of the diffuse emission.We averaged the radio brightness in concentric annuli, centred on the peak of the radio halo and chose the width of the annuli to be half of the FWHM of the beam of the image.We considered only annuli with an average surface brightness profile higher than three times the uncertainty associated with the annuli surface brightness.In the images we show the detection limits for each annulus calculated as rms× √ N beam , where N beam is the number of beams in the annulus.All the profiles show a discontinuity.The central annuli before the discontinuity follow an exponential profile, similar to other classical radio halos 13;9 .This first component can be fitted with an exponential law in the form: where I 0 is the central surface brightness and r e is the e-folding radius.
To perform the fit, we first generated a model using Eq. 1 with the same size and pixel size of the radio image and we convolved it with a Gaussian with a FWHM equal to the beam of the image.Then, we azimuthally averaged the exponential model with the same set of annuli used for the radio halo.The resulting surface brightness profile is the fit that takes into account the resolution of the image and the uncertainties associated with the sampling of the radial profile.
The discontinuity in surface brightness is less pronounced in Abell 665 than in the other three cases.Still, there is a second component that is not consistent with the radio halo profile.In addition, the large-scale diffuse emission in Abell 665 shows clear differences to the central radio halo, also in terms of spectrum and emissivity (see below).Hence, we consider it a megahalo.
In two cases, ZwCl 0634.1+4750 and A697, the diffuse emission beyond the classical radio halos is not symmetric with respect to the centre of the radio halo.Therefore, we performed the fit also in the semi-annuli on the left hand side of the red line shown in Extended Data Fig. 3.The discontinuity is present also in this case.

Spectral index analysis
We measured the spectral index of the radio halo in ZwCl 0634.1+4750 in the central region shown in Extended Data Fig. 1 (a).We obtained a flux density of 39.6±4.0 mJy at 144 MHz and 138.3±15.9mJy at 53 MHz, corresponding to α = −1.25±0.15.The uncertainties on the integrated flux densities in this section take into account the systematic error of shock and turbulent acceleration could produce a flatter and less uniform spectrum than expected from turbulence alone 51 .Hence, we focus on the southern part of the cluster, marked by the regions in Extended Data Fig. 1 in order to derive the integrated spectral index.The choice of the regions is based on the surface brightness radial profile (Extended Data Fig. 3, c, d).In particular, the limit between the regions where we measure the flux densities corresponds to the annulus where the surface brightness profile flattens with respect to the classical radio halo exponential function.In these regions we measure a flux density for the radio halo of 120.3 ± 12.1 mJy at 144 MHz and 614.6 ± 62.9 mJy at 44 MHz, corresponding to α = −1.39 ± 0.12.In the region beyond the radio halo, we measure 58.2 ± 6.1 mJy at 144 MHz and 398.± 45.8 mJy at 44 MHz, obtaining a spectral index α = −1.64 ± 0.13.
Emissivity We calculated the volume-averaged emissivity at frequency ν in radio halos and megahalos by assuming that their radio power, P ν comes from a sphere of radius R: Extended Data We estimated the source radii via √ R min × R max 52 , where R min and R max are the minimum and maximum radii of the 3σ contours, respectively, in Fig. 1 (we used the 30 resolution images for radio halos and the 2 resolution images for megahalos).We have subtracted the extended radio galaxies in the fields to the best of our abilities given current techniques.However, we are aware that the central regions of the 2 resolution images shown in Fig. 1 may be affected by residuals from the subtracted radio halos.It seems reasonable to assume that the megahalo also permeates the region of the classical radio halos.Hence, in order to estimate the flux density of the new type of emission, we measured the mean surface brightness in a region that excludes the central halos and then multiply it with the total area within the 3σ contours.A direct measurement of the flux density inside the area delimited by the 3σ contours would give marginal differences, of the order of 5 − 15%.In order to compare these emissivities to the typical emissivity of radio halos 24 , we estimated the flux density at 1.4 GHz, by assuming a conservative spectral index of -1.3, for both radio halos and megahalos.If the spectral index of the newly detected emission is steeper, the emissivity at 1.4 GHz would be even lower.The emissivities of radio halos and megahalos for each cluster are listed in Extended Data Table 2, together with the size and the spectral index (when available) of the sources.The uncertainties include the systematic errors given by the uncertainty on the flux scale, the statistical error associated with the rms noise of the image in the integration area and the subtraction error related to the uncertainty on the subtracted flux.For the latter we used the approach described in Botteon et al. (2022) 42 .We do not take into account the uncertainty in the estimated size of the diffuse emission.The emissivity of megahalos is a factor of ∼ 20 − 25 lower than the emissivity of the radio halos in the same clusters.For comparison, the typical emissivity of classical radio halos at 1.4 GHz ranges between 5 × 10 −43 to 4 × 10 −42 erg s −1 cm −3 Hz −1 in clusters with similar masses 24 .Simulations of turbulence in galaxy clusters High-resolution hydrodynamic, cosmological simulations of the ICM can shed some light on megahalos.We analysed a set of 20 (M 500 ≥ 3 × 10 14 M at z = 0.0) galaxy clusters 22 and used small-scale filtering to calculate the turbulent kinetic energy flux, F turb = ρ σ 3 v /L, within simulated cells (each cell having a 32 3 kpc 3 volume), where ρ is the local gas density and σ v is the dispersion of the velocity field measured within the scale length L 53 .Then we measured the average distribution of F turb in relaxed and disturbed clusters, and in the inner as well as outer regions, where megahalos are observed (approximately 0.4 R 500 < r ≤ R 500 ).
Extended Data Fig. 5 shows that the turbulent kinetic energy flux, F turb , in the central regions of post-merger clusters is elevated by a factor ∼ 10 compared to relaxed clusters.This is consistent with the idea that radio halos are produced by the dissipation of turbulence following a merger 54;18;19;55;56 .However, in the outer regions F turb is remarkably similar in relaxed and disturbed clusters.This shows the presence of a baseline level of turbulence that is induced by the continuous accretion of matter in cluster outskirts.Hence, this level of turbulence is likely to be more common to all clusters.In this picture, megahalos may also be generated via Fermi II re-acceleration in more relaxed clusters without central radio halos.Future deeper observations, such as those that will be made with LOFAR 2.0, will allow us to test this scenario.

Figure 1 :
Figure 1: Radio images of the four galaxy clusters.a) ZwCl 0634.1+4750,b) Abell 665, c) Abell 697, d) Abell 2218.Black contours: LOFAR 144 MHz at 30 resolution, showing the location of the radio halos.Only sources in the high-resolution image have been subtracted.Contours at (4, 8, 16, 32) × σ.Orange image and white contours: LOFAR 144 MHz image at 2 resolution.All sources, including the radio halos, have been subtracted.Contour at 3σ confidence level.The radius of the light blue circle corresponds to R 500 .

Figure 2 :
Figure 2: Large-scale emission from the intracluster medium in ZwCl 0634.1+4750.In blue, the gas density distribution in a simulated massive galaxy cluster is shown 22 only for illustrative purposes.In orange the LOFAR 144 MHz 2 resolution image of the cluster ZwCl 0634.1+4750 is shown after the subtraction of all the sources in the field, including the radio halo.The central region, outlined in black, shows the LOFAR 144 MHz 30 resolution image of the classical radio halo.Inset: Surface brightness radial profile of the diffuse radio emission in ZwCl 0634.1+4750.On the right we show the 60 resolution 144 MHz image and the annuli we used to measure the surface brightness reported in the plot on the left (see Methods section for more details).

Figure 3 :
Figure 3: Mass-redshift diagram for the Planck clusters in LoTSS.Stars mark clusters where we found the megahalos.Solid lines define the region of the plot where we expect to be able to detect emission beyond radio halos with current LOFAR observations.Possible sources lying in the grey shaded regions would be resolved out (left) or unresolved (right).The dashed line marks the regions where we expect to be able to detect megahalos with LOFAR 2.0 ultra-low frequency observations (the 1σ rms noise for sources with spectral index -1.6 is expected to decrease by a factor of two with respect to current observations).Error bars represent the 1σ uncertainty on the mass estimate6 .

1 :
LOFAR LBA images.a) ZwCl 0634.1+4750,b) Abell 665.Orange images and white contours: LOFAR data at 2 resolution after subtraction of all sources, including the radio halo.The contour is at 3σ. Black contours: high-resolution data (18.1 × 18.1 for ZwCl 0634.1+4750 and 23 × 14 for Abell 665), contours at (3, 6, 12, 24) × σ.The magenta and green regions mark the area used to calculate the integrated spectral index of the radio halos and of the newly discovered emission, respectively.The radius of the light blue circle corresponds to R 500 .

Extended Data Figure 3 :
Surface brightness radial profiles.From top to bottom: ZwCl 0634.1+4750,Abell 665, Abell 697, Abell 2218.a, c, e, g: surface brightness radial profile measured on the regions shown on the maps on the right.The blue line is the exponential fit to the radio halo data.The dashed lines mark the 1σ detection limit for each annulus (or semi-annulus).The vertical line marks R 500 .Shaded regions represent the 1σ uncertainty.For ZwCl 0634.1+4750 and Abell 697 we derived the profiles in the entire annuli (blue shaded region) and only in the semi-annuli on the left hand side of the red line shown on the right panels (orange shaded region).b, d, f h: 144 MHz images at 60 resolution after subtraction of the sources visible in the high-resolution images.

Extended Data Figure 5 :
Simulated turbulence.Probability distribution function of the turbulent kinetic energy flux in post-merger or relaxed clusters simulated by Vazza et al. (2011)22 , either within the region of classic radio halos (dashed line), or within the outer region (solid line) where the newly discovered megahalos are located.

Table 1 :
Properties of the LOFAR images.
Extended Data

Table 2 :
Properties of radio halos and megahalos.