Efficiency of Turbulent Reacceleration by Solenoidal Turbulence and Its Application to the Origin of Radio Megahalos in Cluster Outskirts

Recent radio observations with the Low Frequency Array (LOFAR) discovered diffuse emission extending beyond the scale of classical radio halos. The presence of such megahalos indicates that the amplification of the magnetic field and acceleration of relativistic particles are working in the cluster outskirts, presumably due to the combination of shocks and turbulence that dissipate energy in these regions. Cosmological magnetohydrodynamical (MHD) simulations of galaxy clusters suggest that solenoidal turbulence has a significant energy budget in the outskirts of galaxy clusters. In this paper, we explore the possibility that this turbulence contributes to the emission observed in megahalos through second-order Fermi acceleration of relativistic particles and magnetic field amplification by the dynamo. We focus on the case of A2255 and find that this scenario can explain the basic properties of the diffuse emission component that is observed under assumptions that are used in previous literature. More specifically, we conduct a numerical follow-up, solving the Fokker–Planck equation by using a snapshot of an MHD simulation and deducing the synchrotron brightness integrated along the lines of sight. We find that a volume-filling emission, ranging between 30% and almost 100% of the projected area, depending on our assumptions on the particle diffusion and transport, can be detected at LOFAR sensitivities. Assuming a magnetic field B ∼ 0.2 μG, as derived from a dynamo model applied to the emitting region, we find that the observed brightness can be matched when ∼1% of the solenoidal turbulent energy flux is channeled into particle acceleration.


INTRODUCTION
Galaxy clusters are filled with a hot intra-cluster medium (ICM), that has a characteristic temperature similar to the the cluster's virial temperature.This suggests that the ICM is heated by the gravitational energy released in the hierarchical merger and accretion processes of clusters (e.g., Press & Schechter 1974;Kravtsov & Borgani 2012).A fraction of the energy can also be channeled into non-thermal components, such as relativistic particles and magnetic fields.Shocks and turbulence could be favorable sites for the particle acceleration and the amplification of the field (see Brunetti & Jones 2014, for review).
Radio observations of galaxy clusters probe those nonthermal components by studying diffuse synchrotron emission of relativistic (cosmic-ray) electrons (CRe).Radio halo is a diffuse emission with an extension of ∼1 Mpc, often found in the central region of merging clusters (see van Weeren et al. 2019, for review).The radiative cooling time of CRe is significantly shorter than the time required for diffusion or advection over ∼ 1 Mpc, implying that there is an in situ mechanism that produces CRe.Reacceleration by merger-induced turbulence is the most plausible scenario (e.g., Brunetti et al. 2001;Petrosian 2001;Fujita et al. 2003;Cassano & Brunetti 2005), although cosmic-ray protons (CRp) in the ICM may be important ingredients in the physics of those phenomena (e.g., Dennison 1980;Blasi & Co-lafrancesco 1999).For example, they can provide seed CRe to re-accelerate through the hadronic pp collision with the thermal protons in the ICM (e.g., Brunetti & Lazarian 2011a;Brunetti et al. 2017;Pinzke et al. 2017;Nishiwaki & Asano 2022).It has been shown that the observed statistical properties of radio halos are in line with the reacceleration model considering the resonant interaction between compressible turbulence (Cassano & Brunetti 2005;Nishiwaki & Asano 2022;Cassano et al. 2023).Reacceleration by turbulence is also proposed to explain radio emission detected on larger scales, such as radio bridges (Brunetti & Vazza 2020) that are radio filaments connecting massive pairs in the early stage of mergers discovered by the Low Frequency Array (LO-FAR) (Govoni et al. 2019;Botteon et al. 2020a).
More recently, Cuciti et al. (2022) reported the existence of radio "mega halos" in four clusters, using the LOFAR observation.The volume of the mega halos is almost 30 times larger than that of radio halos, suggesting that the entire volume of the cluster is filled with CRe and magnetic field.The radio power of mega halos is comparable to or even larger than that of classical halos.The detection of synchrotron radiation at the large distance (1 -2 Mpc) from the cluster center also constrains the magnetic field strength in this region.Since the pressure of non-thermal components, including magnetic field and CRe, should be smaller than the thermal one as indicated from the observations and numerical simulations (e.g., Vazza et al. 2016;Eckert et al. 2019), the magnetic field should be in the range of 0.1µG ≲ B ≲ 1.7µG (Botteon et al. 2022).The lower bound (0.1µG) is at least one order of magnitude larger than the value expected by the compression of primordial fields.One possible mechanism of this nonlinear amplification of the field is dynamo in a turbulent medium.
Cosmological simulations of galaxy clusters suggest that turbulence and shocks driven by continuous accretion of matter fill the entire volume of the cluster up to the virial radius (e.g., Vazza et al. 2011;Nelson et al. 2014;Miniati 2015;Steinwandel et al. 2023).As in the cluster center, the ICM in the outskirts is a weaklycollisional plasma, and the perturbations would cause instabilities that effectively reduce the mean free path (mfp) of thermal protons (e.g., Schekochihin et al. 2005;Kunz et al. 2011;Brunetti & Lazarian 2011b), potentially establishing a well-developed inertial range (e.g., Schekochihin et al. 2009).
Indeed, the infall of the clumps of mass drives the turbulence with a typical scale of a few 100 kpc in the cluster outskirts (e.g.Vazza et al. 2017).The timescale of the turbulent cascade, t cas ∼ 600 Myr (e.g., Brunetti & Lazarian 2007), is much shorter than the Hubble time, which allows the turbulent dynamo to work in that region.
Nowadays, the best studied case of the cluster hosing diffuse radio emission on the entire cluster volume is the case of Abell 2255 (hereafter A2255) (Botteon et al. 2022).A2255 is a nearby (z = 0.0806) cluster, which shows a complex dynamical state in optical and X-ray observations (e.g., Burns et al. 1995;Yuan et al. 2003;Golovich et al. 2019;Feretti et al. 1997;Akamatsu et al. 2017).The cluster is known to host a radio halo and radio relics, and they have been studied in wide range of frequencies (e.g., Jaffe & Rudnick 1979;Feretti et al. 1997;Govoni et al. 2005;Pizzo et al. 2008;Pizzo & de Bruyn 2009;Botteon et al. 2020b).
More recently, LOFAR observations found that the cluster hosts diffuse emission extending in very large scales and enveloping the classical halo and relics (Botteon et al. 2022).Such emission is complex, showing a number of relic-like features embedded in a truly diffuse component, and a spectral index distribution between 40 -144 MHz ranging from 0.6 to 2.5.The flat spectrum emission is coincident with the radio relics located in the north and southwest parts of the cluster, while the steep emission is associated with the diffuse component.
In this paper, we attempt to explore the possibility that turbulence contributes to the observed emission via magnetic field amplification and particle reacceleration.We use cosmological MHD simulation of a cluster to examine the turbulent and magnetic fields in the cluster outskirts.In Sect.2, we describe the setup of the MHD simulation.In Sect.3, we review the reacceleration model and claim that that is compatible with the observed spectrum of the mega halo.In Sect.4, we numerically solve the Fokker-Planck (FP) equation, considering the distribution of turbulence and the projection along the line of sight.The limitations of our models are discussed in Sect. 5. Finally, we summarize the results in Sect.6.

COSMOLOGICAL MHD SIMULATION
To examine the property of the ICM and turbulence in cluster outskirts, we use a snapshot of a high-resolution cosmological ideal MHD simulation of a massive galaxy cluster, produced with grid code ENZO (Bryan et al. 2014).We use the same simulated cluster as in Botteon et al. (2022), which has a mass of M 200 = 8.65×10 14 M ⊙ at z = 0.This simulation includes eight levels of Adaptive Mesh Refinement (AMR) to increase the spatial and force resolution within the virial volume of the cluster, reaching a peak spatial resolution of ∆x = 3.95 kpc/cell (comoving) (Vazza et al. 2018).
The simulation was started assuming a uniform "primordial" seed magnetic field of B 0 = 0.1 nG (comoving) at z = 40.The low redshift properties of the magnetic field in the cluster volume are found to be fairly independent of the exact origin scenario, due to the effect of efficient small-scale dynamo amplification (Vazza et al. 2018).At the late stage of the cluster evolution, most of the central volume within < 1 Mpc is simulated with the finest resolution ∆x = 3.95 kpc/cell.The virial volume of clusters is refined at least up to ∆x = 15.8 kpc/cell (Domínguez-Fernández et al. 2019).At 1-2 Mpc from the center, the spatial resolution is comparable to, or coarser than, the MHD scale, l A ∼ 1 kpc, where the turbulent velocity becomes comparable to the Alfvén velocity, so the field amplification may be underestimated in the simulation.Thus, we evaluate the field strength in this region in a post-process (see Sect. 3 for the detail).
The turbulent kinetic flux is calculated with smallscale filtering explained in Vazza et al. (2017), which allows us to reconstruct (and remove) the contribution from shock waves to the total turbulent energy budget.For each cell of the simulation, the dispersion of the velocity field σ v (L) is measured within a scale L. The outer scale of the turbulence Λ is defined as the scale where the change in σ v (L) with increasing L becomes sufficiently small (see Vazza et al. 2017, for the detail).The turbulent kinetic energy flux in each simulated cell can be calculated as: where ρ is the gas mass density in the cell and ∆x is the cell size.Based on previous works, we assume that the turbulent spectrum in the inertial range roughly follows the Kolmogorov scaling (e.g.Vazza et al. 2011Vazza et al. , 2017)), and so the value of F turb is insensitive to the specific value of L as long as it is in the inertial range.To study the quantities in the cluster outskirts, we extract a box of a region of 1.6 Mpc on a side located at 1.2 Mpc from the center of the simulated cluster.Each cell has a volume of 16 3 kpc 3 , and the box is composed of 10 6 cells.In this region, the outer scale Λ is found to be typically Λ ≈ 200 kpc.In the following, we use the turbulent velocity measured at L = 160 kpc.Fig. 1 gives the visual impression of the line of sight gas velocity and of the gas temperature in the volume of the simulated clusters, in particular in the box used to extract the physical parameters used in our Fokker-Planck calculations.The selected box is at the cluster periphery, in a region only partially detectable through typical X-ray observations, and along the direction of a prominent filamentary accretion, but with-out the presence of massive substructures, or prominent shock waves.
Fig. 2 (left panel) shows the histogram of turbulent kinetic flux in the extracted region.We decompose the turbulent velocity into solenoidal (∇ • v = 0) and compressible (∇ × v = 0) modes, using the procedure of Vazza et al. (2017).We find that the solenoidal mode typically has a larger kinetic flux than the compressible mode, in line with previous simulations (e.g., Miniati 2015;Porter et al. 2015;Vazza et al. 2017).The mean value of the solenoidal kinetic flux per unit volume is F turb = 1.1 × 10 44 erg/s/Mpc 3 , which is almost 10 times larger than the compressible component.In the right panel, we show the cumulative number of cells that have turbulent kinetic flux larger than F turb .More than 30% of the cell has the solenoidal turbulent flux larger than 10 44 erg/s/Mpc 3 , while the fraction decreases to ≈ 3% in the case of the compressible mode.Note that we are considering the sector where the feature of mass accretion from the nearby cluster can be seen (Fig. 1) and it is more turbulent than other sectors of cluster outskirts (see also Sect.5).For comparison, we study other regions with the same volume and the distance from the cluster center and find a factor ∼5 smaller turbulent flux, whereas the solenoidal mode dominates the compressible one in every region.In the following sections, we focus on the solenoidal velocity to calculate the acceleration efficiency and the dynamo field.Note that F turb denotes only the volumetric turbulent kinetic flux of the solenoidal component.

CR REACCELERATION AND FIELD AMPLIFICATION BY SOLENOIDAL TURBULENCE
Turbulence driven through the formation process of galaxy clusters is typically sub-sonic (M s < 1) and super-Alfvénic (M A > 1) (Brunetti & Lazarian 2007).A fraction of the turbulent energy can be converted into non-thermal components such as cosmic rays and the magnetic field through stochastic acceleration and turbulent dynamo.Turbulence in astrophysical environments may accelerate particles through different mechanisms, including resonant and non-resonant mechanisms (e.g., Ptuskin 1988;Schlickeiser & Miller 1998;Cho & Lazarian 2006;Brunetti & Lazarian 2007;Lynn et al. 2014;Brunetti & Lazarian 2016;Bustard & Oh 2022;Lemoine 2021;Lazarian & Xu 2023).Recent MHD simulations of galaxy clusters suggested that the turbulence in the ICM is dominated by the solenoidal (incompressible) mode (e.g., Miniati 2015;Vazza et al. 2017).One possible acceleration mechanism working in incompressible turbulence is the acceleration due to the interaction   between magnetic field and particles diffusing in super-Alfvénic turbulence (Brunetti & Lazarian 2016;Brunetti & Vazza 2020).
In turbulent reconnection theory (Lazarian & Vishniac 1999), the Alfvén scale l A ≡ LM −3 A (for the Kolmogorov scaling) is the dominant scale, where the reconnection speed may reach v rec ∼ v A .In the regions where the magnetic field dissipates due to the reconnection, the particles trapped in contracting islands gain energy through a mechanism that is similar to the firstorder Fermi mechanism (Kowal et al. 2012;del Valle et al. 2016).On the contrary, the particles are expected to cool in the regions where the dynamo is efficient and the magnetic field lines diverge.
According to Brunetti & Lazarian (2016) particles diffusing through this complex pattern experience a second-order acceleration mechanism.In this mechanism, the diffusion coefficient in the momentum space is where p is the momentum of the particle, A is the plasma beta with γ ad = 5/3 is the adiabatic index, and ψ ≡ λ mfp /l A with λ mfp is the mean free path (mfp) of the particle.The mfp is an important parameter in the model as it determines the spatial diffusion coefficient, and consequently the efficiency of particle crossing regions of the size l A .Combining the requirement that the fractional change of momentum in each scattering should be ∆p/p ≪ 1 and the effect of particle scattering due to mirroring in a super-Alfvénic flow, ψ should satisfy 0.01 ≪ ψ ≲ 0.5; ψ ∼ 0.5 is the reference value that has been motivated and used in the previous studies (Brunetti & Lazarian 2016;Brunetti & Vazza 2020).
We consider that the solenoidal turbulence with super-Alfvénic injection velocity shows a power spectrum with a well-established inertial range.In such a situation, a fixed fraction (η B ≈ 0.05) of the turbulent energy flux is consumed through the amplification of the magnetic field (e.g., Cho et al. 2009;Beresnyak 2012;Xu & Lazarian 2016).MHD simulations of galaxy clusters succeed in resolving the small-scale turbulent dynamo in the central region (e.g., ZuHone et al. 2011;Vazza et al. 2018;Domínguez-Fernández et al. 2019;Steinwandel et al. 2022), but this effect is quenched in the cluster periphery due to the limited resolution.Thus, we adopt the same procedure as Brunetti & Vazza (2020) and estimate the dynamo-amplified field in post-processing, in the output of the MHD simulation (Sect.2).Assuming that a fraction η B of the turbulent kinetic energy is channeled into the field amplification, the field strength can be estimated as where t eddy = L/δv turb is the eddy turn-over time.In this case, M A and η B are related as . In the simulations of galaxy clusters, the value of η B is as small as η B ≈ 0.05 (Beresnyak & Miniati 2016;Vazza et al. 2018).Using η B = 0.05, we obtain the mean field strength of ⟨B⟩ ∼ 0.24µG in the simulated box, which is almost 2 times larger than the mean value of the field strength directly measured in the simulation (⟨B sim ⟩ = 0.11µG).The plasma beta is β pl ≈ 200 − 500, which is slightly larger than that in the cluster center (e.g., Brunetti & Lazarian 2007).Adopting the dynamo field (Eq.( 4)), D pp (Eq.( 2)) can be expressed as a function of η B as (Brunetti & Vazza 2020) (5) The spectral index map of A2255 in Botteon et al. ( 2022) is a mixture of flat (α syn ∼ 1.0) and steep (α syn ∼ 2.0) regions.The flattest index is coincident with the arc-shaped structure, suggesting the particle energization by shocks.On the other hand, the emission with steep radio indices comes from the diffuse envelope permeating the cluster volume (see Fig. S2 of Botteon et al. 2022) .The mean value of the index is around α syn ≈ 1.6.If we assume a scenario where particles emit in a homogeneous field, those steep indices suggest that the observing frequency is close to the steepening frequency, at which the spectrum starts to decline rapidly due to the cooled spectrum of CRe.
In the turbulent reacceleration model of CRe, one can define the break energy γ b as the energy where the cooling timescale t cool becomes comparable to the acceleration timescale t acc .As shown in Cassano et al. (2010), the steepening frequency ν s in the synchrotron spectrum appears at a factor ξ ∼ 5 − 7 times larger than the break frequency ν b corresponding to γ b , i.e., ν s = ξν b .The cooling time of ultra-relativistic CRe due to synchrotron and inverse-Compton radiation can be written as t cool = 6πm e cγ −1 /(σ T (B 2 + B 2 CMB )), where γ is the Lorenz factor of the particle, σ T is the Thomson cross section and B CMB = 3.25(1 + z) 2 µG is the inverse-Compton (IC) equivalent field.Note that 100 MHz emission is mainly produced by the synchrotron radiation of CRe with p = γ 2 − 1 ∼ 10 4 , and the Coulomb loss is negligible for those CRe.Using the magnetic field estimated from the dynamo model of Eq. ( 4), the cooling timescale at ν b can be estimated as where µ is mean molecular weight, ν o = ν s /(1 + z) is the observing frequency.The cooling is dominated by the radiation through the IC scattering (B 2 ≪ B 2 CMB ) in the cluster outskirts, so we neglected the B 2 term in the denominator.The acceleration timescale should be comparable to this value to sustain the emission observed at the LOFAR frequency.
Comparing Eq. ( 6) with the acceleration timescale t acc = p 2 /(4D pp ) obtained from Eq. ( 5), one can re- late two fundamental parameters, ψ and η B .Fig. 3 shows the relation between η B and ψ in the case of ν o = 50 MHz and z = 0. We adopted the typical value for the ICM density and temperature in cluster outskirts; n ICM = 3 × 10 −5 cm −3 and T ICM = 4 keV.Adopting η B ∼ 0.05 and the typical turbulent flux found in MHD simulation, F turb ≈ 10 44 erg/s/Mpc 3 (Fig. 2), we find that the observed steep spectrum at LOFAR frequencies can be explained with the mfp of CRe ψ ∼ 0.5 of the Alfvén scale.This value is consistent with that adopted by Brunetti & Lazarian (2016) and Brunetti & Vazza (2020) to explain radio halos and bridges.The mfp comparable to l A in turbulence reconnection is supported by the studies of tracers in MHD simulations (Kowal et al. 2011(Kowal et al. , 2012)).1 One can calculate the number of cells that are important for the radio emission by comparing t acc and t cool measured in the extracted cubic volume.Using Eqs. ( 5) and ( 6), the scaling relation for the fraction of those timescales can be expressed as In Fig. 4, we show the cumulative fraction of t acc /t cool .When t acc < t cool , the acceleration mechanism of Eq.( 2) can efficiently re-accelerate CRe that radiate synchrotron emission in the LOFAR frequencies.
The fraction of cells is ∼ 15% adopting a reference value ψ = 0.5, while it increases to ≳50 % for ψ < 0.2−0.3.The dependence on η B is rather weak (Eq.( 7)).The fraction ranges form 25% to 10% for η B = 0.01−0.1,for a fixed ψ = 0.5.Such a significant fraction implies a volume-filling emission.One also needs to take care about the duration of reacceleration, since the reacceleration is important only when it lasts for several times longer than t acc .The filling factor of the radio-emitting cells will be revisited in Sect. 4.

FOKKER-PLANCK SIMULATION
Next, we study the synchrotron spectrum in our dynamo reacceleration model with a numerical calculation.We consider the situation that the pre-existing population of seed CRe is re-accelerated by incompressible turbulence and produces the observed radio emission by the synchrotron radiation.As in Sect.2, we consider the volume of 1.6 3 Mpc 3 located 1.2 Mpc from the center, and the distribution of the seed CRe is proportional to the ICM number density.We solve the Fokker-Planck equation of CRe of the following form (e.g., Cassano & Brunetti 2005): where N e (p, t) is the spectrum of CRe at time t, ṗ represents the momentum loss per unit time, D pp is the momentum diffusion coefficient, and Q e is the injection spectrum of CRe.The momentum loss rate ṗ includes the effect of radiative(synchrotron and IC) and Coulomb losses, and D pp is calculated with Eq. ( 2).We neglect the injection of the secondary electrons from inelastic pp collisions of cosmic-ray protons.Due to the low ICM density in the cluster outskirts, the contribution of the secondary electrons is expected to be smaller than that in classical radio halos (see Appendix A for further discussion).
We neglect the term related to the spatial transport of CRs in Eq. ( 8), since it cannot be properly followed with a snapshot of the MHD simulation.However, we are considering the case where the mfp of CRe is comparable to the Alfvén scale (l A ∼ 1 kpc), and this leads to a large value of the spatial diffusion coefficient D ∼ 10 31 cm 2 /s, as shown in Eq. ( 3).The detailed investigation of the effect of spatial transport in the reacceleration model is left for future works using a Lagrangian tracer method.In this exploratory work, we adopt following two approaches: a one-zone approximation (Sect.4.1), where we adopt the average value of the physical quantities found in the snapshot, and the calculation considering the variation of F turb in each cell and the integration along the line of sight (LOS) (Sect.4.2).

One zone calculation
We first discuss the energy density of CRe required to reproduce the observed emission under the onezone approximation.We adopt the mean values in the simulation box for the physical quantities; n ICM = 6 × 10 −5 cm −3 , T ICM = 4.5 keV, v sol = 320 km/s (M s ≈ 0.4), and t eddy = 0.53 Gyr.We adopt η B = 0.05, which corresponds to B = 0.24µG (Eq.( 4)).We summarize the values of parameters in Tab.4.1.As an initial condition, a cooled spectrum of seed CRe is calculated by integrating Eq.( 8) for 2 Gyr with D pp = 0 and Q e (p, t) ∝ p −αinj δ(t), where α inj = 2.2 and δ(t) is the Dirac delta function.Due to the radiative and the Coulomb cooling, the seed CRe has steep cut-offs around p min ∼ 10 and p max ∼ 10 3 (Fig. 5 right panel).The normalization of the initial spectrum is treated as a free parameter.
As seen from Fig. 3, for M s ∼ 0.4 − 0.7, a value ψ ∼ 0.3 − 0.5 is compatible with the observed emission.This is confirmed by numerical calculation in Fig. 5 (left panel), where we plot the synchrotron spectrum in the one-zone model with a dashed line.While the spectral shape is mainly determined by the parameters η B and ψ (Eq.( 7)), the normalization (brightness) depends on  b The magnetic field is calculated with Eq. ( 4).
the initial spectrum N e (p, 0) and the duration of reacceleration T dur .In this section, we assume that T dur is sufficiently long and the CRe spectrum reaches a steady state due to the balance between the reacceleration and the radiative cooling (Eq.( 8)).We find that this steady state is achieved at T dur ≥ 3 Gyr.The duration of T dur = 3 Gyr is comparable to the dynamical timescale of the ICM in the simulated region, so the mean value of the turbulent flux F turb can evolve in this timescale.Thus, the assumption of constant F turb during the FP calculation should be considered as a rough approximation.
We determine the normalization of the initial spectrum as the synchrotron brightness at 49MHz matches the observed value.The data points in Fig. 5 show the range of the brightness of the diffuse envelope measured in 1.5-2 Mpc from the center of A2255 (Botteon et al. 2022).
We define the efficiency of the reacceleration, η acc , as the fraction of the turbulent kinetic flux that turns into the increase of the CRe energy density and the radiation: where ϵ CRe is the CRe energy density and ε rad = ε syn + ε IC is the sum of the frequency-integrated emissivities of synchrotron and IC radiations.In the steady state of N e , the first term on the right-hand side is negligible.The synchrotron emissivity ε syn is calculated in the range of 10 4 − 10 10 Hz, and ε IC is calculated with In Fig. 5 (dashed line), we show the CRe energy spectrum in the steady state.The energy is dominated by 100 MHz-emitting particles and the particles with p < 10 3 are not energetically important.We find ϵ rad ≈ 7×10 −30 erg/s/cm 3 and η acc ≈ 2%, which is consistent with the estimate in Botteon et al. (2022).Due to the assumption of stationary conditions, this result is independent of the initial spectrum.
We find that the spectrum in Fig. 5 is reasonably reproduced by considering a range of values of ψ around the reference value ψ ∼ 0.5.On the other hand, the assumption of a mfp that is significantly reduced, ψ < 0.2, generates spectra that are too hard compared to the observation.Note that this result is based on the assumption that the turbulence and physical parameters in the simulated ICM are representative of the external regions in A2255.

projection along the LOS
Next, we consider the distribution of turbulent energy in the box and study how the emission from highly turbulent cells affects the spectrum integrated along the LOS.To consider the CRe spectrum in each simulated cell, we calculate the FP equation for various F turb found in the simulated box.We make a histogram of F turb with 160 bins equally spaced in the logarithmic scale in the range of 10 38 erg/s/Mpc 3 < F turb < 10 46 erg/s/Mpc 3 (Fig. 2).For simplicity, we adopt the same values of dynamo B (Eq. ( 4)) and D pp (Eq.( 5)) for the cells in the same F turb by calculating the mean values of ρ ICM and v turb in each F turb bin.The synchrotron brightness is calculated by integrating the emissivity of 100 cells (corresponds to 1.6 Mpc) along the axis of the simulation.In one projection, there are 100×100 LOS.
We assume the same initial spectrum as in Sect.4.1, i.e., the spectrum with cut-offs at p min ∼ 10 and p max ∼ 10 3 .As in the previous section, we consider a sufficiently long T dur = 3 Gyr to ensure that the steady state due to the balance between cooling and reacceleration is achieved in many of the cells.Note that the turbulent flux in each cell would change in a few eddy turnover time, t eddy ∼ 0.3 − 1 Gyr.When T dur is short and the steady state is not achieved, the result may depend on the initial condition.This point is further discussed in Appendix B, where we show that the observed emission can be explained with a different combination of T dur and the initial condition.
The non-dimensional parameters in our model, ψ (Eq.( 2)) and η B (Eq. ( 4)), are assumed to be constant over cells.As a test case, we adopt the same η B = 0.05 as the one zone calculation.We use a snapshot of the simulation and do not consider the evolution of the background fluid.We assume that the initial energy density of the CRe in each cell is proportional to the thermal energy density i.e., ϵ CRe (t = 0) ∝ ϵ ICM .We neglect the particle transport between cells during the calculation.Those simplification reduces the computational cost and enables us to calculate the spectrum in every cell in the box for several Gyrs.The study on the impact of CR transport is left to future works featuring the Lagrangian tracer approach (Sect.5).
We calculate the distribution of the spectral index and compare it with the observation by Botteon et al. (2022).The spectral index and its statistical properties would depend on the beam size and sensitivity of the observation, so we consider those of LOFAR.The beam size in Botteon et al. (2022) is 35", which corresponds to almost 3 times the cell size at the redshift of A2255, so we calculate the spectrum of one beam by summing up the intensities of 3 × 3 neighboring LOS.Each simulated beam consists of the emission from 900 cells.Following Botteon et al. (2022), we introduce a cut in our simulated data, neglecting the beams that has 145 MHz brightness below 2σ HBA , where σ HBA = 200µJy is the r.m.s.noise per beam of the LOFAR HBA band.
In Fig. 5, we show the typical synchrotron spectrum of beam detectable with the LOFAR sensitivity (black); here we specifically use ψ = 0.6.Since more turbulent cells contribute more than less turbulent cells the results are not exactly consistent with those in the onezone model (Sect.4.1) or with Fig 3 .In fact, the model considering the LOS integration of ∼ 10 3 cells (in each beam) requires a slightly larger value of ψ and consequently a slightly less efficient acceleration mechanism.
We decompose the spectrum into the contributions from F turb > 5 × 10 44 erg/s/Mpc 3 (orange) and F turb < 5 × 10 44 erg/s/Mpc 3 (blue) cells.We confirm the trend seen in Brunetti & Vazza (2020); at higher frequencies (ν > 100MHz), the highly turbulent cells dominate the emission, which occupies only a small fraction (4.5%) of the volume, while the emission is more volume-filling at lower frequencies.In Fig. 6, we show the contribution from the turbulent cells as a function of frequency.CRe spectra are compared in Fig. 5. Unlike one-zone calculation, the overall spectrum has a bump around p ∼ 100.CRe in the cells with smaller turbulent energy dominates the CRe energy, although they do not significantly contribute to the emission at ∼ 100 MHz (Fig. 6).The typical momentum of CRe that corresponds to 100 MHz emission shifts to p ∼ 2 − 5 × 10 3 , since the magnetic field in F turb ≈ 5×10 44 erg/s/Mpc 3 cells is a factor ∼2 stronger than B = 0.24 µG adopted in the one-zone model.As discussed in Appendix B, the contribution from low F turb can be larger and the emission becomes more volume-filling under different assumption on the initial condition.
In reality, the CR distribution would be smoothed by diffusion and/or streaming, and CRs would experience multiple reacceleration within the dynamical time.In such a situation, the gap between cells with large F turb For comparison, the spectrum in the one-zone model is shown with the dashed line.The data points show the typical value of the brightness measured at 1.5-2 Mpc from the cluster center (Botteon et al. 2022).The error bars show the ranges of the brightness in this region including 1σ error.We adopt ψ = 0.3 for the one-zone model, while ψ = 0.6 and for the cell-wise calculations.In both models, we assume T dur = 3 Gyr.The beam spectrum is decomposed into the contributions from the cells with F turb < 5 × 10 44 erg/s/Mpc 3 (blue) and F turb > 5 × 10 44 erg/s/Mpc 3 (red).Right: CRe spectra in the same calculations.Spectra before and after the reacceleration is shown with thin and thick black lines, respectively.The dashed lines show the spectra in the one-zone model.and small F turb would be reduced and the emission becomes more volume-filling.This point would be further studied in future studies with a Lagrangian tracer method (see also Sect.5).
We find that 375 beams out of 1089 satisfy the criterion of detection, i.e., S HBA > 2σ HBA , so almost 34% of the area of emission region can be covered by the LOFAR sensitivity.This fraction increases with the efficiency of reacceleration (smaller ψ), as the contribution by the less turbulent cells (F turb < 5×10 44 erg/s/Mpc 3 ) becomes more significant.However, we find that, for ψ ≤ 0.3, the typical spectrum starts to become flatter than the observed one.
The value of η acc (Eq.( 9)) differs in cell to cell, and it becomes η acc ≈ 0 when F turb is very small and the effect of the reacceleration is negligible.To obtain a typical value of η acc in this model, we calculate where the sum is taken for the cells in each beam (∼ 10 3 cells).We calculate the mean value of η acc for the 375 beams and find ⟨η acc ⟩ ≈ 1%.Although there is a plenty of turbulent energy in F turb > 5 × 10 44 erg/s/Mpc 3 cells, the occurrence of such cells is small (4.5%).The combination of those two result in η acc ≈ 1%, which is slightly smaller than that found in the one zone model (Sect.4.1).Fig. 7 shows the distribution of spectral index in each beam calculated with ψ = 0.6 and η B = 0.05.The black dashed line shows the result reported in Botteon et al. (2022).The mean value of the index in our calculation is ⟨α⟩ ≈ 1.7, and the difference in the results of three different projections is marginal.Our results are in line with the observation of diffuse radio emission enveloping A2255.

LIMITATIONS
Clearly, the complex morphology observed by LOFAR suggests that at least shocks and turbulence contribute   2022).In our calculation, the total number of observable beams is ≈ 300 (30% of the area), while the observational data shows the distribution of 836 pixels.For the model parameters, we adopted ηB = 0.05 and ψ = 0.6.The reacceleration is calculated for T dur = 3 Gyr, and the beam size is assumed to be 35".We tested the projection along three different axes (x, y, and z, distinguished by colors).
to the emission.A full modeling of A2255 is however beyond the aim of the present work.We have explored the possibility of turbulent reacceleration in the context of the specific model Brunetti & Lazarian (2016).We find that the steep-spectrum diffuse emission observed at 1-2 Mpc distance from the center of A2255 can be explained using parameters (η B and ψ) that are in line with previous literature (Brunetti & Lazarian 2016;Brunetti & Vazza 2020).In this case, the spatial diffusion coefficient of the CRe becomes D ∼ 10 31 cm 2 s −1 (Eq.( 3)), i.e., CRe can diffuse over 450 kpc within 1 Gyr.Within the acceleration time of radio-emitting CRe t acc ∼ 500 Myr (Eq.( 6)), the diffusion length is ∼ 300 kpc.In addition, CRe can be spread and mixed by turbulent motion.This implies that CRe can travel ∼30 cells during reacceleration, making the distribution smoother.Since the cells that efficiently accelerate CRs appear in every ∼ 3 3 cell (i.e., ∼ 4%), the diffusion is fast enough to fill the space between turbulent cells.We tested two cases: the one-zone model in Sect.4.1 and the cell-wise calculation in Sect.4.2.In the former case, we adopt the average values of the physical quantities found in the simulated box to calculate the FP equation.In the latter case, we neglected the transport within a short duration of the reacceleration, and cal-culated the reacceleration using the local value of the turbulent flux in each cell.A future study with a Lagrangian tracer method will be important to discuss the evolution of the CR spectrum and spatial distribution due to the combination of the reacceleration, the spatial diffusion, and advection.An observation with higher angular resolution would also be important to study the correlation between the flat spectrum and large turbulent kinetic flux predicted in Sect.4.2, or the gradient of the spectral index around the turbulent region due to the diffusion.
Concerning the initial condition, we have assumed that there is a cooling phase with D pp = 0 before the reacceleration starts, as in many models of the giant radio halo (e.g., Brunetti et al. 2001;Nishiwaki & Asano 2022).However, CRs can be gently re-accelerated for several Gyrs by the modest level of turbulence continuously driven by mass accretion.This would imply that there is no clear onset of the mega halo emission, unlike the classical halos which are supposed to be driven by the mergers of clusters (e.g., Cassano et al. 2016).In Appendix B, we assume a different initial condition, assuming that the reacceleration was working before the epoch of the snapshot.In this case, the observed emission can be explained by a an efficiency that is slightly reduced, η acc ≈ 1.7%.
The tracer approach is also important for that issue.Following the energy gain and loss of CRe with a tracer method in a simulated galaxy cluster, Beduzzi et al. (2023) demonstrated that ≈ 22 − 57% of the mega halo region (0.4R 500 < r < R 500 ) is filled with radio-emitting CRe.Although the details of the simulation setup and the definition of the volume filling factor are different from our study, their result suggests that the clusterscale diffuse emission is produced through the multiple episodes of turbulent reacceleration.
When extracting the peripheral region of the simulated cluster in Sect.2, we choose the particular sector where a large-scale accretion can be seen (Fig. 1), motivated by the fact that the most turbulent region dominates the emission in our model (Sect.4).We note that the typical turbulent energy can be smaller in other sectors without such an accretion feature.Considering that the properties of the gas and turbulence seen in the MHD simulation can be different from those in A2255, the parameters reported in this study are basically indicative values.

CONCLUSIONS
Recent LOFAR observations reported the presence of diffuse radio emission permeating the volume of galaxy clusters up to the virial radius.This emission is termed a mega halo, and its mechanism is still unclear.
The diffuse radio emission in A2255 has a very large extension, enveloping the example of classical halo and relics reported in previous observations (e.g., Botteon et al. 2022).The complex radio morphology is likely generated through a mix of different processes, e.g., particle acceleration by shocks and turbulence.The diffuse emission permeating the cluster volume on large scales is a candidate for turbulent reacceleration.In this paper, we have explored this hypothesis.
Recent MHD simulations observe that the turbulent energy is dominated by the solenoidal component in the cluster outskirts, so we focus on the acceleration mechanism presented by Brunetti & Lazarian (2016), where particles are accelerated by the interaction between particles and magnetic field lines diffusing in a super-Alfvénic solenoidal turbulent flow.The acceleration efficiency D pp depends on the mfp of the CR particle ψ ≡ λ mfp /l A , which is treated as a free parameter (Eq.( 2)).
We use a snapshot of a simulated cluster mimicking A2255 and study the property of turbulence and magnetic field in the cluster outskirts.We extract a cubic volume of a region in the cluster periphery, where the prominent feature of mass accretion can be seen (Fig. 1).Since the spatial resolution of the simulation is not sufficient to resolve the small-scale dynamo at the Alfvén scale, the magnetic field is estimated in a post-process.We assume that η B ≈ 0.05 of the turbulent flux is consumed as the dynamo and obtain B ≈ 0.2 µG, which is roughly two times larger than the field found in the simulation.Comparing the efficiency of reacceleration and radiative cooling, we find that the diffuse emission in A2255 can be explained with a fiducial value of ψ(∼ 0.5) that was derived from considerations based on mirroring of electrons in a super-Alfvénic flow (Brunetti & Lazarian 2016;Brunetti & Vazza 2020).
The emission spectrum is calculated by numerically integrating the FP equation (Eq.8) with the quantities found in the MHD simulation.Before the reacceleration, CRe is accumulated around p ∼ 10 2 due to the radiative and Coulomb cooling before the reacceleration.In Sect.4.1, we adopted the one-zone approximation and used the average quantities in the simulation box.We find that the reacceleration by the solenoidal turbulence is efficient enough to produce a CRe spectrum peaked at p ∼ 10 4 , corresponding to the synchrotron frequency of ∼100 MHz in the ∼ 0.1µG field.The acceleration efficiency is η acc ≈ 0.02, in line with the estimate in Botteon et al. (2022).
We also calculate the FP equation in 10 6 cells in the simulation box, assuming that the seed CRe is uniformly distributed.We find that 100 MHz emission is dominated by the cells with large turbulent kinetic flux (F turb > 5×10 44 erg/s/Mpc 3 ), which fill only a few % of the volume.Considering the LOS integration within the beam size of LOFAR, we find that a large fraction of the beams (≳ 30%) can be detected with the LOFAR sensitivity.In this model, only 1% of the turbulent energy needs to be consumed for the particle acceleration in the turbulent cells.Our model predicts that the emission will be more volume-filling when observed with higher sensitivity at lower frequencies.
The reported parameters and the derived efficiencies are indicative values, as the simulations do not necessarily reproduce the turbulent properties of A2255.For this reason, our study simply suggests that there is room for turbulent reacceleration to contribute to the observed emission.
In reality, CRe in each cell would be mixed by the streaming and/or diffusion, although we neglect those effects in the current study.As discussed in Sect.5, the mfp comparable to the Alfvén scale leads to a strong diffusion with D ∼ 10 31 cm 2 /s, and CRe can diffuse over a few 100 kpc within ∼1 Gyr.A method that can incorporate those effects, such as the application of a Lagrangian tracer method (e.g.Vazza et al. 2023), will be important for more accurate modeling.
In addition to the turbulent reacceleration, the injection of secondary leptons from the the decay of pions produced by pp collision between CRp and thermal protons is often invoked for the mechanism for the diffuse radio emission in galaxy clusters (e.g., Dennison 1980;Blasi & Colafrancesco 1999;Keshet & Loeb 2010;Enßlin et al. 2011).The limits on the gamma ray emission from ICM (e.g., Jeltema & Profumo 2011;Ackermann et al. 2014Ackermann et al. , 2016) ) suggest that the contribution of secondary particle production through pp collision to the observed emission is sub-dominant (Brunetti et al. 2017;Adam et al. 2021).However, the secondary leptons can contribute to the diffuse radio emission by providing the seed population for the reacceleration (e.g., Brunetti et al. 2017;Pinzke et al. 2017;Nishiwaki et al. 2021).Also, the mini halos in the dense core of cool-core clusters can originate from the injection of secondaries (e.g., Pfrommer & Enßlin 2004;Fujita et al. 2007;ZuHone et al. 2015;Ignesti et al. 2020).In this section, we explore the contribution of the secondary electrons to the diffuse emission found in cluster outskirts.
We calculate the injection of secondary CRe for a given distribution of CRp, using the procedure of Nishiwaki et al. (2021).We assume a single power-law spectrum of CRp, N p ∝ p −αp , with α p = 3.2 to fit the radio synchrotron spectrum with α syn ≈ 1.6.We set p/(m p c) = 1 for the minimum momentum of CRp.The energy density of CRp ϵ CRp is treated as a free parameter.For simplicity, we neglect the re-acceleration of both CRe and CRp.We adopt the one-zone approximation explained in Sect.4.1.
In Fig. 8, we show the energy densities of CRp and the magnetic field required to explain the observed brightness at 49 MHz in the pure hadronic model without reacceleration.We find that ϵ CRp should be as large as ϵ CRp ≈ 10 3 ϵ ICM to explain the observed brightness when we adopt the same magnetic field as Sect.4.1, i.e., B = 0.24 µG.The energy density of secondary CRe is negligible compared to that of CRp.The synchrotron brightness scales as S syn ∝ K e B (δ+1)/2 , where K e ∝ ϵ CRp /(B 2 + B 2 CMB (z)) and δ ≈ α p + 1 are the normalization and the power-law index of secondary CRe spectrum, respectively (e.g., Brunetti & Jones 2014).We find that the minimum value of the non-thermal energy density (ϵ CRp + ϵ B ) is 3 times larger than ϵ ICM , and it appears when the magnetic field is as large as Energy densities of CRp and magnetic field required in the pure hadronic model without reacceleration as a function of the magnetic field.The red and blue lines show ϵCRp and ϵB, respectively.The thick black line shows the sum of those two.The energy are normalized by that of the ICM in the simulated box (ϵICM = 1.3 × 10 −13 erg/cm 3 ).The spectral index of CRp is assumed to be αp = 3.2.The vertical dotted line shows the magnetic field calculated with ηB = 0.05 (Sect.3).B = 6 µG.Thus, the classical "pure-hadronic" model, which does not include any reacceleration, is incompatible with the observed mega halo.
Note that the above discussion does not exclude the possibility that the diffuse emission is produced by the reacceleration of the secondaries injected through the pp collision.A comprehensive study of this "secondary re-acceleration" model should incorporate the re-acceleration of CRp, which is out of the scope of this work.It is worth noting that turbulence might be significantly damped by the back reaction of the reacceleration of CRp when ϵ CRp ≳ ϵ turb (e.g., Brunetti & Lazarian 2007).Such back reaction is not considered in the derivation of Eq. (2).

B. DEPENDENCE ON THE INITIAL CONDITION
In Sect.4, we assumed that the initial spectrum before the onset of the reacceleration, N e (p, 0), is determined by the radiative and the Coulomb cooling.The spectrum has a peak around p min ∼ 10 as shown in Fig. 5.However, it is also possible that the "initial spectrum" is affected by the turbulent reacceleration working before the time of the snapshot.Recent simulation by Beduzzi et al. (2023) observed that the radio-emitting CRe in the ICM experience multiple episodes of reacceleration.In such a case, the peak of the seed CRe spectrum should appear at a larger momentum p min ∼ 10 2 − 10 3 .In this section, we discuss how the results in Sect.4.2 are modified when we adopt a different initial spectrum.We adopt the same method as Sect.4.2, taking into account the LOS integration.We fix η B = 0.05 also in this section.
To consider the reacceleration before the epoch of the snapshot, we assume the initial spectrum for the calculation with p min = 10 3 .The shape of the initial spectrum is assumed to be the same in all cells.The initial energy density of the CRe follows ϵ CRe (t = 0) ∝ ϵ ICM , as in Sect.4.2.
The difference of the initial spectrum is not important when T dur (calculation time of the FP equation) is very long and the steady state due to the balance between the cooling and reacceleration is achieved in most of the cells.In Sect.4, we find that the steady state is achieved at T dur ≥ 3 Gyr when t acc ≈ 0.5 Gyr.
In this Section, we limit T dur by ≈ 2t eddy , where t eddy is the eddy turnover time measured in each cell.For the cells with 2t eddy ≳ 1 Gyr, we terminate the calculation at 1 Gyr.The mean value of t eddy in the simulated box is ≈ 0.5 Gyr, so T dur = 2t eddy is shorter than 3 Gyr in most of the cells.
With the above conditions, we find that the typical spectral index,α syn ≈ 1.6, can be reproduced with ψ = 0.5.Fig. 9 (right) shows the CRe spectra before and after the calculation.Summing up the CRe spectra over cells included in one beam (∼ 10 3 cells), we find that p min after the calculation does not change much from the initial value, p min = 10 3 .This is consistent with the assumption that p min = 10 3 is caused by the turbulent acceleration prior to the epoch of the snapshot.We calculate the acceleration efficiency using Eq. ( 10) and find a typical value of ⟨η acc ⟩ ≈ 1.7%.
Although we assumed a shorter T dur than that in Sect.4.2, the contribution of the cells with smaller turbulent energy (F turb < 5 × 10 44 erg/s/Mpc 3 ) increases, and almost 60% of the area of emission region can be covered at the LOFAR sensitivity.
If we adopt the initial spectrum with p min ∼ 10 and calculate reacceleration for T dur = 2t eddy , the number of CRe with p ∼ 10 4 at the final state becomes smaller than the p min ∼ 10 3 case and that model cannot explain the observed brightness unless η acc ≳ 10%.

Figure 1 .
Figure 1.Maps of the X-ray weighted gas temperature along the line of sight (left) and of the gas velocity along the line of sight (right) for the simulation snapshot analyzed in this work.The additional white contours show the regions which can be approximately detected with X-observations in the 0.5-2 keV energy band, while the green square shows the location of the box used to extract the turbulent flow properties used in Sec.4.Each image is made of 1024 × 1024 pixels.

Figure 2 .
Figure 2. Left: Histogram of turbulent energy flux per unit volume for solenoidal (blue) and compressible (red) modes in the cluster peripheral region shown in Fig. 1.The vertical axis shows the number of cells in the simulated box.The dotted vertical lines show the mean values for each mode.The total number of cells in the extracted box is N = 10 6 .Right: Cumulative number of cells with the turbulent kinetic flux larger than F turb .The vertical dashed line shows F turb = 10 44 erg/s/Mpc 3

a
We use the values found each simulated cell.

Figure 5 .
Figure5.Left: The thick solid line shows the synchrotron spectrum in a single beam, typically seen in our 10 6 cell calculation.For comparison, the spectrum in the one-zone model is shown with the dashed line.The data points show the typical value of the brightness measured at 1.5-2 Mpc from the cluster center(Botteon et al. 2022).The error bars show the ranges of the brightness in this region including 1σ error.We adopt ψ = 0.3 for the one-zone model, while ψ = 0.6 and for the cell-wise calculations.In both models, we assume T dur = 3 Gyr.The beam spectrum is decomposed into the contributions from the cells with F turb < 5 × 10 44 erg/s/Mpc 3 (blue) and F turb > 5 × 10 44 erg/s/Mpc 3 (red).Right: CRe spectra in the same calculations.Spectra before and after the reacceleration is shown with thin and thick black lines, respectively.The dashed lines show the spectra in the one-zone model.

Figure 7 .
Figure 7. Distribution of spectral indices of the beams observable with the LOFAR sensitivity.The vertical axis shows the fraction of beams with each αsyn.The black dotted line shows the observation of Botteon et al. (2022).In our calculation, the total number of observable beams is ≈ 300 (30% of the area), while the observational data shows the distribution of 836 pixels.For the model parameters, we adopted ηB = 0.05 and ψ = 0.6.The reacceleration is calculated for T dur = 3 Gyr, and the beam size is assumed to be 35".We tested the projection along three different axes (x, y, and z, distinguished by colors).
Figure 8.Energy densities of CRp and magnetic field required in the pure hadronic model without reacceleration as a function of the magnetic field.The red and blue lines show ϵCRp and ϵB, respectively.The thick black line shows the sum of those two.The energy are normalized by that of the ICM in the simulated box (ϵICM = 1.3 × 10 −13 erg/cm 3 ).The spectral index of CRp is assumed to be αp = 3.2.The vertical dotted line shows the magnetic field calculated with ηB = 0.05 (Sect.3).

Table 1 .
Parameters adopted in the FP calculations