The Planck clusters in the LOFAR sky IV: LoTSS-DR2: statistics of radio halos and re-acceleration models

Context. Di ﬀ use cluster-scale synchrotron radio emission is discovered in an increasing number of galaxy clusters in the form of radio halos, probing the presence of relativistic electrons and magnetic ﬁelds in the intra-cluster medium (ICM).


Introduction
Galaxy clusters are the largest gravitationally bound structures in the Universe.They contain the mass of ∼10 14 −10 15 M in regions with a size of 2−3 Mpc.They consist in particular of dark matter (70−80% of the cluster mass), a hot and low-density gas (for a ∼15−20%), the intracluster medium (ICM), and a few percent of cluster galaxies (e.g.Kravtsov & Borgani 2012, for a review).Observations of galaxy clusters in the radio band show the presence of diffuse synchrotron radiation in merging galaxy clusters in the form of radio haloes (RHs) at the cluster centres and radio relics at the cluster outskirts.Mini-haloes (MH) are also observed as less extended objects surrounding the central radio galaxy, which is associated with the brightest cluster galaxy typically in relaxed clusters (e.g.Feretti et al. 2012;van Weeren et al. 2019, for reviews).Recent evidence is the detection of radio bridges connecting pairs of massive clusters (e.g.Govoni et al. 2019;Botteon et al. 2020;Brunetti & Vazza 2020).All these sources prove the existence of cosmic-ray electrons and magnetic fields in the ICM and pose fundamental questions about the origin of these components, their impact on the thermal ICM (microphysics), and their connection with the cluster dynamics and evolution (e.g.Brunetti & Jones 2014, for a review).
RHs are diffuse, megaparsec-sized, synchrotron radio sources with steep radio spectra (α > 1, with f (ν) ∝ ν −α ) that are observed in the central regions of a fraction of galaxy clusters.A statistical connection between RHs and the cluster dynamical status has been found: clusters with RHs often show signs of merger activity from the analysis of X-ray observations (e.g.Buote 2001;Cassano et al. 2010b;Wen & Han 2013;Cuciti et al. 2015;Giacintucci et al. 2017).These observations provided support to the hypothesis that the turbulence that is generated during cluster mergers re-accelerates pre-existing fossil and/or secondary electrons in the ICM to the energies that are required to produce the observed radio emission (e.g.Brunetti & Jones 2014).
According to turbulent re-acceleration models (e.g.Brunetti et al. 2001;Petrosian 2001;Fujita et al. 2003;Cassano & Brunetti 2005;Brunetti & Lazarian 2007, 2011, 2016;Beresnyak et al. 2013;Donnert et al. 2013;Miniati 2015;Pinzke et al. 2017;Nishiwaki & Asano 2022), the formation history of RHs depends on the cluster merging rate throughout cosmic epochs and on the mass of the hosting clusters themselves, which ultimately sets the energy budget that is available for the acceleration of relativistic particles.A statistical model based on semi-analytic calculation of the galaxy clusters formation history (Press & Schechter 1974;Lacey & Cole 1993) and on a simple prescription (i.e.homogeneous conditions) to estimate the turbulence injected during mergers and the synchrotron spectra that are generated by the turbulent (re)acceleration process has been developed and tested in the past decade (e.g.Cassano & Brunetti 2005;Cassano et al. 2006Cassano et al. , 2010aCassano et al. , 2012)).A key expectation of this scenario is that the statistical properties of the population of RHs should depend on the frequency of the observations, since the synchrotron spectra of the haloes are characterised by a steepening frequency ν s .Above this frequency, the spectrum of RHs gradually steepens, and thus the halo could be difficult to observe at ν > ν s .Since the value of ν s is connected to the energetics of the merger events that generate the haloes, observations at ∼1 GHz frequency should discover RHs preferentially in massive objects undergoing energetic merging events.RHs in less massive merging systems should be difficult to observe at frequency ≥1 GHz (Cassano et al. 2006(Cassano et al. , 2010a)).Statistical studies of clusters, such as the GMRT Radio Halo Survey (Venturi et al. 2007(Venturi et al. , 2008;;Kale et al. 2013Kale et al. , 2015)), and its extension to mass-selected clusters (from the Planck SZ cluster catalogue; Planck Collaboration XXIX 2014) led to the first statistical evidence that RHs are predominately found in merging clusters, whereas clusters without diffuse emission are typically relaxed (e.g.Brunetti et al. 2007;Cassano et al. 2010bCassano et al. , 2013)).They have also shown that RHs become progressively less common in less massive clusters (e.g.Cuciti et al. 2021a,b), thus indicating an observational connection between RH formation and the energetics of the hosting systems.
However, the key expectation of the model that a large fraction of RH has very steep spectra (α > 1.5 at ∼1 GHz frequency) and glows preferentially at lower frequencies remains poorly explored so far.A number of haloes with very steep spectra were detected at lower frequencies (e.g.Brunetti et al. 2008;Macario et al. 2013;Wilber et al. 2018;Duchesne et al. 2021;Bruno et al. 2021;Rajpurohit et al. 2021), thus supporting the idea that such a population might exist.Evidence of clusters also emerges that exhibit some level of dynamical disturbance hosting multi-component haloes or a central MH, such as emission surrounded by a diffuse radio emission on a larger scale (e.g.Savini et al. 2018Savini et al. , 2019;;Biava et al. 2021;Riseley et al. 2022).In this respect, LOFAR has recently enabled observations of galaxy clusters at frequencies <200 MHz with unprecedented high-sensitivity and resolution.More specifically, LOFAR is carrying out wide and deep surveys of the entire northern sky at 120−168 MHz and 42−66 MHz in the context of the LOFAR Two-meter Sky Survey (LoTSS; Shimwell et al. 2017) and LOFAR LBA Sky Survey (LoLSS; de Gasperin et al. 2021), respectively.One of the main goals of these surveys is the discovery of new diffuse megaparsec-scale radio sources in galaxy clusters, providing samples that are suitable for testing the formation models.A first step in this direction has been carried out by van Weeren et al. (2021), who performed a first statistical investigation of diffuse emission in galaxy clusters selected from the second Planck catalogue of SZ sources (PSZ2; Planck Collaboration XXVII 2016) that have been covered by the first LoTSS Data Release (LoTSS-DR1; Shimwell et al. 2019).A further important step to constrain the spectrum of RH has been carried out by Di Gennaro et al. (2021a,b), who observed a small sample of massive clusters at high redshift (z ≥ 0.6) with LOFAR and then followed up the detected haloes with the uGMRT.In line with the models, about 50% of these RHs exhibit a very steep spectral index (i.e.α ≥ 1.5 between 150−650 MHz) and are found among the less massive clusters in that sample.
More recently, we have started a large project1 to study diffuse radio emission in the ICM of the galaxy clusters selected from the PSZ2 clusters that have been covered by the second LoTSS Data Release (LoTSS-DR2; Shimwell et al. 2022), covering 5634 square degrees (27% of the northern sky).In Botteon et al. (2022, Paper I), we presented the sample, described the methods and data we used, classified the cluster radio sources, and provided measurements of different quantities.In Bruno et al. (2023) we presented the procedures and derived the upper limits to the radio flux density and power of clusters in which we do not observe diffuse emission.In Zhang et al. (2023) we derived the X-ray properties for clusters in the sample with available Chandra, and/or XMM-Newton, archival data.In Cuciti et al. (in prep.), we study the scaling relations of RHs in the sample by comparing the distribution of RHs and upper limits in the cluster mass versus radio power diagram.In Jones et al. (2023) we present the statistical analysis of radio relics in the same sample.In this paper, we focus on RHs and investigate their flux density and redshift distribution, their occurrence as a function of the cluster mass and redshift, and their connection with the cluster dynamics, and we compare these findings with the expectations from theoretical models.

Cluster sample
The PSZ2 catalogue (Planck Collaboration XXVII 2016) contains 1653 SZ-sources detected over the entire sky.Three hundred and nine of these sources lie in the LoTSS-DR2 footprint and are the subject of Paper I, and 281 out of the 309 sources have mass and redshift information.In Paper I we classified the diffuse radio emission in these clusters by visually inspecting a set of LOFAR images at different resolutions (with and without individual point source subtraction) together with the optical/X-ray overlay images.To make this classification most objective and easily reproducible, we made use of a decision tree to classify the diffuse emission.Briefly (see Paper I for details), the radio sources were classified as RHs, which are extended sources that occupy the region in which the bulk of the X-ray emission from the ICM is detected, radio relics (RR), which are elongated sources whose position is offset from the bulk of the X-ray emission from the ICM, candidate radio haloes/relics (cRH/cRR), in which rather clear RH/RR emission is present, but the absence of Chandra or XMM-Newton X-ray observations prevents them from being firmly claimed, uncertain (U), which signifies that the emission was either significantly affected by calibration or subtraction artefacts or that it did not fall easily in the categories of RHs and RRs, no diffuse emission (NDE), which applies to objects that lack diffuse emission that is not associated with an AGN, and not applicable (N/A), which means that the emission cannot be adequately classified because of poor data quality.It is important to note that in these studies, we preferred to consider MHs together with RHs, we therefore refer generically to both as RHs (see Paper I).
As described in Paper I, we used the Halo-Flux Density CAlculator2 (Halo-FDCA; Boxelaar et al. 2021) to measure the integrated flux density from the observed RHs assuming exponential profiles for the fitting.This model has two main parameters: the central surface brightness (I 0 ), and the e-folding radius (r e ).These were determined for each RH.As suggested by Murgia et al. (2009), when calculating the Halo-FDCA derived flux densities, we integrated the best-fit models up to a radius of three times the e-folding radius.This choice led to a flux density that was ∼80% of the density that would be obtained by integrating the model up to infinity.It is motivated by the fact that haloes do not extend indefinitely.Because the integrated flux density measurements were obtained at 144 MHz, the k-corrected radio powers at 150 MHz were derived accordingly to the usual formula and assuming α = 1.3 (see Eq. ( 5) in Paper I); typical values are indeed in the range α = 1−1.5, and thus our measurements are only marginally affected by the adopted (unknown) radio spectral index (see e.g.Feretti et al. 2012;van Weeren et al. 2019).For galaxy clusters without detected diffuse radio emission, we derived upper limits to the radio power of a possible halo by injecting the visibilities of simulated RHs in the LOFAR data.The injection technique and the resulting upper limits to the flux density are reported in Bruno et al. (2023).
The radio power-mass correlation for cluster RHs and their comparison with the upper limits are discussed in Cuciti et al. (in prep.).Here we used the P 150 MHz −M 500 best-fit relation in the form log P 150 MHz 10 24.5 W Hz −1 = B log M 500 10 14.9 M + A, with A = 1.1 ± 0.1 and B = 3.59 ± 0.48, which has a measured scatter σ raw ∼ 0.4 (see Cuciti et al., in prep.).
For the statistical analysis of the RH properties, we focused on cluster sub-samples spanning a redshift range of 0.07 < z < 0.5.The distribution in redshift and mass of the clusters we studied is shown in Fig. 1.The lower redshift cut was made to exclude the very nearby Universe, for which the cluster statistics is low because of the limited volume.The high-z cut was imposed because the Planck selection is such that only the most massive clusters are detected, which are rare at these redshifts, especially in the relatively small DR2 area (see also the discussion of the incompleteness of the redshift information in Planck Collaboration XXVII 2016).To select our sample, we also considered the Planck completeness.Specifically, we converted the selection function (available in the Planck archive for the full survey region and for a signal-to-noise threshold of 4.5), originally defined in the SZ signal -cluster size plane, into the M−z plane as described in Planck Collaboration XXVII (2016), and determined the boundary line above which the detection probability is higher than a given percentage.As reference, we used the 50% Planck completeness line in (M 500 , z), (see Fig. 1).This means that we may have missed 50% of the clusters whose mass is close to the 50% completeness line.However, as soon as the mass increases, the completeness increases as well, so that we have, for instance, an 80% completeness for clusters with M 500 5 × 10 14 M at z ∼ 0.3.We do not expect this choice to introduce significant biases in our results because there are no significant differences between the completeness functions for regular and disturbed clusters, as shown by Planck Collaboration XXVII (2016).
We refer with RH to both RH and cRH, unless indicated differently.This also includes cRH* and RH*, for which sources the radio power could not be measured (see Paper I for details).They were therefore not used to estimate the RH flux density distribution.We also point out again that we did not distinguish between MHs and RHs (see Paper I), which means that some MH might be included in the RH category.The U cases were treated with caution because they might host an RH or different diffuse sources that are not necessarily associated with the ICM.When RR clusters with U emission at their centre are considered as well, the total number of U cases in the sample is 30.

Detecting radio haloes in LoTSS DR2
In Fig. 2 we show either the radio power (red points) or upper limits (green triangles) of the haloes as a function of redshift.
A43, page 3 of 12 This plot allows determining the minimum power of a RH that could be detected in LoTSS-DR2.We also report an analytic expression that we used in previous papers (Cassano et al. 2010a(Cassano et al. , 2012) ) to estimate the minimum flux density of an RH that can be detected in a given survey by assuming that the halo is detectable when the integrated flux within 2×θ e (θ e is the angular size corresponding to the e-folding radius r e ) gives a signal-tonoise ratio ξ, i.
where F rms is the rms noise in µJy, and θ b is the beam angular size in arcsec.The corresponding minimum radio power P min (z) is reported in Fig. 2 as the black line assuming F rms = 200 µJy beam −1 , θ b = θ b (z) depending on redshift with a fixed linear size of 150 kpc (see the data reduction strategy described in Paper I) and θ e corresponding to r e = 170 kpc (which is about the median values of r e in our sample).With this parameter choice, Eq. ( 2) with ξ = 5 roughly describes the behaviour of the upper limits as a function of redshift.The blue line in Fig. 2 was obtained by applying the P 150 MHz −M 500 best-fit relation to the 50% Planck completeness line reported in Fig. 1.It indicates the minimum power of RHs in PSZ2 clusters under the assumption that they follow the radio power-mass correlation.The fact that this line is always above the line traced by the upper limits indicates that LOFAR would be able to detect RHs in clusters with a mass above the 50% completeness line.As a consequence, to compare model expectations with observations, we used the blue line to determine the minimum power of a detectable RH in PSZ2 clusters that lie above the 50% completeness line at each redshift (see Sect. 5, for details).

Merger-driven turbulent re-acceleration scenario
Our LOFAR sample spans a range of cluster masses and redshifts that was not probed by observations so far (Fig. 1).Furthermore, it is the first large statistical sample that was observed at low frequencies, which makes it ideal for testing models.The main goal of the paper is to determine whether the observed properties of the RH population in this sample are consistent with model expectations.To explore this point, we considered a scenario in which RHs form in galaxy clusters during cluster-cluster mergers due to the turbulent re-acceleration of relativistic electrons.

Basic model
A detailed description of the model that we used can be found in Cassano & Brunetti (2005) and Cassano et al. (2006), while applications to RH predictions for future surveys (with LOFAR, Apertif on WSRT, and ASKAP) can be found in Cassano et al. (2010aCassano et al. ( , 2012)).In this section we provide a summary of the theoretical framework and of the most important implications for RH statistical properties and connection with the host clusters.
We modelled the properties of the RHs and their cosmic evolution by means of a Monte Carlo approach, which is based on the semi-analytic model of Lacey & Cole (1993, i.e., the extended Press & Schechter 1974) to describe the hierarchical process of formation of galaxy cluster dark matter haloes.The merger history of a synthetic population of galaxy clusters is followed back in time, and the generation of turbulence in the ICM is estimated for each merger identified in the merger trees.In these calculations, turbulence is assumed to be injected in the cluster volume that is swept by the sub-clusters, which is bound by the effect of the ram-pressure stripping.The turbulent energy is calculated as a fraction η t (∼0.1−0.3) of the PdV work done by the sub-clusters falling into the main cluster.In these models, the turbulent energy, acceleration rate, and magnetic field per unit volume are considered constant (i.e.homogeneous models, Cassano et al. 2010a).
The most important expectation of turbulent re-acceleration scenarios is that the synchrotron spectrum of RHs gradually becomes steeper above a frequency, ν s , that is determined by the competition between acceleration and energy losses and is connected to the energetics of the merger events that generate the haloes (e.g.Fujita et al. 2003;Cassano & Brunetti 2005).At higher frequencies, the synchrotron spectrum of haloes steepens.Following Cassano et al. (2010a), we adopted the convention that RHs have a spectral index α = 1.9 between ν s /2.5 and ν s .In homogeneous models, the frequency ν s depends on the acceleration efficiency χ and on the mean magnetic field strength in the RH volume B as ν s ∝ B χ 2 /( B 2 +B 2 cmb ) 2 , (e.g.Cassano et al. 2006Cassano et al. , 2010a)), where B cmb = 3.2(1 + z) 2 µG is the equivalent magnetic field of the cosmic microwave background (CMB) radiation.Monte Carlo simulations of cluster mergers allow us to evaluate χ from the estimated rate of turbulencegeneration and the physical conditions in the ICM.We can then derive the dependence of ν s on cluster mass, redshift, and merger parameters in a statistical sample of synthetic clusters (see e.g.; Cassano & Brunetti 2005, for details).
The three main model parameters are η t , the typical radius of RH, R H , and the magnetic field in the RH volume B .The dependence of the model expectations on the parameter values is explored in a number of papers (e.g.Cassano et al. 2006Cassano et al. , 2008)).In this paper we adopted a reference set of model parameters, namely a volume-averaged magnetic field strength B = 2 µG (in line with e.g., Bonafede et al. 2010), independent of cluster A43, page 4 of 12 redshift, η t = 0.2, and an RH size R H 400 kpc (which is in line with the median size of the haloes in our sample).This parameter set has been routinely used in recent papers (Cassano et al. 2019;Botteon et al. 2021;Di Gennaro et al. 2021b) and was found to reproduce the observed RH statistics at low (LOFAR) and high (uGMRT) radio frequencies in a sample of high-z clusters (Di Gennaro et al. 2021b).Furthermore, based on our previous works on the RH statistics at 1.4 GHz (e.g.Cassano et al. 2006) and on the comparison between LOFAR and uGMRT RH statistics for high-z clusters (e.g.Di Gennaro et al. 2021b), we expect that the general results of this paper are independent of the adopted parameter values.An exploration of the full range of model parameters will be performed upon completion of the LoTSS.

Occurrence of radio haloes
As a consequence of the adopted scenario, the population of RHs is expected to be made of a complex mixture of sources with different spectra.Massive (and hot) clusters show the tendency of generating haloes with spectra that are flatter than those in less massive systems.
In order to estimate the occurrence of RHs in a survey at a given observing frequency ν o , we assume that only haloes with ν s ≥ ν o can be observed.Energy arguments imply that RHs with ν s ≥ 1 GHz are generated in connection with the most energetic merger events in the Universe.Only these mergers can generate enough turbulence on megaparsec scales and potentially produce the acceleration rate that is necessary to maintain the relativistic electrons emitting at frequencies higher than 1 GHz (Cassano & Brunetti 2005).In general, the fraction of clusters with RHs in this model increases with the cluster mass because more massive clusters are more turbulent (e.g., Vazza et al. 2006;Hallman & Jeltema 2011) and thus are more likely to host an RH.This agrees with the fact that current surveys carried out at ν o ∼ 1 GHz detect RHs only in the most massive and merging clusters (e.g., Cassano et al. 2013;Cuciti et al. 2015Cuciti et al. , 2021b)).For similar energy arguments, RHs with lower values of ν s , i.e. or RHs with ultra-steep radio spectra (USSRHs, defined as haloes with α > 1.9 between 250−600 MHz), must be more common because they can be generated in connection with less energetic phenomena, e.g.major mergers between less massive systems or minor mergers in massive systems, that are more common in the Universe.
In Cassano et al. (2012) we showed that the fraction of clusters with haloes increases at lower values of ν o .The size of this increment depends on the considered mass and redshift of the parent clusters.It is larger for lower cluster masses and at higher redshifts.
In For the same masses, the percentages of RHs with a very steep spectrum increase considering higher redshifts.

Radio halo luminosity function
The luminosity functions of RHs (RHLFs) with ν s ≥ ν 0 (the expected number of haloes per comoving volume and radio power that can be observed at a frequency ν 0 ) can be estimated by where dN H (z)/dM dV is the theoretical mass function of radio haloes with ν s ≥ ν 0 , which is obtained by combining Monte Carlo calculations of the fraction of clusters with RHs and the Press & Schechter (PS) mass function of clusters (e.g.Cassano et al. 2006).We estimated dP(ν 0 )/dM from the radio power-mass correlation, but with respect to previous papers, this was taken directly from the new correlation obtained at 150 MHz for the LOFAR discovered RH (Eq.( 1) and Cuciti et al., in prep.).
A43, page 5 of 12 The RHLFs at three different redshifts are reported in Fig. 3 (right panel).As already discussed in Cassano et al. (2006Cassano et al. ( , 2010a) ) and, the shape of the RHLF flattens at low radio powers because of the expected decrease in efficiency of the particle acceleration in the case of less massive clusters.However, the flattening at low power is less relevant than that expected considering higher values of ν 0 (see Cassano & Brunetti 2005) because RHs with lower values of ν s contribute to the low-power end of the lowfrequency luminosity function.Finally, we note that the normalisation of the RHLFs decreases with increasing redshift (Fig. 3, right panel) due to the evolution with z of both the cluster mass function and the fraction of galaxy clusters with RHs (Fig. 3, left panel; see also Cassano et al. 2006).

Comparison between the observed radio haloes in LoTSS DR2 and model predictions
In this section we compare model expectations and the data from LoTSS DR2.As already mentioned, we used reference values of model parameters ( B = 2 µG, η t = 0.2, R H 400 kpc) that were used in several previous papers (Sect.4).
The Planck selection function (see Fig. 1, Planck Collaboration XXVII 2016) implies that our measures arise from the combined effect of a simultaneous mass and redshift selection.In order to proceed with a comparison between observations and models, model predictions must be calculated including the same selection effects as in the observed sample.
The number of RHs with flux ≥ f min (z) in the redshift interval, ∆z = z 2 −z 1 , can be obtained by integrating the RHLF (Eq.( 3)) above a given f min (z) 3 , The estimate of f min (z) was the more critical aspect in previous papers, where we derived the number of RHs that is expected to be discovered in future radio surveys.The strength of the current approach is that now we know the sensitivity of LOFAR to detect RHs in PSZ2 clusters (see Sect. 2, for details) because we derived upper limits to the radio flux density of NDE clusters (see Sect. 2 and Bruno et al. 2023).In Sect. 2 we showed that the radio power of these upper limits as a function of redshift is well described by Eq. ( 2) and that this power is always lower than the radio power expected for galaxy clusters with M 500 ≥ M 50%,Planck (z) hosting an RH on the P 150 −M 500 correlation (the blue line in Fig. 2).In other words, the Planck selection function (M 500 ≥ M 50%,Planck (z)) combined with the P 150 −M 500 correlation (Eq.( 1), see also Cuciti et al., in prep.) would determine the minimum radio power of a halo on the correlation we can have in our sample for each redshift.
To derive the expected number of LOFAR-detectable RHs in PSZ2 clusters in the LoTSS-DR2 footprint, we used Eq. ( 4) with f min (z), derived as explained above (blue line in Fig. 2).In order to compare the predicted number of RHs with the observed number, we first need to normalise the total number of clusters we have in the model (from the Press & Schechter mass function) so that this matches the number of PSZ2 clusters in our observations.To do this, we divided the predicted number of RHs by the ratio of the number of clusters in the model by the observed number of PSZ2 clusters above the 50% completeness line.
3 For a similar approach to radio relics, see Brüggen & Vazza (2020).The agreement between the observed and expected flux density distribution of the RHs is very good.The deviation of the observed points at the highest flux density could be due to the insufficient statistics (just one observed RH).
The number density of RHs in Fig. 4 is due to the contribution of RHs observed at different redshifts.Both IC losses and the merging rate for a given mass depend on redshift.It is therefore important to compare model predictions with the observed distribution of RH with z.By using the same normalisation procedure as described above, we derived the redshift distributions.In Fig. 5 we show the cumulative number of RHs within a given redshift N H (<z), (left panel) and the number of RHs per redshift bins N H (z, ∆z), (right panel).In the left panel of Fig. 5, we also report the upper boundaries to the total number of RHs within a given redshift bin that were obtained by simply adding the A43, page 6 of 12   1).numbers of RHs and U.In both plots, the observed and expected distribution of RHs agree very well.

Radio haloes in different redshift and mass bins: Mass dependence
As already outlined in Sect.4, the fraction of clusters with RHs is expected to increase with cluster mass in the framework of the particle acceleration scenario, independently of the selected redshift range (see Fig. 3).In this section, we investigate the dependence of the RH occurrence on the cluster mass in the PSZ2 clusters in the LoTSS-DR2 area.This is a difficult task because we need to preserve a high statistics and at the same time consider that the Planck selection (see the dotted black line in Fig. 1) leads to higher masses at higher redshifts; this implies that our data provide a constraint that is derived from a combination of mass and redshift.We divided the cluster sample into three redshift bins.In each redshift bin, the sample was divided into two sub-samples with a similar number of clusters: a high-mass sample, and a low-mass sample (see Fig. 6, left panel).Note that, because of the horizontal cuts in mass selection, some clusters just below the 50% completeness line can be included in the sample to derive the halo fraction in the various (M, z) bins.In Table 1 we report the observed fraction of clusters with RHs, f RH (obs), in the three redshift ranges for the highmass and low-mass subsamples.We show the cases in which   For each bin, the reported observed minimum and maximum value of f RH is f RH (MC) with U = no RH case and f RH (MC) with U = RH case, respectively (see Table 1).
U clusters in the sample are considered as clusters without RH (U = no RH, Table 1, upper panel) and as RH clusters (U = RH, Table 1, lower panel).
To take the effect on the RH fraction due to the uncertainty associated with the statistical error on the masses into account, we ran a Monte Carlo routine.We randomly extracted the mass of each cluster from a Gaussian distribution with a median value µ = M 500 and standard deviation σ = σ M 500 , where σ M 500 is the value of the uncertainty on the mass, as reported in the Planck catalogue (Planck Collaboration XXVII 2016).Then we split the clusters into two sub-samples according to their mass and calculated the fraction of RHs in each sub-sample.The derived fraction (i.e. the mean of the distribution) and its uncertainty (i.e. the standard deviation) are reported in Table 1 as f RH (MC).In the right panel of Fig. 6, we show the fraction of clusters with RHs, f RH , derived in the high-mass (red shadowed region) and low-mass (black shadowed region) samples in the three redshift bins.For each mass and redshift bin, the minimum and maximum reported values are f RH (MC) with U = no RH case and f RH (MC) with U = RH case, respectively (see Table 1).In all three considered redshift ranges, f RH clearly increases from the low-to the high-mass sample.
We used the same model (as described in Sect.4) and predicted the fraction of clusters with RH in the same mass and redshift bins of the observations.We report in Fig. 7 the comparison between the expected and observed values (the same as reported in Fig. 6) in the low (left panel) and in the high (right panel) mass bins.The expected and derived RH fractions agree fairly well.The model roughly under-reproduces the RH occurrence at low redshift in the high-mass subsample.However, this might be explained by the lack of a representative number of massive clusters at low redshift (due to the limited volume of the Universe).

Quest for haloes with a very steep spectrum
Overall, we have shown that a simple version of the reacceleration models that is based on homogeneous conditions in the ICM and Monte Carlo simulations of merger-turbulent connection, and that uses reference parameters already adopted in the past, provides an excellent description of the LOFAR observations (Sect.5).
Model predictions are anchored on the assumption that RH can be observed at ν o only if ν s ≥ ν o .In particular, we find that the model predicts that most haloes (70−80%) detected in PSZ2 clusters in LoTSS-DR2 have a steep spectrum (ν s ∼ 150−600 MHz); this is shown in Fig. 5 (right panel).Unfortunately, these predictions cannot be tested with the current data as A43, page 8 of 12 Based on current data, we can compare the fraction of clusters with RHs observed in our LOFAR sample with that observed with GMRT at 610 MHz in PSZ2 clusters with M 500 ≥ 6 × 10 14 M and z ≤ 0.35 (Cuciti et al. 2021b).In this mass and redshift range, we measure a fraction of clusters with RHs ∼70% in our LOFAR sample, compared to ∼45% that is measured in the GMRT sample (see Table 2 for details).
We derived the expected fraction of clusters with RHs in the cluster population with M 500 ≥ 6 × 10 14 M in the redshift range z 0.08−0.35by assuming our statistical model (see Sect. 4) and requiring ν s > 610 MHz and ν s > 150 MHz, to be compared with the GMRT and LOFAR fractions, respectively.We report these fractions in Table 2, together with those derived from the observations.This can be considered indirect proof of the presence of USSRH in the PSZ2-DR2 sample.The observed fractions of clusters agree well with RHs and those derived from the model.The model predicts that the difference between the low-and high-frequency fractions is caused by the intervening population of very steep spectrum haloes (∼33% of haloes at these cluster masses and redshifts) that become visible preferentially at the lower frequencies.

Radio haloes and connection with cluster dynamics
Quantitative measurements of the morphology of the X-ray emission have proved to be an effective way to characterise the dynamical state of large samples of galaxy clusters (e.g.Buote 2001;Santos et al. 2008;Cassano et al. 2010b;Rasia et al. 2013;Rossetti et al. 2017;Lovisari et al. 2017, and references therein).
The combination of two morphological parameters as the concentration parameter, c sensitive to the core robustness, and the centroid shift, w, sensitive to the presence of substructures, is an optimal choice to characterise the merger status of galaxy clusters in relation to their diffuse cluster-scale emission (Cassano et al. 2010b;Cuciti et al. 2015Cuciti et al. , 2021a)).In Zhang et al. (2023) we derived the morphological parameters over a physical scale of R ap = 500 kpc centred on the X-ray emission peak.The measurements are reported in Paper I.Here we briefly review the definition of the parameters and then show the results in relation to the analysis of the RH versus NDE clusters.
The concentration parameter was introduced by Santos et al. (2008) as the ratio of the photon flux within two circular apertures to effectively identify cool cores even at high redshift.
Here we adopt the choice of apertures made by Cassano et al. (2010b), where F(r < 100 kpc) and F(r < R ap ) are the exposure-corrected counts within the apertures of 100 kpc and 500 kpc, respectively.The centroid shift (Mohr et al. 1993;Poole et al. 2006) is defined as the variance of the separation between the X-ray peak and the centroid of the emission obtained within a number of apertures of increasing radius out to R ap , where ∆ i is the distance between the X-ray peak and the centroid of the ith aperture.It traces the variation in the position of the centroid introduced by the substructures in the X-ray emission.The detailed X-ray data reduction and analysis processes are described in Paper I and Zhang et al. (2023).
In the left panel of Fig. 8 we report the mass-z distribution of the clusters above the 50% Planck completeness line with information about their dynamical status (coloured dots) together with the distribution of clusters in the LoTSS-DR2/PSZ2 sample without X-ray observations (grey points)4 .Of the 164 clusters within z ∼ 0.07−0.5 and above the 50% line, 93 have observations in the XMM-Newton and/or Chandra archives; here we refer to these clusters as to the morphological sample.These data were used to derive the morphological parameters c and w (Paper I, Zhang et al. 2023).When XMM-Newton and Chandra measurements were both available, we used a value obtained from the combination of the two (see Paper I). Limitations on the measure of c and w parameters are discussed in Zhang et al. (2023).We conclude that the values of c and w are fairly unbiased since our clusters are below z = 0.5 and their number of X-ray counts is greater than 1000.In the right panel of Fig. 8, we report the c−w distribution of these 93 clusters together with the reference dividing lines c = 0.2 and w = 0.012 (see Cassano et al. 2010b, for details).The fraction of clusters with RHs has a tendency to increase from the top left corner (high c, low w) to the bottom right corner (high w, low c), in other words, from more relaxed to more disturbed systems.In contrast, the fraction of clusters with NDE tends to increase towards the more relaxed systems.We also note that, as expected, relic clusters are all in the region of the most disturbed systems (see Jones et al. 2023, for details).Unfortunately, the information about the dynamics is incomplete, especially for NDE clusters.Morphological parameters are indeed available for 32% of NDE, for 85% of RH, for 60% of U, and for 62% of RR.This implies that the relative fractions of RH and NDE clusters in the morphological sample are quite different (50% and 24%, respectively) from those in the full sample (31% and 43%, respectively; see Table 3).This prevents us from drawing firm conclusions about the relative fraction of RH and NDE versus cluster dynamics.However, since the dynamical information about RHs is quite complete, we can confidently say that RHs are found preferentially in merging systems (see also the histogram in Fig. 9), in fact, ∼77% of them live in clusters with c < 0.2, and ∼80% of them live in clusters with w > 0.012 (∼72% of the RHs have both c < 0.2 and w > 0.012).In addition, the   10) increases from the bottom right side (merging) to the upper left side (more relaxed).This indicates that LOFAR is able to detect RHs in less disturbed systems.Although we do not have spectral information about these newly detected RHs, we speculate based on our framework that the majority of them are likely characterised by very steep spectra (Sect.7; Fig. 5).A possible contamination can be due to the presence of MHs in relaxed clusters.However, to our knowledge, only Abell 1068 (the cluster with the higher value of c in the right panel of Fig. 8) is an MH (Biava et al., in prep.).Additional investigation is required for the other clusters.

Summary
The most commonly accepted scenario for an explanation of the origin of RHs in galaxy clusters assumes that they are the result of particle acceleration due to turbulence produced in the ICM during cluster-cluster mergers.In addition to a connection with cluster mergers, this scenario predicts a heterogeneous population of RHs with synchrotron radio properties, spectral shapes, and luminosities that are correlated with the energetics of the merging events.For this reason, the study of the statistical properties of the RH populations has the power to constrain their origin and evolution and thus to test theoretical models.
In this paper, we used the Planck clusters in LoTSS-DR2 (Botteon et al. 2022) to carry out the first statistical study of RHs at low radio frequencies.We selected a sub-sample of 164 PSZ2 clusters above the 50% (M, z) Planck completeness line and spanning a wide redshift 0.07 < z < 0.5 range.This sample contains ∼50 RH.For the first time, this allows a statistical study of RHs in an unprecedented range of cluster masses, including clusters down to M 500 ∼ 2.5−3 × 10 14 M .It overcomes the wall of M 500 ∼ 6 × 10 14 M that limited previous statistical studies (e.g.Cassano et al. 2013;Cuciti et al. 2021a).The statistics is sufficient to start investigating the dependence of the RH properties on cluster mass and redshift.However, due to the (M, z) dependence of the Planck selection function, our measurements always entail mixed information of mass and redshift.
We compared the occurrence of RHs in this sample with that derived from a sample of PSZ2 clusters with M ≥ 6 × 10 14 M and z ≤ 0.35 observed with GMRT at 610 MHz (Cuciti et al. 2021b).In the same mass and redshift range, the fraction of clusters with RHs increases from ∼45% in the GMRT sample to ∼70% in the LOFAR sample (Sect.7).The observed increase is in line with predictions of the re-acceleration scenario, which implies that more RHs should be visible at lower frequency because their spectra are very steep.This sample offers the unique opportunity to test theoretical models in an uncharted range of mass and redshift and of doing this at low radio frequency for the first time.For this reason, we quantitatively compared the model expectations with our LOFAR observations.We used semi-analytic models developed in the framework of the merger-driven turbulent re-acceleration scenario (Cassano & Brunetti 2005;Cassano et al. 2006Cassano et al. , 2010a) ) to derive the expected properties of the RH population in the PSZ2 clusters.We adopted the set of values of model parameters that was previously assumed in a series of papers (see Sect. 4).By using the observed mass and redshift limit of the observed sample and by normalising the number of clusters in the theoretical model to match the observed number of clusters, we showed that we can reproduce the integral number of RHs (∼40−70 expected RH), their flux density, and their redshift distributions (Sect.5).Using the same modelling, we predict that about 100−200 RH may be detected in PSZ2 clusters (above the 50% completeness line and z ∼ 0.07−0.5)by the full LoTSS.
A clear expectation of this model is that the fraction of clusters with RHs increases with the cluster mass.Although this has already been tested at higher frequencies for massive clusters (see e.g.Cuciti et al. 2021b), here we are in the position to investigate for the first time the occurrence of RHs with cluster mass as observed at low radio frequencies and in an unprecedented range of cluster masses.We divided our sample into three redshift bins.For each redshift bin, we measured the fraction of A43, page 10 of 12  clusters with RHs, f RH , in the high-and low-mass sub-samples (Table 1).The RH occurrence clearly increases with the cluster mass at a fixed redshift.This agrees with expectations derived using the re-acceleration scenario.
Although the statistical model we used is rather simplified, we showed that it provides an excellent description of the observed properties of the RH population in the PSZ2 clusters in the LoTSS DR2.It reproduces their observed integral number, redshift, and flux density distributions very well, and it also explains the increase in RH fraction with cluster mass and at low radio frequencies.We note that the model, with the same set of values of parameters as was adopted in the current paper, has been shown to explain the RH fraction measured in a sample of massive high-z clusters (z 0.6−0.9)observed with LOFAR and followed-up at higher frequencies with the uGMRT (Di Gennaro et al. 2021b).About 50% of the LOFAR-detected RHs were found to be characterised by very steep radio spectra, in line with model expectations.Future follow-up observations at higher frequency of all RHs of the LoTSS-DR2 PSZ2 sample would definitively test this expectation, making the golden test of the scenario.Exploiting the full range of model parameters may also provide interesting indications.In the near future, as soon as the full LoTSS is completed, extensive simulations will allow us to identify a meaningful range of values of model parameters that match the data best.
We used the current sample to investigate the RH-cluster merger connection at low radio frequencies for the first time.Observations at low radio frequency are expected to find RH also in less strongly disturbed systems.Although we showed that the morphological information is not complete for our sample (especially for NDE clusters), the fraction of clusters with RHs has a tendency to increase towards more disturbed systems.In addition, the fraction of newly detected RHs by LOFAR increases from merging to more relaxed systems, which indicates that LOFAR starts to observe RHs in less strongly disturbed systems and might detect RHs with very steep radio spectra.

Fig. 1 .
Fig. 1.Redshift-mass distribution of PSZ2 clusters within the DR2 area up to z = 0.5.The Different colours show the radio classification: N/A (open black dots), NDE and RR with no RH (green dots), RH, cRH, RH* e cRH* (red dots), U (yellow dots).The 50% Planck completeness line is also show (black dotted line).

Fig. 2 .
Fig. 2. Radio power of haloes (red points) and upper limits (green triangles) as a function of redshift.The minimum radio power derived from Eq. (2) is shown as a black line (parameters are the same as in the figure panel).The blue line has been obtained by combining the P 150 MHz −M 500 best-fit correlation and the 50% Planck completeness (M, z) line in Fig. 1.

Fig. 3 .
Fig. 3. Model expectations at 150 MHz.Left panel: expected fraction of clusters with RHs with ν s ≥ 150 MHz as a function of the cluster mass.Middle panel: expected fraction of RH with very steep radio spectra (ν s < 600 MHz) as a function of the cluster mass.Right panel: RH luminosity function at ν o = 150 MHz.In all panels, the lines refers to z = 0−0.1 (black line), 0.2−0.3(blue line), and 0.5−0.6 (red line).
Fig. 3 (left panel) we plot the expected fraction of RHs with ν s ≥ 150 MHz (black upper line) as a function of the cluster virial mass and at different redshift (see the figure legend and caption).This is obtained by assuming the reference set of model parameters (i.e.B = 2 µG, η t = 0.2).At each redshift, the fraction of clusters hosting RHs with ν s ≥ 150 MHz increases with the cluster mass.In Fig. 3 (central panel), we report the fraction of RHs with ν s ≥ 150 MHz that in homogeneous models would have very steep radio spectra, specifically those that have 150 < ν s < 600 MHz.The percentage of RHs with a very steep spectrum is a strong function of the cluster mass.It decreases rapidly for high-mass clusters and also depends on the cluster redshift.For instance, we found that in clusters with virial mass ∼8 × 10 14 M at z 0.05 (i.e.M 500 ∼ 4 × 10 14 M ), ∼80−90% of the clusters with ν s ≥ 150 MHz RH have a steep radio spectrum, while at the same redshift, this percentage becomes 40% for clusters with M v ∼ 1.4 × 10 15 M (i.e.M 500 ∼ 7 × 10 14 M ).

Fig. 4 .
Fig. 4. Expected (red line and region) and observed (black dots) number of RHs in PSZ2 clusters with a flux density greater than f 150 and within z ≤ 0.5 in the LoTSS-DR2 area.The black arrow shows the upper boundary to the total number of RHs (due to classification uncertainty), and the bars on the black dots are obtained by a Monte Carlo procedure and represent the errors due to the uncertainties on the flux densities of RHs.

Fig. 5 .
Fig. 5. Number of RH as a function of redshift.Left panel: integral number of RHs within a given redshift N H (<z).The black arrows show the upper boundaries to the total number of RH within that redshift bin that are obtained by adding RH and U cases.Right panel: number of RHs in redshift bins N H (z, ∆z); the percentage of RHs with ultra-steep spectra (i.e.ν s < 610 MHz) is reported in each redshift bin.In both panels, red lines and regions are the expected values, and black line and dots are the observed values.

Fig. 6 .
Fig. 6.Left panel: M 500 vs z distribution of the cluster sample.The rectangular regions show the three redshift bins, and within each of them, the dashed line divides the sample into high-mass (above) and low-mass (below) sub-samples.Right panel: fraction of clusters with RHs, f RH , derived in the high-mass sample (red shadowed region) and in the low-mass sample (black shadowed region).For each bin, the reported minimum and maximum values of f RH are f RH (MC) with U = no RH case and f RH (MC) with U = RH case, respectively (see Table1).

Notes.
From left to right, the diving mass between the low-mass and high-mass sub-samples is in each redshift bin M 500 = 3.71 × 10 14 M , 4.81 × 10 14 M , and 5.68 × 10 14 M (in the = no RH cases); M 500 = 3.71 × 10 14 M , 4.81 × 10 14 M , and 5.54 × 10 14 M (in the U = RH cases).The main uncertainty for f RH (obs) is driven by the RH classification (upper vs. lower panel).

Fig. 7 .
Fig. 7. Comparison between expected fractions of clusters with RH (black lines) and observed f RH (red regions) in the three (redshift, mass) bins outlined in Fig. 6.Low-mass bin values are reported in the left panel, and high-mass bin values are plotted in the right panel.For each bin, the reported observed minimum and maximum value of f RH is f RH (MC) with U = no RH case and f RH (MC) with U = RH case, respectively (see Table1).

Fig. 8 .
Fig. 8. Left panel: mass-z distribution of the clusters with z = 0.07−0.5 above the 50% Planck completeness line: coloured dots show the radio classification of clusters (see figure legend in the panel) that have information about their dynamical status; grey dots are the other clusters of the sample without X-ray observations.Right panel: c−w morphological diagram for the clusters with available X-ray Chandra and/or XMM-Newton data (represented by coloured dots in the right panel).Vertical and horizontal dashed lines are adopted from Cassano et al. (2010b) and are c = 0.2 and w = 0.012.

Fig. 10 .
Fig. 10.c−w morphological diagram for the clusters with RHs, new RHs, i.e.RHs detected for the first time by LOFAR are shown with a different colour (cyan dots).Fractions of newly detected RHs are reported in different c−w regions.Vertical and horizontal dashed lines are adapted from Cassano et al. (2010b) and are c = 0.2 and w = 0.012.

Table 1 .
Observed fraction of clusters with radio haloes.

Table 2 .
Fraction of clusters with radio haloes in high-mass clusters.Ranges in the observed f RH values were obtained considering the U cases once as RH clusters and once as clusters without RH.many RHs discovered by LOFAR do not have follow-up observations at other frequencies.At the same time, in the framework of this simplified model, NDE clusters are the systems with ν s < ν o .As a consequence, RHs should become visible in a fraction of these systems in sensitive observations at even lower frequencies, for instance in future LOFAR surveys, such as LoLSS (de Gasperin et al. 2021) and LoDeSS (the LOFAR Decameter Sky Survey, van Weeren et al., in prep.).

Table 3 .
Relative fraction of NDE and RH clusters in the morphological sample and in the full sample.