High $\beta$ Effects on Cosmic Ray Streaming in Galaxy Clusters

TL;DR

This paper investigates how high plasma beta values in galaxy clusters enhance MHD wave damping, leading to faster cosmic ray streaming and reduced radio halo brightness, revising previous models of CR transport.

Contribution

It introduces a revised treatment of MHD wave damping in high beta plasmas, showing increased CR streaming speeds and lower radio halo surface brightness predictions.

Findings

High beta enhances MHD wave damping rates by a factor of beta^{1/2}.

Faster CR streaming results in decreased radio halo brightness.

Re-examination of previous simulations highlights the importance of high beta effects.

Abstract

Diffuse, extended radio emission in galaxy clusters, commonly referred to as radio halos, indicate the presence of high energy cosmic ray (CR) electrons and cluster-wide magnetic fields. We can predict from theory the expected surface brightness of a radio halo, given magnetic field and CR density profiles. Previous studies have shown that the nature of CR transport can radically effect the expected radio halo emission from clusters (Wiener et al. 2013). Reasonable levels of magnetohydrodynamic (MHD) wave damping can lead to significant CR streaming speeds. But a careful treatment of MHD waves in a high $\beta$ plasma, as expected in cluster environments, reveals damping rates may be enhanced by a factor of $\beta^{1/2}$. This leads to faster CR streaming and lower surface brightnesses than without this effect. In this work we re-examine the simplified, 1D Coma cluster simulations (with radial magnetic fields) of Wiener et al. (2013) and discuss observable consequences of this high $\beta$ damping. Future work is required to study this effect in more realistic simulations.

Figures

Equations

$${\Gamma }_{cr}({\mathbf{k}}_{\parallel })={\omega }_{cp}\frac{n_{CR}(>{\gamma }_R)}{n_i}\left(\frac{v_D({\gamma }_R)}{v_A}-1\right)$$

$$\tan {\theta }_{min}=\left(\frac{{\lambda }_{\parallel }}{{\lambda }_{\perp }}\right)\sim {\left(\frac{r_L}{L_{\text{MHD}}}\right)}^{1/4}\sim {\left(\frac{r_Lϵ}{v_A^3}\right)}^{1/4},$$

$${\lambda }_{\perp }\sim r_L^{3/4}{L}_{\text{MHD}}^{1/4}\sim r_L^{3/4}(v_A^3/ϵ{)}^{1/4},$$

$$v_{{\lambda }_{\perp }}\sim v_A{\left(\frac{{\lambda }_{\perp }}{L_{\text{MHD}}}\right)}^{1/3}\sim v_A{\left(\frac{r_L}{L_{\text{MHD}}}\right)}^{1/4}.$$

$${\Gamma }_{\text{damp}}\sim \frac{v_{{\lambda }_{\perp }}}{{\lambda }_{\perp }}\sim \frac{{ϵ}^{1/3}}{{\lambda }_{\perp }^{2/3}}\sim {\left(\frac{ϵ}{r_L{v}_A}\right)}^{1/2},$$

$$\beta =\frac{P_g}{P_B}=\frac{\rho k_{B}T/\mu m_p}{B^2/8\pi }=\frac{v_i^2}{v_A^2},$$

$$\nu \equiv \frac{\omega }{k_0{v}_A},l\equiv \frac{k\cos \theta }{k_0},k_0\equiv \frac{{\omega }_{ci}{v}_A}{v_i^2},$$

$$\nu =l-\frac{i\alpha }{4l}\pm \frac{1}{4l}{\left(l^6-{\alpha }^2+\frac{4i{\alpha }^3}{l^2}\right)}^{1/2}$$

$${\Gamma }_{\text{damp}}=-k_0{v}_A\text{Im}(\nu )\approx k_0{v}_A\frac{\alpha }{4l}$$

$${\Gamma }_{\text{damp}}\approx \frac{\sqrt{\pi }}{4}{k}_0{v}_A{\beta }^{1/2}l{\tan }^2\theta \approx \frac{\sqrt{\pi }}{4}{\beta }^{1/2}{r}_L^{-1}{v}_A{\left(\frac{r_Lϵ}{v_A^3}\right)}^{1/2}$$

$$\Rightarrow {\Gamma }_{\text{damp}}\approx \frac{\sqrt{\pi }}{4}{\beta }^{1/2}{\left(\frac{ϵ}{r_L{v}_A}\right)}^{1/2},$$

$${\Delta }_i\equiv \frac{P_{\perp ,i}-P_{\parallel ,i}}{P_i}=\frac{2.9}{{\nu }_{ii}}\left(\mathbf{bb}:\nabla \mathbf{u}-\frac{1}{3}\nabla \cdot \mathbf{u}\right).$$

$${\Delta }_i>\frac{1}{{\beta }_i},\text{mirror instability}$$

$${\Delta }_i<-\frac{2}{{\beta }_i},\text{firehose instability}$$

$${\Delta }_i=\frac{2.9}{{\nu }_i}S=\frac{2\xi }{{\beta }_i}\to {\nu }_i=\frac{1.45{\beta }_i}{\xi }S,$$

$${\nu }_i\sim \frac{1.45{\beta }_i}{\xi }\frac{U}{L}=\frac{1.45{\beta }_i}{\xi }\frac{U_0}{L_0}{\left(\frac{L_0}{L}\right)}^{2/3}$$

$$\frac{U_d}{L_d}=\frac{U_0}{L_0}{\left(\frac{L_0}{L_d}\right)}^{2/3}\sim \frac{v_i^2}{{\nu }_i{L}_d^2}$$

$$\Rightarrow {\left(\frac{L_d}{L_0}\right)}^{4/3}\sim \frac{v_i^2}{{\nu }_i{U}_0{L}_0}$$

$$\frac{L_d}{L_0}={\left(\frac{v_i^2}{{\nu }_i{U}_0{L}_0}\right)}^{3/4}$$

$${\nu }_i\sim \frac{1.45{\beta }_i}{\xi }\frac{U_0}{L_0}{\left(\frac{{\nu }_i{U}_0{L}_0}{v_i^2}\right)}^{1/2}$$

$$\Rightarrow {\nu }_i\sim {\left(\frac{1.45{\beta }_i}{\xi }\right)}^2\frac{U_0^3}{L_0{v}_i^2}$$

$${\lambda }_i\sim 0.48{\xi }^2\frac{L_0}{{\beta }_i^2}{\left(\frac{v_i}{U_0}\right)}^3\sim 0.48\frac{L_0}{{\beta }_i^2{\mathcal{M}}_0^3}$$

$$L_d\sim L_0{\left(\frac{v_i^2}{{\beta }_i^2{U}_0^4/v_i^2}\right)}^{3/4}\sim L_0{\mathcal{M}}_0^{-3}{\beta }_i^{-3/2}$$

$$U_d=U_0{\left(\frac{L_d}{L_0}\right)}^{1/3}\sim U_0{\mathcal{M}}_0^{-1}{\beta }_i^{-1/2}=v_A.$$

$${\lambda }_i\sim L_0{\beta }_i^{-2}{\mathcal{M}}_0^{-3}\sim L_d{\beta }_i^{-1/2}\sim 14\text{kpc}\left(\frac{L_{\text{MHD}}}{100\text{kpc}}\right){\left(\frac{{\beta }_i}{50}\right)}^{-1/2}$$

$${\lambda }_{ii}=v_i{t}_i=\sqrt{\frac{kT}{m_p}}\frac{\sqrt{m_p}(kT{)}^{3/2}}{4\sqrt{\pi }\ln \Lambda e^4{n}_i}\approx 5\times 10^{20}\text{cm}\frac{T_7^2}{n_{i,-3}}$$

$$t_{\text{st}}\sim D/v_{\text{st}},$$

$$\kappa \sim L{v}_{\text{st}}.$$

$$t_{\text{diff}}\sim \frac{D^2}{\kappa }\sim \frac{D}{L}\frac{D}{v_{\text{st}}}\sim \frac{D}{L}{t}_{\text{st}}.$$

$$\frac{1}{2}\rho v_{\text{t}}^2=\frac{B^2}{8\pi }=\frac{1}{2}\rho v_{\text{A}}^2$$

Full text

††pagerange: 000–000††pubyear: 0000

1 Introduction

Galaxy clusters are the largest gravitationally bound objects in the universe, and are host to processes described by every area of physics. For some clusters, this includes particle acceleration to very high energies, generating cosmic rays (CRs). Cosmic ray electrons (CRe) at these energies will emit synchrotron emission in the presence of magnetic fields which is observable as radio emission.

Galaxy clusters can exhibit radio emission of a few morphological types. We focus here on radio haloes, which refers to extended, diffuse radio emission. Studies have shown that where radio haloes are detected, the total radio luminosity of the halo correlates well with total X-ray luminosity of the host cluster. But many clusters do not have radio haloes at all. Thus there is a peculiar bimodality in the presence of radio haloes in galaxy clusters that must be explained (Brunetti et al. (2007), Enßlin et al. (2011), Basu (2012); Sommer & Basu (2014)).

Related to this problem is the source of the CRe which produce the emission. CRe at the required energies ($\sim$1-10 GeV) have short cooling times ( \sim$$< 100 Myr) for typical cluster magnetic field strengths ($\sim$3 $\mu$G) and densities ($\sim 10^{-3}\ \textrm{cm}^{-3}$). The presence of a radio halo in a given cluster therefore indicates relatively recent injection of CRe. The ‘hadronic model’ claims that this is provided by hadronic collisions of CR protons (CRp) with ambient thermal nuclei that result in pions, which then decay into secondary CRe. A rival model, ‘turbulent reacceleration’, hypothesizes that low-energy CRe undergo Fermi-II reacceleration by gas motions during mergers (Brunetti & Lazarian, 2007). However, the abundance and distribution of CRp is still an important ingredient in this model, as secondary CRe produced during hadronic collisions provide seeds for reacceleration (Brunetti & Lazarian, 2011; Pinzke et al., 2017).111Acceleration from the thermal pool is precluded by strong Coulomb losses (Petrosian & East, 2008). An excellent summary of the state of the field can be found in Brunetti & Jones (2014).

An issue with the hadronic model is that high energy CRp have much longer ($\sim$10-100 Gyr) cooling times. If a cluster undergoes some merger or structure formation process that may be responsible for accelerating CRs, then we may naively expect the resulting radio halo to last for as long as or longer than the age of the universe. If this were the case, we would struggle to explain clusters that don’t exhibit radio halo emission, since every cluster is expected to accelerate CRp at some point in its history through such a mechanism. A lack of hadronic emission implies either that the assumed acceleration efficiencies are too high or that the CRp are quickly diluted by transport effects. We investigate this second possibility. Note that CRp are frequently invoked in turbulent reacceleration models as well, to generate secondary seed electrons for reacceleration. Observations are best fit by a flat CRp profile, which can be explained by CR streaming (Pinzke et al., 2017).

Building on the work of Enßlin et al. (2011), Wiener et al. (2013) attempted to resolve this issue by suggesting that, for favorable magnetic field orientations, CRp can travel at high speeds away from cluster centers, quickly reducing the CRp density and thus the expected radio luminosity on appropriate time scales. Bulk CR transport is limited by the streaming instability, which amplifies hydromagnetic waves traveling in the same direction as the mean cosmic ray velocity if the cosmic rays are anisotropic in the frame of the wave (Kulsrud & Pearce, 1969). However, if other plasma processes damp the wave then a larger streaming speed is required to excite the instability, resulting in faster cosmic ray transport. In Wiener et al. (2013) it was shown that under galaxy cluster plasma conditions the strongest damping mechanism is the turbulent damping process proposed by Yan & Lazarian (2002); Farmer & Goldreich (2004).

Galaxy cluster plasmas are characterized by large ratios of thermal to magnetic pressure, a parameter usually denoted as $\beta$. In this work we show that, in high $\beta$ environments, wave damping is actually stronger than the rate quoted in Wiener et al. (2013). Therefore, CRs can travel even faster than previously thought and radio halos turn off faster than we previously estimated for the same density and magnetic field profile222As we discuss in §2.3, we use a very simple 1D halo model with purely radial magnetic fields for comparison with Wiener et al. (2013). More sophisticated simulations are needed in future work to examine the effects of this damping in real clusters.. While this may explain the bimodality of radio haloes in the hadronic model, it cannot reconcile other issues with the hadronic model. Among others, gamma-ray non-detections in the Coma cluster put strict limits on the amount of CRp present (Brunetti et al. (2012); Zandanel & Ando (2014); Ackermann et al. (2016); Brunetti et al. (2017)). As such, this work should not be seen as an attempt to strengthen the hadronic model of radio haloes so much as a more accurate theory of CR transport in clusters in general, which is still an important ingredient in the reacceleration model (Pinzke et al. (2017)).

Also, although the time evolution of radio halos is the main focus of the paper, we note that faster streaming may also help to explain the puzzling lack of diffuse $\gamma$-rays from galaxy cluster cores (Pfrommer (2008); Pinzke & Pfrommer (2010); Enßlin et al. (2011); Ahnen et al. (2016)). It may also be germane to understanding radio mini-halo luminosities (Jacob & Pfrommer, 2017). While the hadronic model is problematic for understanding giant radio halos, it remains the chief contender for understanding kinematically quiescent radio mini-halos (Zandanel et al., 2014), where the CRp could be sourced by a central AGN. Note that modulo loss processes, the contribution of the AGN jet to CRe could be non-negligible as well.

In §2 we review the theory of MHD wave damping in collisionless, high $\beta$ plasmas and show how it modifies previously estimated turbulent damping rates. An important part of this calculation is an estimate of the effective collisionality of the plasma under the assumption that microturbulence is present due to pressure anisotropies (Schekochihin & Cowley, 2007); we argue that the level of microturbulence is low enough that Alfvén waves at the characteristic wavelengths excited by cosmic ray steaming instabilities are nevertheless collisionlessly damped. We also discuss the possible effects of the magnetic field topology, which is not taken into account in our simulations. In §3 we describe the simulation setup, which is similar to the setup we used previously for the Coma cluster (Wiener et al., 2013). Sections §4 and §5 gives the results and conclusion, respectively.

2 MHD Waves in Turbulent, High Beta Plasma

2.1 Wave damping

As mentioned in §1, the streaming speed as a function of energy is determined by balancing the growth rate of the streaming instability (Kulsrud & Cesarsky, 1971)

[TABLE]

against the mechanism(s) that damp the wave. In (1), $\omega_{cp}$ is the proton gyrofrequency (for proton cosmic rays), $\gamma_{R}\equiv\omega_{cp}/ck_{\parallel}$ is the minimum Lorentz factor of a resonant cosmic ray, and $v_{D}$ is the energy dependent streaming speed.

Equation (1) is computed assuming a perfectly straight and uniform background magnetic field $\mathbf{B}_{0}$. Under this assumption, the fastest growing waves propagate parallel to $\mathbf{B}_{0}$, with waves of any given wavenumber $k_{\parallel}$ being driven primarily by cosmic rays with gyroradii $r_{L}\sim k_{\parallel}^{-1}$.

If $\mathbf{B}_{0}$ has curvature or perpendicular structure, strict parallel propagation is impossible: there is a minimum angle $\theta_{min}$ which depends on the background field structure. Farmer & Goldreich (2004) evaluated $\theta_{min}$ for Alfvén waves excited by cosmic ray streaming in a background field upon which an anisotropic MHD turbulent cascade is imposed, and found it to be

[TABLE]

where $L_{\textrm{MHD}}$ is the length scale at which the characteristic bulk speed is the Alfvén speed $v_{A}$, and the turbulent energy cascade rate is $\epsilon\equiv v_{A}^{3}/L_{\textrm{MHD}}$. The perpendicular lengthscale $\lambda_{\perp}$ and turbulent velocity $v_{\lambda_{\perp}}$ corresponding to $\theta_{min}$ are

[TABLE]

and

[TABLE]

According to the Farmer & Goldreich (2004) scenario, MHD waves are progressively sheared by the turbulent cascade, eventually transferring their power to the dissipation scale. The corresponding “turbulent damping rate” is therefore on the order of the eddy turnover rate at the perpendicular scale

[TABLE]

where in the last equality we have used eqns. (3) and (4). This was the damping rate used in the ZEUS simulations in Wiener et al. (2013). The Farmer & Goldreich (2004) calculations have recently been generalized by Lazarian (2016) to a variety of scenarios.

A key assumption in deriving eqn. (5) is that the dissipation scale is much shorter than the wavelength of the wave. If this is not the case, the dissipation rate can exceed the shearing rate.

In collisionless, high $\beta$ plasmas, MHD waves of even small obliquity are subject to strong ion Landau damping (Foote & Kulsrud, 1979), a process whereby ions with parallel velocity $v_{\parallel}$ which satisfies the resonance condition $v_{\parallel}=\omega/k_{\parallel}$ absorb energy from the wave due to acceleration by its parallel electric field 333This process is not to be confused with nonlinear Landau damping, in which the thermal ions resonate with the low frequency beat wave generated by the interaction of two higher frequency waves (Kulsrud, 2005). NLLD was first invoked for galaxy cluster plasmas by Loewenstein et al. (1991) and was argued to be subdominant in turbulent galaxy cluster plasmas in Wiener et al. (2013)..

We briefly review the Foote & Kulsrud calculation here. They defined

[TABLE]

and framed their analysis in terms of the rescaled wave frequency $\omega$ and wave number $k$:

[TABLE]

where $\omega_{ci}$ is the thermal ion cyclotron frequency. They also introduce a parameter $\alpha=\pi^{1/2}\beta^{1/2}l^{2}\tan^{2}\theta$.

Foote & Kulsrud (1979) show that in the small $l$ limit, which applies to most of the waves amplified by streaming cosmic rays, the dispersion relation can be approximated asymptotically by

[TABLE]

Expanding eqn. (7) for small $\alpha$, the damping rate is given by

[TABLE]

Using (2) for $\theta$ and restoring dimensional units, we arrive at the result

[TABLE]

[TABLE]

where we have used $lk_{0}=k\cos\theta=r_{L}^{-1}$.

Comparison with (5) shows that with this treatment the wave damping rate is enhanced by a factor of $\beta^{1/2}$. In the cluster environments simulated in Wiener et al. (2013), where $\beta$ can be of the order of 100 in the cluster centers, this factor can lead to significantly higher CR streaming speeds than those predicted by (5).

The factor of $\beta^{1/2}$ can be heuristically understood as follows. The Landau damping rate is given by the rate at which resonant ions absorb energy from oblique waves; $\Gamma_{\rm damp}^{\rm Landau}\sim kv_{i}{\rm tan}^{2}\theta$. The turbulent damping rate is given by the rate at which a pair of interacting Alfvén waves cascade, $\Gamma_{\rm damp}^{\rm turb}\sim kv_{\rm A}{\rm tan}^{2}\theta$. The geometrical factor of ${\rm tan}^{2}\theta$ is the same for these pairwise interactions. Thus, the Landau damping rate is larger by a factor $\sim v_{i}/v_{\rm A}\sim\beta^{1/2}$.

2.2 Collisionality

The analysis by Foote & Kulsrud (1979) assumes a collisionless plasma. We must verify that the plasmas in galaxy cluster environments we will be simulating are sufficiently collisionless. In particular, we must determine whether microinstabilities driven by pressure anisotropy can render the thermal ions essentially collisional. The relevant comparison to make here is the ion mean free path to the typical wavelength of the Alfvén waves in question.

The long mean free paths, relatively weak magnetic fields, and pervasive large scale turbulence in galaxy cluster plasmas make them attractive candidates for shear driven pressure anisotropy (Schekochihin & Cowley, 2006) through distortion of the magnetic field and preservation of the particles’ adiabatic invariants. We follow the notation of Kunz et al. (2011) in the following discussion.

Shear in a fluid with velocity field $\mathbf{u}$ and magnetic field direction $\mathbf{b}$ will drive pressure anisotropy while collisions will oppose it. Balancing the two yields an equilibrium anisotropy given by Braginskii (1965)

[TABLE]

In the above, the $i$ subscript indicates we are referring to ions and $\nu_{ii}$ indicates the frequency of Coulomb collisions. These collisions serve to isotropize the pressure.

A complication arises if collisions are so infrequent that the pressure anisotropy exceeds the threshold for microinstabilities:

[TABLE]

[TABLE]

If Coulomb collisions alone are not sufficient to prevent these instabilities, we may expect two possible responses of the fluid. In one scenario, the shear $S\equiv\mathbf{bb}:\nabla\mathbf{u}-1/3\nabla\cdot\mathbf{u}$ will adjust itself to reduce the pressure anisotropy until marginal stability is achieved. In this case the collision frequency does not change - it is simply the Coulomb collision frequency and the mean free path is the Coulomb mean free path. This is the working hypothesis in Kunz et al. (2011).

In the other scenario, the shear $S$ remains unchanged and instead the resulting magnetic fluctuations driven by the instabilities will themselves scatter the ions. This increases the effective total scattering frequency $\nu_{i}$, which now includes these magnetic scatterings in addition to Coulomb scattering, until marginal stability is achieved:

[TABLE]

where $\xi$ is -1 for the firehose instability, and 1/2 for the mirror instability. In this case, the total collision frequency is enhanced, and so the mean free path is shorter than the Coulomb mean free path.

The above marginal stability criterion therefore provides a lower limit on the ion scattering frequency $\nu_{i}$, and thus an upper limit on the ion mean free path $\lambda_{i}=v_{i}/\nu_{i}$. This limit depends on the shear forcing $S$.

Let us approximate this forcing as $S\sim U/L$ for some characteristic speed $U$ and length scale $L$. If we assume a turbulent cascade with Kolmogorov scaling, $U^{3}/L=\textrm{const.}$, we can relate any $U$ and $L$ to the outer scales $U_{0}$ and $L_{0}$. The marginal stability criterion then becomes

[TABLE]

The highest collision frequency will therefore be dictated by the smallest scale of turbulence, the dissipation scale $L_{d}$. This scale is determined by balancing the cascade rate $U_{d}/L_{d}$ with the viscous dissipation rate, which itself depends on the collision frequency, and is of order $v_{i}^{2}/(\nu_{i}L_{d}^{2})$. We obtain

[TABLE]

[TABLE]

[TABLE]

Combining (12) and (13) we arrive at an expression for the ion collision rate at marginal stability

[TABLE]

[TABLE]

indicating an ion mean free path of

[TABLE]

where we have defined the Mach number at the driving scale $\mathcal{M}_{0}\equiv U_{0}/v_{i}$.

We can plug this solution for $\nu_{i}$ back in to (13) to find the dissipation length:

[TABLE]

Note that the turbulent velocity at this scale is

[TABLE]

That is, under the assumption of marginal stability, the characteristic turbulent velocity at the dissipation scale is about the same as the Alfvén speed. This means the dissipation scale $L_{d}$ is about the same as the length scale $L_{\textrm{MHD}}$ defined in Wiener et al. (2013)’s streaming simulations as a measure of the wave damping rate, where it is assumed to be of order 100 kpc.

Now we can make an estimate of the mean free path as a function of $L_{\textrm{MHD}}$:

[TABLE]

This is much longer than the several AU gyroradii of the $\sim$100 GeV CR protons that generate the secondary $e^{\pm}$ which produce the observed radio emission. The assumption of collisionless wave damping is therefore sound if we assume marginal stability to MHD microinstabilities.

We can also compare the above result to the Coulomb mean free path $\lambda_{ii}$. If $\lambda_{ii}$ is significantly shorter, then Coulomb collisions alone are enough to keep the pressure anisotropy (11) low enough to avoid MHD microinstabilities. We have

[TABLE]

where we approximate $\ln\Lambda\approx 30$ and we adopt the subscript notation $Q_{x}=Q/10^{x}$ in cgs units.

In the central regions of our simulated Coma cluster (see §3), $T_{7}=9.5$ and $n_{i,-3}=3.4$, so $\lambda_{ii}\approx 4.4$ kpc. This is less than the mean free path derived from marginal stability above, implying Coulomb collisions sufficiently isotropize the pressure to avoid instabilities, at least in the cluster center. But this is still well above the AU scale gyroradii, implying we are safely in the collisionless damping regime.

The fact that $L_{\rm d}\sim L_{\rm MHD}$ when the plasma is marginally stable to microinstabilities raises potentially serious issues. It could imply that there is no MHD inertial range, since the Reynolds number at the Alfvén scale is of order unity, and parallel trans-Alfvénic motions simply dissipate. Indeed, recent analytic and numerical work finds an upper limit on shear Alfvén fluctuations of $\delta B_{\perp}/B_{0}\sim\beta^{-1/2}$ (Squire et al., 2016, 2017), above which the perturbation is rapidly quenched by the firehose instability. While CR-streaming driven turbulence lies below this limit, the background turbulence does not, and it is unclear whether equation (2) (which relies on canonical Goldreich-Sridhar theory) still applies. Such issues lie well beyond the scope of this paper, but they raise important caveats to keep in mind.

2.3 Magnetic Topology and Other Unmodeled Factors

As will be explained in §3, our numerical simulations are 1D spherically symmetric. This is primarily for two reasons - first, this simplifies the physics for ease of computation. Second, we want an apples-to-apples comparison with the simulations from Wiener et al. (2013), which also employed these simplifications for the first reason. We discuss here the possible implications of these simplifications.

Real galaxy clusters are, of course, not spherically symmetric. The Coma cluster, which we use as our characteristic cluster, has noticeable azimuthal dependence in the radio surface brightness, as reported by Brown & Rudnick (2011). We compare their 1.4 GHz observations of different quadrants of the Coma cluster using the Green Bank Telescope (provided by Larry Rudnick, personal communication) with the azimuthally averaged observations of Deiss et al. (1997) as well as the initial surface brightness in our model cluster in figure 1. There is clearly non-symmetric structure in the radio signal. The comparison between the two sets of data is further complicated by subtraction of the foreground signal from our own Galaxy, which is handled differently in each paper. As far as they relate to the construction of our initial CR profile, these differences in the data are of little consequence - streaming times in our simulations are not sensitive to small changes in the initial CR profile.

However, the departure from spherical symmetry has other, more significant implications for the evolution of out model CRs. The most important of these is the structure of the cluster’s magnetic field. Our 1D simulations necessarily have perfectly radial magnetic fields, which is a best-case scenario for streaming. In the limit of small cross-field diffusion, CRs can only stream along magnetic field lines. They can thus leave the cluster most quickly if the field lines are radial.

But in a real cluster the magnetic field may be significantly tangled. Indeed, the very basis of our streaming model is that the presence of MHD turbulence provides a wave damping mechanism. This same turbulence will also tangle the magnetic field on some length scale, increasing the escape time of CRs in our streaming model in some potentially complicated way.

Fortunately we can make reasonable estimates of how the escape time is affected by such tangling. Suppose we want to know how long it takes for a CR to travel a distance $D$ away from the center of our cluster. In the 1D symmetric case with perfectly radial magnetic field lines, this time is just

[TABLE]

where we simplify to a case where the streaming speed $v_{\textrm{st}}$ is roughly constant in space and time. If instead, the field is tangled on some length scale $L$, we can treat the transport of CRs as a random walk with steps of length $L$, travel at speed $v_{\textrm{st}}$. CR transport is then effectively diffusive with diffusion coefficient

[TABLE]

The time for CRs to diffuse out to a distance $D$ is then

[TABLE]

In other words, tangling of the field on scale $L$ increases the escape time of CRs by roughly $D/L$.

Is this a large factor? We can approximate the tangling length scale consistently by finding the scale where the kinetic energy density in turbulence equals the magnetic energy density. Below this scale, turbulent motions are too weak to bend the field lines. This occurs when

[TABLE]

Using our definition of $L_{\textrm{MHD}}$ as the length scale where the turbulent speed $v_{\textrm{t}}$ is equal to the Alfvén speed $v_{\textrm{A}}$ (see §2.2), we see that the tangling length scale described above is just of order $L_{\textrm{MHD}}$. Our fiducial value is 100 kpc, and our simulated cluster has radius 1 Mpc. We may then expect the effective escape speed to exceed the streaming times in our simulation by about a factor of 10.

There are also some observational constraints on the magnetic field structure in galaxy clusters which rely on Faraday rotation measure (RM) observations of radio sources in or behind the clusters. The theoretical background for this technique is described in detail by many authors (see Murgia et al. (2004) and Feretti et al. (2012) for just two examples). In this framework, the structure of the magnetic field is described by a simple power law in Fourier space, $|B_{k}|^{2}\propto k^{-n}$, between two length scales $\Lambda_{\textrm{min}}$ and $\Lambda_{\textrm{max}}$.

Bonafede et al. (2010) use RM observations in this way to constrain the magnetic field of the Coma cluster in particular. They claim that the field which best fits the RM data is tangled on scales ranging from $\Lambda_{\textrm{min}}\sim 2$ kpc to $\Lambda_{\textrm{max}}\sim 34$ kpc, with a steep Kolmogorov-like spectrum of $n=11/3$. This suggests that most of the power is on large scales, i.e. most of the “steps” in the field line random walk are of length $\Lambda_{\textrm{max}}\sim 34$ kpc. This is reasonably close to our above estimate of $L_{\textrm{MHD}}$. We note, however, that in general the magnetic field structure of galaxy clusters is not well constrained by observations. So while it is hard to say with any certainty, we may reasonably expect our simplified model to underestimate streaming times by a factor of 10 - 30.

In the context of our question about radio halo turnoff times, this is not a cripplingly large factor. It extends the turnoff time from hundreds of Myr to a few Gyr, still in the range of reasonable turnoff times for the hadronic model. More importantly for the context of this work, this factor is independent of the high beta effects considered here. If our older simulations from Wiener et al. (2013) are off by a factor of 10, then the new simulations presented here are also off by this same factor, and our main point is unchanged.

Of course there are non-linear effects which complicate this - as the CR density drops the streaming speed increases. If escape speeds are initially slower by a factor of 10, the overall shutoff time may be changed by an entirely different factor. Determining the effects of this non-linearity would require non-1D simulations which are beyond the scope of this work. However, we are encouraged that in the apples-to-apples comparison we make here, the high-beta effect causes an increase in the initial streaming speeds by a factor of $\beta^{1/2}\sim 8$ in the central regions of the cluster. So differences brought about by these high beta effects should at least be comparable to the effects of field tangling.

Still, it is difficult to assess the effects of field tangling and non-linear evolution other than the above speculation. If we have overestimated the field tangling scale by a factor of about 10, which is possible considering the spectrum in Bonafede et al. (2010) extends down to 2 kpc, then our escape speed becomes of the order of the age of the universe even without accounting for the non-linearity of the evolution. Future observational constraints on magnetic field structure in galaxy clusters may therefore prove critical to the question of long range CR streaming.

There is another point to be made about the magnetic field topology which may be important. We have talked about the effects of a tangled field, but what if the field lines never leave the cluster? Consider a single field line which extends out to some radius $R$ in the cluster before folding back on itself and returning towards the center. Then along this field line, streaming will even out the CR density out to $R$, but CRs will be unable to leak out past $R$. If all the field lines worked this way, no CRs would be able to escape, although they could spread out evenly within radius $R$. However, if some percentage of field lines leave the cluster, then CRs will always be able to leak out along these lines. The halo dropoff time would then be primarily determined by two factors - the percentage of field lines which exit the cluster, and the rate of cross field diffusion which allows CRs on trapped lines to migrate to neighboring lines which escape.

3 Simulation Setup

We reproduce here the spherically symmetric ZEUS hydrodynamic simulations of the Coma cluster from Wiener et al. (2013), but with the damping rate enhanced by the factor of $\beta^{1/2}$ from the analysis in §2. As mentioned in §2.3, spherical symmetry is clearly a drastic simplification (L. Rudnick; personal communication) but we assume it here to keep the problem tractable and to focus on the difference between the present treatment and Wiener et al. (2013). As in Wiener et al. (2013), we model the Coma cluster with density and temperature profiles given by

[TABLE]

[TABLE]

These profiles are taken from Pinzke & Pfrommer (2010), with the density inferred from X-ray observations (Briel et al. (1992)).

The magnetic field is assumed to scale with gas density as

[TABLE]

with $B_{0}=5\ \mu$G and $\alpha_{B}=0.3$, as suggested by constraints from Bonafede et al. (2010) (Wiener et al. (2013) also considered $\alpha=0.5$ which is more in line with Bonafede et al. (2010); the choice of $\alpha=0.3$ here is to properly compare with the results from Wiener et al. (2013). It is conservative in the sense that it implies a lower $\beta$ and therefore a weaker Landau damping effect). The above density, temperature, and magnetic field profiles are fixed for these simulations - only the CR distribution is evolved in time.

The initial CR distribution is of the form used in Pinzke & Pfrommer (2010), motivated by cosmological hydrodynamic simulations of galaxy clusters where cosmic rays are accelerated via diffusive shock acceleration:

[TABLE]

[TABLE]

Here and in the equations to follow, the momenta $p$ are expressed in units of $mc$ for the appropriate $m$. Namely, for actual momentum $P$ we will work in terms of $p_{p}=P/m_{p}c$ and $p_{e}=P/m_{e}c$. In this framework we have to be careful how we define our distribution functions $f_{p}$ and $f_{e}$. To be unambiguous, let us define them as such:

[TABLE]

[TABLE]

where d$n_{p}(r,p_{p})$ is the differential number density of CRp at radius $r$ and unitless momentum $p_{p}$, and similarly for the electrons.

We choose the normalization $C$ to be of the form 444Although this resembles the normalization in Pinzke & Pfrommer (2010), please note the difference between our $C(r)$ (which have dimensions of number density) used here and the $\tilde{C}(r)=C(r)\rho(r)/m_{p}$ (unitless) used by Pinzke & Pfrommer (2010). We utilize this formula as a convenience, but our CR density profile is actually much flatter than in their simulation.

[TABLE]

In Pinzke & Pfrommer (2010), the parameters $C_{\textrm{vir}}$, $C_{\textrm{center}}$, $r_{\textrm{trans}}$, and $\beta_{\textrm{C}}$ are determined from scaling relations. We instead choose values that roughly reproduce the observed synchrotron radiation. For Coma, these are $C_{\rm center}=6\times 10^{-11}\,{\rm cm^{-3}}$, $C_{\rm vir}=5.2\times 10^{-11}\,{\rm cm^{-3}}$, $r_{\rm trans}=55$ kpc, $\beta_{\rm C}=1.09$.

With the above initial distribution we evolve the CRp distribution function forward in time according to

[TABLE]

as in Skilling (1971). In the above, $\mathbf{u}$ and $\mathbf{v_{\textrm{A}}}$ are the gas and Alfvén velocities, $\Omega_{0}$ is the non-relativistic gyrofrequency, and $\mathbf{n}$ is a unit vector which points along the magnetic field. The last term represents diffusion with respect to the Alfvén wave frame due to wave damping, and is directly affected by the discussion in §2. More details of this equation and its numerical evolution can be found in Wiener et al. (2013). We set $\mathbf{u}\equiv 0$ in the simulations discussed here, but include it in eqn. (25) for completeness.

As the CRp distribution function evolves, we derive from it at every time step a steady-state secondary CRe distribution function according to

[TABLE]

which balances the source function from the CRp hadronic collisions

[TABLE]

with the losses from synchrotron emission and inverse Compton (IC) scattering

[TABLE]

In (27), $n_{N}$ is the number density of target nucleons and we assume we are far above the pion production threshold.

This steady state model includes two implicit assumptions. First, the energy losses of the secondary CRe are dominated by synchrotron and IC losses. Namely this means we assume the CRe don’t stream on time scales faster than their loss times ( \sim$$<  100 Myr), a condition we will have to check a posteriori. This assumption will turn out to be well satisfied. Second, we assume that the source function $s_{e}$ is not changing significantly on these same time scales. This will turn out not to be very well satisfied, as we discuss below in §4.

Once we have the secondary CRe distribution function at each time step, we can then determine the synchrotron emissivity at frequency $\nu$. From Rybicki & Lightman (1979),

[TABLE]

In the above, $\nu_{c}=3eBp_{e}^{2}/4\pi m_{e}c$ and the function $F$ is an integral of a modified Bessel function, $F(x)=x\int_{x}^{\infty}K_{5/3}(x^{\prime})\textrm{d}x^{\prime}$. With the emissivity in hand, we determine surface brightness and total luminosity from simple spatial integrals555We reiterate that the hadronic model for giant radio haloes is disfavored by gamma-ray non-detections. We include a radio surface brightness prediction for direct comparison with Wiener et al. (2013). The dimming of the model radio halo can also be thought of as a simple proxy for the evolution of CRp..

4 Results

The total 1.4 GHz luminosity of our simulated Coma cluster as a function of time is shown in Figure 2. Three CR transport models are compared: in the first, there is no wave damping and CRs stream at the Alfvén speed. In the second, there is wave damping according to (5). In the third, there is wave damping with the high-$\beta$ correction factor, (9). We see that the increased factor of $\beta^{1/2}$ in the damping rate allows CRs to stream out even faster than in our original simulation, causing the radio luminosity to also drop on faster time scales. The streaming speeds for the relevant 100 GeV CRp are shown in figure 3. This figure shows the streaming speed at a fixed radius of 300 kpc as a function of time. Initially the streaming speeds are larger with the high-$\beta$ correction in comparison to without, and this enhancement grows with time due to the non-linear evolution of the streaming speeds.

We may wonder if these enhanced streaming speeds interfere with any of our assumptions, namely our assumption that CRe secondary losses are dominated by synchrotron and inverse Compton losses. CRe will stream at the same speeds as CRp of the same energy, and if the CRe stream out on time scales comparable to their loss times, our steady state treatment is incorrect. Looking at Figure 3 we may suspect that this is the case, as the $\sim 10^{4}$ km/s streaming speeds imply streaming times of $\sim$100 Myr across our 1-2 Mpc cluster.

However, for the electrons, the relevant energies for 1.4 GHz emission are 1-10 GeV. As shown in Figure 4, CRp at these energies stream out much more slowly, even in our high-$\beta$ model666Note that most CR energy resides at $\sim$GeV energies, and that even with the high $\beta$ correction these CRs remain more or less locked to the wave frame. Thus, CR heating in cluster cores (Loewenstein et al., 1991; Guo & Oh, 2008; Pfrommer, 2013; Ruszkowski et al., 2017) remains viable.. We can expect the streaming times of our CRe secondaries to be $\sim$1-2 Gyr, much longer than their loss times. Note that since the streaming instability growth rate is proportional to $n_{\rm CR}$ (equation (1)), and $n_{\rm CR,p}\gg n_{\rm CR,e}$, the CRe actually scatter off the wave field created by the much more abundant CRp. If they manage to achieve spatial separation, the CRe could in principle stream much faster.

Still, the streaming speeds of the 100 GeV CRp which source these electrons are high, implying a short crossing time. Our steady state model assumes that the CRp distribution changes on time scales much longer than the CRe loss times, such that the production of secondaries ‘quickly’ reaches an equilibrium before the CRp density changes too much. This assumption is no longer satisfied, so a more complete treatment of secondary CRe production and transport may be necessary. However, we have used somewhat artificial initial conditions which do not reflect the cosmological build up of CRp. As long as the injection time is longer than the streaming time (likely to be true for CRs sourced by structure formation (giant radio halos), less so for CRs sourced by AGN activity (radio mini-halos)), CRs will have a flat distribution, which slowly increases in normalization.

5 Conclusion

Despite significant advances in observations and modeling, the origin and transport of cosmic rays in galaxy clusters remains incompletely understood. As long as this is the case, radiative diagnostics based on cosmic rays will be provisional, and a full assessment of the role of cosmic rays in the dynamics and energy balance of clusters will elude us.

In this paper, we have added a new ingredient to the propagation problem: the role of ion Landau damping in suppressing the growth of Alfvén waves excited by cosmic ray streaming. Although waves which propagate exactly along the background magnetic field $\mathbf{B}_{0}$ are undamped, inhomogeneity in $\mathbf{B}_{0}$ precludes perfectly parallel propagation. At the large $\beta\equiv 8\pi P_{g}/B^{2}$ values characteristic of galaxy clusters, oblique waves are subject to strong Landau damping due to their parallel electric fields.

When the propagation angle is calculated assuming $\mathbf{B}_{0}$ is a uniform field with a superimposed anisotropic MHD turbulent cascade, the resulting damping rate (9) is of the same form as the damping rate estimated from the turbulent shearing rate (5), which was taken as the dominant damping mechanism in earlier work (Wiener et al., 2013), but exceeds the turbulent damping rate by a factor of $\beta^{1/2}$. It is worth noting that the assumption of an anisotropic turbulent cascade is key to this result: if the scale of variation of $\mathbf{B}_{0}$ were taken to be a global scale, the resulting damping rate would be very slow (Zweibel, 2003).

Since Landau damping is a collisionless process, it is important to check whether the effective ion mean free path is much longer than the wavelengths of the Alfven waves. Estimating the mean free path $\lambda_{ii}$ to Coulomb scattering is straightforward and the result (18) shows that the waves are indeed collisionless with respect to Coulomb interactions. However, we also estimated for the first time the mean free path $\lambda_{i}$ due to scattering from microscale instabilities driven by pressure anisotropy (17), and found that it too is much larger than the wavelengths of the Alfvén waves in question, fully validating the role of Landau damping in suppressing the Alfvén waves that scatter cosmic rays in turbulent galaxy cluster plasmas. Note that there are potentially important caveats about the nature of MHD turbulence in a high $\beta$ plasma (Squire et al., 2016, 2017), which are beyond the scope of this paper.

We then reconsidered the evolution of an idealized version of the radio halo of the Coma cluster following Wiener et al. (2013), but including Landau damping. In this model, the halo is produced by secondary electrons and positrons created in hadronic interactions between cosmic ray protons and thermal cluster gas. As expected, transport of the cosmic ray protons is much faster, quickly depleting them and turning off the halo on even faster timescales than predicted in the original model. This strengthens the case for rapid evolution of radio halos, at least in this idealized (1D, radial magnetic field) model cluster, unless the cosmic ray primary proton source is continuously replenished. This setup is the most optimistic case there is for evolution via streaming, an important caveat of this work. Future simulations with higher dimensionality and more complicated field topologies are necessary to study this effect in real clusters.

The problem of cosmic ray electron transport on galaxy clusters was first raised by Jaffe (1977), who pointed out that the radiative loss time of cosmic ray electrons is much shorter than their transport time at the Alfvén speed from the core of the cluster. If the primary CRe could stream out to large distances in less than one cooling time, they could supply the radio emission in the outskirts without being reaccelerated. However, even in our ideal scenario this does not seem to be the case. Although we have not undertaken a detailed comparison of Jaffe’s model with ours, the transport speeds of 5 GeV electrons shown in Figure 4 are well below the $\sim 2000$ km/s transport speed that Jaffe estimated was necessary to form the halo with primary electrons.

The rapid transport speed of cosmic ray protons due to Landau damping of Alfvén waves could help to explain the stringent upper limits on diffuse $\gamma$-ray emission from galaxy cluster cores reported by Pinzke & Pfrommer (2010); Enßlin et al. (2011); Ahnen et al. (2016). It also has important consequences for hadronic production of CRe ‘seeds’ for turbulent reacceleration - rapid streaming produces a flat CRp profile, giving a seed population which better fits observational constraints (Pinzke et al., 2017). It is also relevant to the tension between observed radio halo luminosities and models of CR heating of cluster cores (Jacob & Pfrommer, 2017), since high energy CRs responsible for the former stream at much higher velocities than the lower energy CRs relevant for the latter. These issues are beyond the scope of the current paper, but a topic for future work.

Acknowledgements:

We are happy to acknowledge discussions with Matt Kunz, Eliot Quataert, and Larry Rudnick. JW and EGZ acknowledge support by NSF Grant AST-1616037, the WARF Foundation, and the Vilas Trust. SPO acknowledges support from NASA grant NNX15AK81G.