Turbulence and Particle Acceleration in Giant Radio Haloes: the Origin of Seed Electrons
About of X-ray-luminous clusters show smooth, Mpc-scale radio emission, known as giant radio haloes. One promising model for radio haloes is Fermi-II acceleration of seed relativistic electrons by compressible turbulence. The origin of these seed electrons has never been fully explored. Here, we integrate the Fokker-Planck equation of the cosmic ray (CR) electron and proton distributions when post-processing cosmological simulations of cluster formation, and confront them with radio surface brightness and spectral data of Coma. For standard assumptions, structure formation shocks lead to a seed electron population which produces too centrally concentrated radio emission. Matching observations requires modifying properties of the CR population (rapid streaming; enhanced CR electron acceleration at shocks) or turbulence (increasing turbulent-to-thermal energy density with radius), but at the expense of fine-tuning. In a parameter study, we find that radio properties are exponentially sensitive to the amplitude of turbulence, which is inconsistent with small scatter in scaling relations. This sensitivity is removed if we relate the acceleration time to the turbulent dissipation time. In this case, turbulence above a threshold value provides a fixed amount of amplification; observations could thus potentially constrain the unknown CR seed population. To obtain sufficient acceleration, the turbulent magneto-hydrodynamics cascade has to terminate by transit time damping on CRs, i.e., thermal particles must be scattered by plasma instabilities. Understanding the small scatter in radio halo scaling relations may provide a rich source of insight on plasma processes in clusters.
keywords:acceleration of particles, cosmic rays, turbulence, gamma-rays: galaxies: clusters, radiation mechanisms: non-thermal, galaxies: clusters: general
About one third of X-ray-luminous clusters show smooth, unpolarised radio emission on Mpc scales, known as giant radio haloes (RHs) (Brunetti & Jones, 2014). They appear only in disturbed, merging clusters and the RH luminosity correlates with the X-ray luminosity (Govoni et al., 2001; Feretti et al., 2012) and the Compton -parameter (Basu, 2012; Planck Collaboration et al., 2013). The RHs show that CR electrons and magnetic fields permeate a large volume fraction of the intra-cluster medium (ICM). The dominant CR source, given the smoothness and enormous extent of RHs, is thought to be structure formation shocks (Miniati et al., 2001a; Pfrommer, 2008). At the same time, plasma processes, the origin of magnetic fields and particle acceleration in a turbulent, high- plasma (in which the thermal pressure predominates the magnetic pressure) like the ICM are not well understood. Radio haloes thus provide an incisive probe of non-thermal processes in the ICM.
There have been two competing models proposed to explain RHs. The radio emitting electrons in the “hadronic model” are produced in inelastic (hadronic) CR proton interactions with protons of the ambient thermal ICM, which generates pions that eventually decay into electrons and positrons, depending of the charge of the initial pion (Dennison, 1980; Blasi & Colafrancesco, 1999; Miniati et al., 2001b; Pfrommer & Enßlin, 2004; Pfrommer et al., 2008; Enßlin et al., 2011). CR protons and heavier nuclei may have been accelerated and injected into the ICM by structure formation shocks, active galactic nuclei and galactic winds. However, the strong bimodality that separates X-ray luminous clusters into radio-active and radio-quite clusters (requiring a fast switch on/off mechanism of the RH emission) and the very extended RH emission at low frequencies in Coma (352 MHz) represent a major challenge to this model class (Brunetti et al., 2012; Zandanel et al., 2014).
The alternative model for RHs is re-energisation of seed suprathermal electrons by Fermi II acceleration when ICM turbulence becomes transsonic during mergers (Schlickeiser et al., 1987; Giovannini et al., 1993; Brunetti et al., 2001, 2004; Brunetti & Lazarian, 2007, 2011; Miniati, 2015). Due to the short radiative cooling time of high-energy relativistic electrons, the cluster synchrotron emission quickly fades away after a merger, which naturally explains the observed bimodality of RHs (see e.g. Donnert et al., 2013; Donnert & Brunetti, 2014).
However, there is a salient piece missing in the turbulent reacceleration model. It relies heavily on the assumption of an abundant, volume-filling population of seed suprathermal electrons; direct Fermi II acceleration from the thermal pool is precluded by strong Coulomb losses (Petrosian & East, 2008; Chernyshov et al., 2012). These seeds are presumed to be either fossil CR electrons (CRes) accelerated by diffusive shock acceleration (DSA) during structure formation (Sarazin, 1999), or secondaries injected by hadronic interaction of CR protons (CRps) with thermal protons (Brunetti & Lazarian, 2011).
While analytic estimates have been made, there has been no ab initio demonstration that structure formation can lead to the required abundance of seed electrons with the correct spatial and spectral characteristics. This is a non-trivial requirement: Coulomb cooling in dense cluster cores is severe, and DSA fossil electrons may not survive. On the other hand, for secondaries to constitute the seed population, the CRp population required in the best-studied case of the Coma cluster must have a very broad and flat (or even slightly inverted) spatial profile (Brunetti et al., 2012), in contrast with the thermal plasma whose energy density declines steeply with radius. In Figure 1 we show that such a distribution is not predicted by cosmological simulations (see also Pinzke & Pfrommer, 2010; Vazza et al., 2014). If CRps are predominantly advected with the cluster plasma, their distribution will be peaked towards the cluster centre and show a similar characteristics as the thermal plasma. As a consequence, the distribution of secondary electrons and the resulting radio synchrotron emission is also peaked since the hadronic reaction is a two-body scattering process. Hence, the simulated emission falls short of the observed extended and flat radio profile of the Coma cluster.
Indeed, arriving at a seed population with the required characteristics is highly constraining, and has the potential to teach us much about the origin of CRps/CRes in clusters. In this work, we use our hydrodynamical zoom simulations of galaxy clusters in a cosmological setting to follow the distribution functions of seed populations for CRps and CRes, and integrate the Fokker-Planck equation of CR transport along Lagrangian particle trajectories. We model diffusive shock acceleration at structure formation shocks, and account for various loss processes of CRs. Utilising new insights from our recent work on DSA generated fossil electrons (Pinzke et al., 2013), we generate the first quantitative calculation of primary and secondary seed electrons.
To compare this to observations, we model second-order Fermi acceleration by CR interactions with magnetised turbulence. However, we assume a simplified and stationary model for magnetic fields and turbulence. We do not account for the time-varying energy density in compressible waves, which are thought to be necessary for the acceleration process (Brunetti & Lazarian, 2007, 2011), as the cluster merger proceeds. So our approach is orthogonal (and complementary) to e.g., simulations of Miniati, 2015 that focus on the time-dependent compressible turbulence while adopting a simplified treatment of CR. Our approach of parametrizing turbulence enables us to vary parameters associated with the spatial profile and the overall amplitude of compressible waves (that can in principle vary depending on the details of a particular cluster merger).
In this paper, we explore how the radio surface brightness profile and spectrum of the best known radio halo, Coma, can be used to constrain the underlying properties of the seed CRs and turbulence. We aim to constrain the normalisation and spatial profile of these two input ingredients in turbulent reacceleration models. The outline of this paper is as follows. In Section 2, we outline the basic physics of turbulent reacceleration of CRs which we use. In Section 3, we use cosmological simulations to generate a seed CR population, and combine it with our parametrized model of turbulence to produce radio surface brightness profiles and spectra of Coma. We find that it is possible to fit the observations using physically motivated modifications of the seed population (rapid streaming; enhanced CR electron acceleration at shocks) or turbulence (increasing the turbulent-to-thermal energy density, with radius), but only at the expense of fine tuning. In Section 4, we explore the reason for this fine-tuning, and seek ways to overcome it. We perform a parameter study on spherically symmetric, static models where we vary properties of the seed population and turbulence. We find exponential sensitivity to the amount of turbulence, which can be eliminated if the turbulent acceleration and dissipation time are linked. In this case, turbulence above a threshold value required to overcome cooling provides a fixed amount of amplification; observations could then potentially constrain the unknown CR seed population. We summarise and conclude in Section 5.
2 Cosmic ray transport
The transport of relativistic electrons and protons across cosmic time into galaxy clusters is a complex problem that depends on the velocity field of the gas (and its thermodynamic properties such as density, temperature, and pressure) as well as non-thermal processes (turbulence, magnetic fields, fossil CRs). We use high resolution galaxy cluster simulations to derive the thermal and fossil CR properties (shock accelerated primary CRes and CRps, as well as secondary CRes produced in p-p collisions, see Pfrommer et al., 2007; Pfrommer, 2008; Pinzke & Pfrommer, 2010; Pinzke et al., 2013).
2.1 Basic equations
As previously noted, secondaries produced by shock accelerated CRp have the wrong spatial profile to explain RH observations. Because they arise from a two body process, they are too centrally concentrated. They also produce gamma-ray emission in excess of Fermi-LAT upper limits (Arlen et al., 2012; Brunetti et al., 2012; Ackermann et al., 2014; Ahnen et al., 2016).
Given a seed population of CRs, we adopt essentially the same set of plasma physics assumptions as the reacceleration model for RHs (Brunetti & Lazarian, 2007, 2011). We solve the isotropic, gyro-phase averaged Fokker-Planck equation (via a Crank-Nicholson scheme) for the time evolution of the CRe distribution in the Lagrangian frame (Brunetti & Lazarian, 2007, 2011):
Here is the one-dimensional distribution in position (suppressed for clarity), momentum and time (which is normalised such that the number density is given by ), is the Lagrangian derivative, is the gas velocity, represents Coulomb (Coul, Gould, 1972) and radiative (rad, Rybicki & Lightman, 1979) losses, respectively,
Here is the dimensionless velocity of CRs, is the Lorentz factor of CRs, is the plasma frequency, is the number density of free electrons, and is the Thomson cross section. The rms magnetic field strength is denoted by and the equivalent field strength of the cosmic-microwave background is given by , where denotes the redshift. In the peripheral cluster regions, where , the CRes loose virtually all their energy by means of inverse Compton emission. is the momentum space diffusion coefficient (see Section 2.2), and denotes the injection rate of primary and secondary electrons in the ICM (see section 3.1). The first term containing the expression represents Fermi-I acceleration and the second term of this form describes adiabatic gains and losses.
During post-processing of our Coma-like cluster simulation, we solve the Fokker-Planck equation over a redshift interval from to 0. The simulated cluster undergoes a major merger over the last 1-2 Gyrs that is thought to inject large turbulent eddies. As is commonly assumed (Brunetti & Lazarian, 2007, 2011; Yan & Lazarian, 2004; Beresnyak et al., 2013) we assume that about one Gyr after core passage the fields have decayed down to the smallest scales , and the radio halo turns on shortly after. We choose this simulation snapshot to analyse. Note that in recent simulations by Miniati (2015), the turbulent reacceleration is strongest around core passage. However, we are not very sensitive to the adopted decay time, since the thermal and CR quantities are very similar a few 100 Myrs before and after , where we have chosen to evaluate the simulations. In all our calculations we assume that turbulent reacceleration efficiently accelerates particles for Myrs (which is roughly the cascade time on which turbulence is damped) and that during this turbulent phase CR streaming and spatial diffusion can be neglected. In Section 4.3, we explore sensitivity to the last assumption.
Thus far, we have ignored CR transport. However, if CRps stream in the ICM, then their spatial profile could potentially flatten sufficiently (Enßlin et al., 2011; Wiener et al., 2013). This scenario is very attractive: it generates seed electrons with the right spatial footprint, and by removing CRps from the core, obeys gamma-ray constraints. Turbulence plays two opposing roles: Alfvénic turbulence damps waves generated by the CR streaming instability (Yan & Lazarian, 2002; Farmer & Goldreich, 2004), thus reducing self-confinement; but compressible fast modes scatter CRs directly. Turbulent damping is still efficient for highly subsonic conditions (Wiener et al., 2013), while we assume compressible fast modes to only provide effective spatial confinement during the periods of transsonic, highly super-Alfvénic () turbulence associated with mergers. Thus, CRs can stream out when the cluster is kinematically quiescent. Furthermore, even Alfvénic streaming timescales are relatively short ( Gyr; Wiener et al., 2013) compared to the timescale on which the CRp population is built up. Based on these findings, we adopt a toy model for our M-streaming scenario in which CR streaming instantaneously produces flat CRp profiles. We assume that CRs cannot stream significantly past perpendicular -fields at the accretion shock, so that the total number of CRs is conserved within the virial radius during the streaming process.
The time evolution of the spectral energy distribution of CRps, , is similarly given by:
where is the streaming velocity in the isotropic transport approximation, is the Alfvén speed, and denotes the injection rate of shock accelerated CRps as a function of momentum and time (see section 3.1). The timescale of hadronic losses that produce pions via CRp collisions with thermal protons of the ICM is given by
where we use the cross-section, , given by the fitting formula in Dermer (1986).
2.2 Turbulent reacceleration
In turbulent reacceleration, particles gain energy via Fermi II acceleration. Wave-particle energy exchange takes place via transit time damping (TTD, Brunetti & Lazarian, 2007, 2011). The TTD resonance requires the wave frequency to obey the resonance condition, , where and are the parallel (projected along the magnetic field) wavenumber of the compressible mode and particle velocity, respectively. This implies that the particle transit time across the confining wave region matches the wave period, (hence the moniker, ’transit time damping’). Note that the CRs’ gyroradius does not enter the resonance condition. Hence the CRs that are in resonance with compressible waves experience Fermi-II acceleration irrespective of the length scale of the perturbation.
However, the resonance changes the component of the particle momentum parallel to the mean magnetic field, which over time leads to increasing anisotropy in the particle distribution that decreases the efficiency of reacceleration with time. As in Brunetti & Lazarian (2011), we assume that there exists a mechanism—such as the firehose instability—that isotropises the CR distribution function at the gyroscale and on the reacceleration time scale, which ensures sustained efficient reacceleration with time.
where averages interaction rates over the CR pitch angle ; for ICM conditions, and
is an energy-averaged wavenumber, are the wavenumbers associated with the outer scale and the cutoff scale respectively, and we have assumed a total energy spectrum (composed of both kinetic and potential energy, where the two are assumed to be in equipartition, Sarkar et al. 1991):
which defines the normalisation (the subscript ’c’ emphasizes that we specialise to compressive modes). Intuitively, we can understand the form of the diffusion coefficient from the fact that for second-order Fermi acceleration, , where is the energy averaged wave-particle interaction rate, and is the typical momentum change during wave-particle scattering. Thus:
Equations (6) and (8) make the important aspects of turbulence clear: the energy in compressive modes , the inner and outer scale ( and ), and the slope of the energy spectrum . All of these can vary spatially and temporally.
We adopt a Kraichnan spectrum () for the fast modes, as seen in simulations (Cho & Lazarian, 2003). This can also be written as:
where is the volumetric injection rate of turbulence at scale (which is assumed to be constant), is the phase speed of waves. Since the Kraichnan spectrum is critical in what follows, it is worthwhile taking a moment to recount the origin of the spectral slope and the cascade rate. In the magnetically dominated regime, two counter-propagating wave packets of scale interact on a wave crossing time , rather than the eddy turnover time . Since , each interaction results in a small velocity change . If these changes behave like a random walk, the cascade time it takes for an eddy to become non-linear ( and cascades to smaller scales occurs with a characteristic velocity , which we can solve to obtain the cascade time:
where we have implicitly assumed isotropy. This implies a dissipation rate , or . Thus, since , this gives a kinetic energy spectrum:
It is important to realise that this Kraichnan scaling only applies at small scales, where turbulence is magneto-hydrodynamical. At large scales, turbulent motions are subsonic () but super-Alfvénic (). Thus, except in the case where motions are transsonic and weak shocks become important, motions are fundamentally hydrodynamic and turbulence follows a Kolmogorov () spectrum. However, while the turbulent energy density decreases at small scales, the magnetic energy density (which is dominated by the large scale mean field) is scale-independent. Thus, at some scale , where , the magnetic field becomes dynamically important, and turbulence transitions to the MHD regime. The cutoff scale in equation (8) can be computed by setting the turbulent cascade rate for fast modes to the transit time damping rate on thermal electrons. This yields (Brunetti & Lazarian, 2007; Miniati, 2015):
where . If there is an unimpeded cascade from the injection to the cutoff scale, and . However, in our case the hydrodynamic (Kolmogorov) cascade transitions to the MHD cascade at the Alfvén scale so that and and hence, we get111Note that this differs from previous work, which assumes Kraichnan turbulence from the injection scale onward and effectively adopts wave number and compressible velocity at the outer scale for this estimate of .
This gives kpc in the ICM. This constitutes an effective mean free path for CRs, unless plasma instabilities can mediate interactions between turbulence and particles on smaller scales (Brunetti & Lazarian, 2011), a possibility we discuss in Section 4.3. Another possibility is that compressible modes dissipate in weak shocks, resulting in Burgers’ turbulence (Kowal & Lazarian, 2010; Porter et al., 2015; Miniati, 2015). If Burgers turbulence dominates, then particle acceleration rates are too slow in the face of cooling processes to explain radio haloes (Miniati, 2015), and an alternative model for radio haloes is required.
The spatial profile of injected turbulence depends on the details of the merger such as time during the merger, the impact parameter, the merger mass ratio, and the degree of cluster anisotropy (Miniati, 2015). We parametrize these uncertainties and assume that volumetric injection rate of turbulent energy, , and determine the normalisation by requiring that the turbulent energy in compressible modes , where is the total thermal energy. Given these definitions, one can show that:
What is a typical energy density in compressive fast modes? The total turbulent energy is typically per cent of the thermal energy in a cluster (Vazza et al., 2011a); it rises rapidly during a merger. The compressible component is per cent of the total turbulent energy, and shows more rapid temporal variations compared to the incompressible component (Beresnyak et al., 2013; Miniati, 2015). Note the compressible component in cluster simulations is much larger than in stirring box simulations with similar Mach numbers (Kowal & Lazarian, 2010; Lynn et al., 2014). This is likely due to the compressive nature of turbulent driving (transsonic infall and merger), whereas idealised simulations tend to use incompressible solenoidal driving and allow compressive fluctuations to develop on their own. Overall, we adopt a compressive energy density which is of the thermal energy as a baseline estimate.
The important effects are best summarised in terms of the acceleration rate which is governed by advection in momentum space:
In the last step we have used that for turbulent reacceleration (equation 9). Hence the acceleration time is independent of momentum. This should be compared against the lifetime of turbulence, , and the cooling time . In Table 1, we show both the thermal quantities and the timescales for CR cooling and (re)acceleration for three different spatial regions of the RH. The densities (Briel et al., 1992) and temperatures (Bonamente et al., 2009; Arnaud et al., 2001) are derived from X-ray observations. To calculate synchrotron cooling times, we use -fields derived from Faraday rotation measurements (Bonafede et al., 2010). To calculate the acceleration time, we need to assume an outer scale. We assume an injection scale , where kpc, which were assumptions adopted in previous work (Subramanian et al., 2006; Brunetti & Lazarian, 2007, 2011). This length scale corresponds to an eddy turnover time on the outer scale of Gyr if , as is characteristic of a merger. Note that hydrodynamical simulations of clusters have sometimes found larger , in some cases comparable to the size of the cluster (e.g. in Vazza et al., 2011a; Miniati, 2015). This choice is degenerate with ; in Section 4.3 we argue that is a more appropriate choice, but also find that when the decay time is appropriately scaled, we are relatively insensitive to the choice of . We present for 3 different models, which we discuss in the next section. The reacceleration timescale is similar between our three models, where the difference comes from turbulent profile parametrized with . This implies that even small differences in the turbulent profile could impact the seed CRs significantly. Finally, we adopt an duration of acceleration Myr, in line with previous assumptions (Brunetti & Lazarian, 2007), roughly corresponding to the turbulent decay time.
|[ g cm]||3.0||1.6||0.15|
(1) Median quantities from our simulated post-merging cluster g72a derived during last 300 Myrs in time.
(2) Radius of the giant radio halo in Coma where .
(3) Fermi-II reacceleration for both electrons and protons at all energies.
(4) Inverse Compton and synchrotron cooling for electrons.
(5) Catastrophic losses for protons.
(6) Coulomb cooling for electrons (protons factor smaller).
3 Cosmological simulations
In this section, we solve the Fokker-Planck equation for CR transport on Lagrangian particle trajectories through cosmic history. We then use our parametrized model of turbulence to apply turbulent re-acceleration, and compare radio halo profiles against observations of the Coma cluster. We follow the philosophy of adopting assumptions roughly in line with previous successful models, which are now confronted with more accurate calculations of the seed CR population, and find the minimal modifications required to fit observations. We re-examine these choices in Section 4.
Since the standard vanilla model requires a CR seed population which is inconsistent with the simulations (Figure 1), some modifications are necessary, either in the seed CR population or in the properties of turbulence. We focus on three scenarios. (i) In model M-primaries, we assume that CR electrons, which have been accelerated by cosmic formation shocks and successively cooled by inverse Compton and synchrotron losses, form a fossil seed population for reacceleration. (ii) In model M-streaming, we account for the outward streaming of central CRps, which produces a flat CR distribution in the ICM and equivalently a flat secondary seed population of CRes for reacceleration. (iii) In model M-turbulence we adopt a spatially flatter turbulent profile than what was adopted before but assume that seed CRps and secondary CRes follow the steep profile that is suggested by structure formation simulations.
3.1 Modelling diffusive shock acceleration
In this paper we focus on our simulated cluster, g72a, which is a massive cluster of mass that experienced a merger about 1 Gyr ago (Dolag et al., 2009). Since the cluster mass, density and temperature profiles are all similar to the well studied Coma cluster (Pfrommer et al., 2007; Pinzke & Pfrommer, 2010), we will compare our calculations to radio and gamma-ray observations of Coma.
We use a simple test-particle model for the CR acceleration and injection, where each shock injects CRs that trace a power-law in momentum,
determined by the normalisation and the spectral index that depends on the adiabatic index and the Mach number of the shock (see also Quilis et al., 1998; Miniati et al., 2001a; Pfrommer et al., 2006). It is given by the ratio of the upstream velocity () and the sound speed (). The CR number density and CR energy density are derived from
where is the kinetic energy of a proton with momentum . We adopt a fit to Monte Carlo simulations of the thermal leakage process that relates the momentum of injected protons () to the thermal energy () of the shocked plasma (Kang & Ryu, 2011):
Here , is the amplitude of the downstream MHD wave turbulence, and is the magnetic field along the shock normal. The physical range of is quite uncertain due to complex plasma interactions. In this paper, we adopt , which – as we will later see – corresponds to a conservative maximum energy acceleration efficiency for protons of . To derive the acceleration efficiency, , we first have to infer the particle injection efficiency, which is the fraction of downstream thermal gas particles which experience diffusive shock acceleration (for details see Pinzke et al., 2013),
The particle injection efficiency is a strong function of that depends on both and . The energy density of CRs that are injected and accelerated at the shock (neglecting the CR back reaction on the shock) is given by
and the CR energy injection and acceleration efficiency is:
The dissipated energy density in the downstream regime, , is given by the difference of the thermal energy densities in the pre- and post-shock regimes, corrected for the adiabatic energy increase due to gas compression.
We limit the acceleration efficiency to by steepening the spectral index of the injected population to so that is always fulfilled. The slope impact via the mean energy per particle, . This procedure conserves energy and is motivated by models of non-linear shock acceleration where a sub-shock with a lower compression ratio (and hence steeper spectral index) forms (e.g., Ellison et al., 2000). Given our assumed , we find that for strong shocks where the spectral slope is steepened by a maximum of per cent in low temperature regimes ( keV), while the steepening is much smaller for high temperature regimes ( keV) that are more relevant for clusters. Since remains fixed, so does the CR number density . Hence we can solve for the renormalised normalisation constant using and Eqn. 21:
where the new distribution function is given by . We set an upper limit on the ratio of accelerated proton-to-dissipated energy in the downstream of strong shocks that varies from , depending on the adopted model (for more details, see section 3).
In our Galaxy, the CRe-to-CRp ratio at a few GeV is . Hence, we adopt this as a fiducial value for the CRe-to-CRp acceleration efficiency in our models M-streaming and M-turbulence (see Pinzke et al., 2013, for more discussion). However, as recent PIC simulations have shown, this is likely very different at weak shocks, with electrons efficiently accelerated at perpendicular shocks (Guo et al., 2014a, b), and ions (and electrons) efficiently accelerated at parallel shocks (Caprioli & Spitkovsky, 2014; Park et al., 2015). Thus, depending on magnetic geometry, could be either larger or smaller. Some observations of radio relics suggest high values of , due to the absence of gamma-ray emission, which probes the CRp population (Vazza & Brüggen, 2014). This suggests primary CRes as a viable alternative scenario to secondary CRes as seeds for the giant RHs. In our M-primaries scenario, the injected distribution of CRes is derived in the same way as for the CRps. Once they have been accelerated to relativistic energies, injected electrons and protons are indistinguishable. We therefore assume that CRp and CRe have the same distribution function , with a different normalisation (due to differing acceleration efficiencies) (which is viable for primarily perpendicular shocks, Guo et al., 2014a). Given the unknown magnetic geometry at cluster shocks, we investigate the consequences of this additional degree of freedom.
3.2 Radio emission profile
In Figure 2, we show radial profiles for the radio emission in all three scenarios in which the seeds undergo Fermi-II reacceleration in turbulent fields that are shaped such as to reproduce the Coma RH profile at 352 MHz. Adopting our parametrization for the volumetric injection rate of turbulent energy, , we find that to fit the observations, we require for M-turbulence, for M-streaming, and for M-primaries. As a result, the ratio of turbulent-to-thermal energy densities slightly increase with radius as shown in Figure 3. The fine-tuning of these exponents is somewhat problematic, as we discuss in Section 4. After turbulent reacceleration, the volume-weighted, relative CRp energy density and relative CRp number density inside the RH for M-turbulence (M-streaming), are found to be 2 (3) per cent and (), respectively.
Figure 2 demonstrates that the modelled radio profiles without turbulent reacceleration are too steep. In the bottom right panel of Figure 2 (labelled with Brunetti et al. 2012), we show that our simulated profiles of reaccelerated CRs, which only take advective CR transport into account, i.e. they neglect CR streaming or a flatter turbulent profile, produce radio profiles that are too steep. Indeed, even using the assumptions of previous work, where complete freedom in the seed population was allowed, it is not possible to reproduce observations in both frequencies in any model222Note that in previous work on the Coma cluster, was adopted which approximately corresponds to (Brunetti et al., 2012) and together with the different distributions for seed CRes constitute the main differences compared to our work. – with or without turbulent reacceleration. Decreasing the acceleration efficiency with radius does not change this conclusion much because of the weak radial dependence of . This signals that the problem is generic and requires either additional modifications to the plasma physics of acceleration or a better understanding of potential observational systematics. In addition there are differences in the simulated density and temperature profiles in comparison to the observed profile in Coma that impact the CR abundance as well as cooling and reacceleration.
In Figure 4 we show how the turbulent reacceleration timescales in our three models scale with radius. As expected, the M-primaries model with has the flattest profile with , where the small dip at large radius driven by the decrease in thermal energy density. The M-turbulence model has the flattest turbulent profile parametrized by the smaller which explains the steepest profile. Note that a fixed reacceleration timescale is required to explain the observations at each radius and for each model (see also Table 1). Since , these two parameters are degenerate and can be traded off one another.
In principle, reacceleration via TTD leads to spectral steepening with particle energy due to the inefficiency of the acceleration process to counter the stronger cooling losses with increasing energy. Since synchrotron emission peaks at frequency , this translates into a spectral steepening of the radio spectrum (see the left panel of Figure 5 where the continuous injection of secondary CRes is absent). A given radio window samples higher energy electrons for a decreasing field strength in the cluster outskirts. Hence, the spectral steepening with energy should translate into a radial spectral steepening (Brunetti et al., 2012). However, because of the weak dependence of the electron Lorentz factor on emission frequency (), this effect is only visible in our simulations for GHz. Most importantly, our simulated fluid elements at a given radius sample a broad distribution of shock history, density and temperature, which implies very similar synchrotron brightness profiles at MHz and 1.4 GHz. The discrepancy of the observed and simulated 1.4 GHz profiles could instead be due to systematic flux calibration error in single dish observations. These could arise, for instance, due to errors in point source subtraction. Interestingly, we can match the 1.4 GHz data if we reduce the zero point by adding 0.1 of the central flux to every data point; this flattens the outer profile333Lawrence Rudnick, private communication.. Alternatively, this may point to weaknesses in the theoretical modelling of the particle acceleration process and may require a stronger cutoff in the particle energy spectrum.
3.3 Radio spectrum
In Figure 5 we show that our three models that include Fermi-II reacceleration can individually reproduce the convexly curved total radio spectrum found in the Coma cluster. Seed CRs in M-streaming and M-turbulence that do not experience turbulent reacceleration have a power-law spectrum in disagreement with observations. In order to match both the spatial and spectral profiles in Coma, we adopt an acceleration efficiency for the strongest shocks in our three models M-primaries, M-streaming, and M-turbulence to , , and , respectively. Following the Mach number ()-dependence of the acceleration efficiency suggested in Pinzke et al. (2013), the efficiency in weak shocks () that dominates the CR distribution function, has an acceleration efficiency for protons , and for electrons .
Interestingly, we find that the radio luminosity from clusters in the OFF-state (DSA only) and ON-state (DSA and reacceleration) differ by about a factor 10-20 in all our three models. This means that the secondary CRes are dominated by the reaccelerated fossil CRes and not from the CRes produced by reaccelerated CRps. However, for high frequencies ( GHz) where synchrotron cooling is more efficient than reacceleration, the emission is dominated by the CRes produced in the continuous injection of electrons from reaccelerated CRps. It is also worth mentioning that the radio emission from secondary CRes are smoothly distributed around the cluster because of the continuous injection, hence it is not dominated by outliers.
However, for M-primaries, the primary CRes that generate most of the radio emission from the cluster in the OFF-state are dominated by only a small fraction of the CRes. These electrons are injected very recently and have not had time to cool yet. Hence we expect there to be a large variance in the OFF-state of different simulated clusters. As mentioned in section 4, combining radio observations with gamma-ray limits allows us to put a lower limit to . If is smaller than in our adopted models (where we assume ), then the efficiency of DSA has to be larger than for the secondary CRes to reproduce the radio observations. However, since the turbulent reacceleration acts on both the secondary CRes and the CRps, while only affects the CRps, M-streaming and M-turbulence would produce too much gamma-rays. Hence we conclude that if all other parameters are kept fixed. Although, we caution the reader to take this limit too stringent because of the uncertainty in and that impact for a fixed . This parameter space needs to be explored further in future work in order to put more stringent limits on the level of turbulence in clusters using radio and gamma-ray observations in combination with turbulent reaccelerated CRs.
3.4 Gamma rays
The gamma-ray emission from CRps that produce decaying neutral pions could be substantial if the CRps are reaccelerated efficiently enough, hence it is interesting to estimate this emission for our models and compare to upper limits. We follow the formalism outlined in Blasi & Colafrancesco (1999) (and references therein) and calculate the gamma-ray emission numerically for our three models. We predict the gamma-ray emission from M-turbulence (M-streaming) with . The fluxes from these models are slightly larger than in Brunetti et al. (2012), where the differences comes from our steeper CRp profiles in addition to the simulation based formalism we rely on that accounts for both Coulomb and hadronic losses during the build up of the CR distribution in contrast to the scaling relations adopted in their paper. Interestingly the gamma-ray flux from both our scenarios are just below recent Fermi-LAT limits derived from a gamma-ray profile similar to M-turbulence 444Fabio Zandanel, private communication; see also Zandanel & Ando (2014); Ackermann et al. (2014); this will be probed in the next few years by Fermi-LAT. where . The spectral index of the CRp distribution is relatively steep () for the CRp energies GeV that are relevant for the injection of radio-emitting secondary CRes. The steep spectrum is ultimately a consequence of the shock history of the simulated cluster, with a weak dependence on our test particle model for Fermi-I acceleration (Pinzke et al., 2013), where we steepen the spectral index to avoid acceleration efficiencies above .
4 Parameter Space Exploration and Overcoming Fine Tuning
In this paper we rely on several critical parameters describing relatively unknown non-thermal physics in the ICM. Here we develop a simplified framework for our reacceleration model of secondary electrons. We will explore how radio emission depends on the parameters describing the spatial profile of CR protons and turbulence. Our fiducial model is meant to be compared against the Coma cluster.
As we have seen before, the most uncertain aspects of radio halo emission models are the profile of compressible turbulence (which determines the amount of Fermi II acceleration) and the distribution of pre-existing CRs. Hence we vary parameters describing the amount of energy contained in turbulence, (defined by ), the spatial profile of turbulence as parametrized by (defined by , where is the injection rate of turbulence), as well as the spatial CR profile that we will parametrize by (see below). We hold fixed thermal plasma properties (temperature and density profiles), -field profiles, total CR energy content, and the turbulent outer scale (corresponding to a wavenumber ). The CR energy content is suggested by our simulations (observations only give an upper bound; Arlen et al. 2012). We focus on the uncertain CR distribution rather than the overall normalisation, since the impact of the latter (an overall linear scaling) is clear.
In order to quickly explore this parameter space, we solve the Fokker-Planck equation in static spherical shells for injection, reacceleration, and losses of the CRs, i.e., we neglect Lagrangian evolution during re-acceleration. This ignores the effect of adiabatic compressive heating of the CRs, though this is generally subdominant (e.g., see Figure 7 of Miniati 2015). Once the CRs have been reaccelerated for , we calculate the resulting radio emission numerically using the formalism outlined in Rybicki & Lightman (1979) and compare the emission profiles and spectra as we vary one parameter at a time relative to our fiducial model.
with . The virial and core radii of Coma are given by (Reiprich & Böhringer, 2002) and , respectively. In accordance with Faraday rotation measure measurements, we assume , where G and (Bonafede et al., 2010).
The bulk of the CRps are injected by relatively low Mach number shocks and parametrized by , where in our simulations. The CRps approximately trace the thermal gas with (Pinzke & Pfrommer, 2010; Vazza et al., 2016), where the normalisation is fixed by the injection rate of CR energy in the last 650 Myrs. Our simulations show that the CR energy approximately amounts to 0.03 per cent of the thermal energy inside the virial radius, i.e. . The spectrum of the initial CRp distribution is determined by the steady state between injection and cooling,
where we fix the normalisation by requiring . Note that the injected CR energy in the last 650 Myr is smaller by about a factor of 10 in comparison to the cumulative CR energy injected over the entire cosmological history of the cluster. However, since the injected CR energy rate averaged over the formation time of the cluster is similar to that during a merger, the CRes injected during the merger are especially important for the radio emission above 1 GHz since these CRs have not yet had time to cool. Similarly, the initial CRe distribution is given by the steady state between cooling (Coulomb, inverse Compton, and synchrotron) and injection of secondary CRes from . The diffusion constant is calculated for each radial bin. All parameters and assumptions are similar to what is used for our simulated cluster (see section 3). Our fiducial model assumes , , , and where kpc. We also assume a fixed acceleration time of .
Figure 6 shows the impact of turbulence and the spatial distribution of CRs on the radio emission.
Impact of overall level of turbulence (). From the top panels of Figure 6, we see that as increases, there are 3 important changes: an exponential increase in radio luminosity, a flattening of the radio surface brightness profile, and an increase in spectral curvature. We discuss these in turn.
As we show in Section 4.3, the exponential increase in radio surface brightness is easily understood from equation 30, since for a power-law initial distribution function , then . This exponential sensitivity is somewhat modified by cooling (which results in a non power-law spectrum; also, the shorter the acceleration time, the larger the pool of seed electrons available which would otherwise cool away), but overall is a good approximation.
An increase in flattens the surface brightness profile, since high acceleration efficiency leads to larger amplification in the cluster outskirts (where cooling is less important) than the centre. In particular, in the cluster outskirts, the reduced impact of Coulomb cooling implies that there is a larger pool of low-energy electrons available for reacceleration (see timescales in Table 1).
Larger levels of turbulence also increase spectral curvature, which might seem puzzling. It can be understood as follows. The pre-acceleration electron distribution function results from the competition between hadronic injection and cooling, which results in a quasi-steady state for the non-thermal (secondary) electrons. At low momenta, when the Coulomb cooling time is short, , at high momenta, when inverse Compton/synchrotron cooling dominates, . In between, there is a quasi-adiabatic regime where electrons accumulate (for more details, and an analytic self-similar solution, see Sarazin 1999; Pinzke et al. 2013). In the absence of cooling, momentum advection in the limit (so that is momentum independent) simply shifts the distribution . When the acceleration efficiency is low, most observable emission corresponds to the power-law tail of the distribution function set by the balance between injection and synchrotron/IC cooling. However, as the acceleration efficiency increases, radio emission starts to probe the ’bump’ around (given by ) where electrons accumulate and the distribution function is curved. This results in a curved emission spectrum. The synchrotron spectrum steepens at the frequency (Brunetti et al., 2001):
where , which increases for shorter .
Impact of spatial profile of turbulence (). From the middle panels of Figure 6, we see that as expected, a flatter profile of the turbulent pressure directly translates into a flatter radio surface brightness profile. Since seed electrons are more concentrated toward the centre (the collisional production of secondaries is more rapid there), and magnetic fields are stronger, concentrating the turbulence toward the cluster centre for fixed results in higher radio luminosities, and slightly more curvature (due to the increased importance of cooling near the centre). Overall, however, the spatial profile of turbulence has a much weaker effect than its overall normalisation.
Impact of spatial profile of seed CRs . The spatial distribution of CRs has an even smaller impact on radio emission. At fixed total CR energy content , concentrating the CRs towards the centre leads to more centrally dominated surface brightness profiles, as expected, and higher radio luminosities (for the same reasons as above: secondaries are more easily produced in the centre, and magnetic fields are stronger). We have also confirmed that radio surface brightness profiles scale linearly with , as expected.
Overall, our results suggest that radio haloes are much more sensitive to the level of turbulence (exponential dependence) rather than CR abundance (linear dependence), and that the spatial distribution of turbulence and CRs, while important, are second-order effects. In our parametrization, the most important controlling variable is . The overall level of turbulence has to be such that (see Table 1), otherwise, too little or too much amplification takes place. For instance, for , changing by a factor of two changes the radio surface brightness by a factor of (see top panels of Figure 6). The required depends only logarithmically on the abundance of seed CRs.
While the requirement of a threshold level of turbulence may explain why radio brightness is bimodal, it also raises a fine-tuning problem: why is the vs. relation in active radio haloes so tight (Brunetti et al., 2009)? Depending on the details of infall or mergers, we would naturally expect fluctuations in , which would translate into large scatter in the – relation. This can be only be understood if the timescale over which acceleration takes place also depends on the properties of turbulence, so that the ratio has relatively little scatter. We address this next.
4.3 Self-limiting turbulent reacceleration
In this section, we explore the physical origin of the high sensitivity of radio emission to turbulence levels (e.g., top right panel of Fig. 6), which thus requires strong fine-tuning to explain observed radio profiles. We will find that by relating the acceleration time to the turbulent decay time, this sensitivity can be eliminated. We also challenge the common assumption of assuming a Kraichnan spectrum, beginning at the outer scale (see Section 2.2). Since turbulence is essentially hydrodynamic at large scales, this is unjustified. Instead, we shall assume that the outer scale of the compressive fast modes is the Alfven scale.
There are 3 important timescales in this problem: the acceleration time , the duration over which turbulence is active and acceleration takes place, , and the cooling time . Only is momentum dependent (and is different for ions and electrons). Thus, the outcome of acceleration depends essentially on two dimensionless numbers, , and .
We have seen that the radio luminosity depends very sensitively on , through the very sensitive dependence on (Figure 6; note from equation (15) that ). Since this raises questions of fine-tuning in to explain the observations, it is worth understanding this property in more detail. To this end, we ignore cooling, which is a good approximation for the CR protons, because hadronic cooling times are long in comparison to the other relevant timescales of the problem.555In principle, CRp can still lose energy via wave heating, at a rate , but we assume CR streaming is suppressed during mergers, when they are spatially confined by scattering. In this case, , and after a time , we have
(where we have used the fact that is momentum independent). For an initial power law distribution function , this momentum increase can be rewritten as a change of normalisation, , where . A slightly more careful derivation by direct solution of the Fokker-Planck equation yields:
with the analytic solution given by
At the CR distribution is exponentially sensitive to , which is the physical reason underlying the extreme sensitivity to in Figure 6.
The natural solution to this puzzle would be a process that couples and . There are two possible limits. Let us suppose that the timescale on which there exists a source of turbulent driving is , which approximately corresponds to the timescale of the merger, or the dynamical time. Let the turbulence decay on a timescale . If , then , which depends on the details of the merger and should result in considerable scatter in in different systems. On the other hand, if , then . Since and are both related to properties of the turbulence, it is conceivable that this would result in much less scatter in .
An obvious candidate for is:
i.e., the eddy turnover time at the outer scale. This is subject to uncertainties about the location of the outer scale ; estimates in the literature range from Mpc. It is also worth remembering that MHD turbulence only applies for , where is the Alfvén scale where . Invoking fast modes, Kraichnan scalings, etc., is only valid below these scales. For , turbulence is basically hydrodynamic.
In the hydrodynamic regime, a standard Hodge-Helmholtz decomposition usually shows that the compressive component of the velocity field is Burgers-like (), while the solenoidal component is Kolmogorov-like (), see e.g., Federrath (2013). The Burgers-like component does not reflect a genuine cascade, but rather the appearance of shocks which directly transfer power from large to small scales. At face value, we should use the Burgers spectrum for compressible modes. However, as already found by previous authors (Miniati, 2015; Brunetti & Lazarian, 2016), and as we shall discuss, this does not produce significant particle acceleration. If this is the correct spectrum, then the paradigm of turbulent particle acceleration is simply flawed, and some other mechanism is necessary to explain radio halos. Alternatively, it is well-known that due to mode-mode coupling, solenoidal modes can give rise to compressive modes, and vice-versa (Kida & Orszag, 1990; Cho & Lazarian, 2003; Kritsuk et al., 2007), even for subsonic turbulence, since pressure fluctuations of order arise. Indeed, note that many numerical studies (such as Cho & Lazarian (2003); Kritsuk et al. (2007)) only have solenoidal driving on the outer scale, but then are also able to study the compressive modes that develop. In hydrodynamic turbulence, the energy in compressible modes which develop in this way scales as , for ; this coupling is strongest at the Alfven scale (Cho & Lazarian, 2003). In MHD turbulence, mode-mode coupling is weak below the Alfven scale and the Alfven, fast and slow modes proceed as separate cascades.
We thus assume that some fraction of the Kolmogorov-like hydrodynamic turbulence which cascades down to the Alfven scale ends up as compressive fast modes with a Kraichnan spectrum, as seen in simulations (Cho & Lazarian, 2003). Henceforth, we can consider the outer scale of the fast modes to be . If some fraction of the turbulent energy density at this scale is in compressible modes, then the energy density at the outer scale is . The turbulent reacceleration time is:
This can be related to the eddy turnover time at the outer scale from to:
where we have used and . The decay time at the Alfven scale is comparable (and could be larger, if is smaller than we have assumed) to the eddy turnover time at the outer scale. Thus, it is appropriate to consider the latter as the decay time for the fast modes, rather than the cascade time at the driving scale. The reason for these comparable timescales despite the disparate length scales is that large scale hydrodynamic modes cascade on a single eddy turnover time, whereas small scale sub-Alfvenic MHD modes are in the weak turbulence limit, and require multiple wave-wave interactions for a mode to cascade.
Remarkably, this expression is independent of properties of the turbulence such as , , or , and depends only on properties of the plasma (). Ultimately, this arises because the timescale on which a fast mode wave transfers energy due to wave particle interactions in second order Fermi acceleration, , is closely related to the timescale on which it cascades due to wave-wave interactions, . This implies , where is the outer scale on which the fast mode cascade begins, and is the characteristic wavelength for wave-particle interactions. For transit time damping and an outer scale of , we have . We thus have . More careful consideration of dimensionless factors boosts this estimate by an order of magnitude to give equation (35).
On the other hand, equation (35) points toward a pessimistic scenario where turbulent reacceleration with TTD on thermal particles is never effective. A key reason is that we assume Kraichnan turbulence only applies for . There is then insufficient separation of scales: the cutoff scale . Although there is a fair large separation of scales between the outer driving scale and (a factor of for ; we have for our fiducial assumptions), Kolmogorov turbulence, with its steeper spectrum, has more energy at large scales ( for Kolmogorov and Kraichnan turbulence, respectively). This implies that the energy-weighted scales at which wave-particle interaction take place are larger in Kolmogorov turbulence, and thus that the wave-particle interaction rate is lower. If (as is frequently seen) we instead assume that Kraichnan turbulence begins at the outer scale , with characteristic decay time , then we obtain a more palatable result:
This arises because the decay time is now times longer, and the acceleration time is now times shorter, boosting by . However, as we have argued, turbulence in the hydrodynamic regime