Constraints on Cosmic Ray Transport in Galaxy Clusters from Radio and Gamma Ray Observations

TL;DR

This paper uses radio and gamma-ray observations of galaxy clusters to constrain cosmic ray transport properties, deriving a minimum diffusion coefficient and limits on proton acceleration efficiency.

Contribution

It provides the first quantitative constraints on cosmic ray transport mechanisms in galaxy clusters based on observational data.

Findings

Minimum diffusion coefficient in Coma cluster is ~10^{31} cm^2/s.

Protons cannot be accelerated more than 15 times efficiently than electrons.

Constraints on magnetic field topology and cosmic ray loss mechanisms.

Abstract

The nature of cosmic rays (CRs) and cosmic ray transport in galaxy clusters is probed by a number of observations. Radio observations reveal the synchrotron radiation of cosmic ray electrons (CRe) spiraling around cluster magnetic fields. $\gamma$-ray observations reveal hadronic reactions of cosmic ray protons (CRp) with ambient gas nuclei which produce pions. To date, no such cluster-wide $\gamma$-ray signal has been measured, putting an upper limit on the density of CRp present in clusters. But the presence of CRe implies some source of CRp, and consequently there must be some CRp loss mechanism. In this paper we quantify the observational constraints on this loss mechanism assuming that losses are dominated by CR transport, ultimately deriving a minimum diffusion coefficient of $\sim10^{31}$ cm$^2$ s$^{-1}$ in the Coma cluster. This lower limit on transport may help illuminate some unknown properties of the cluster field topology. Conversely, measurements of cluster field tangling scales can constrain other model parameters, such as the relative acceleration efficiency of protons to electrons. To be consistent with the Coma observations, protons cannot be accelerated more than 15 times more efficiently than electrons of the same energy.

Figures

Tables

Table 1. Table 1: Results for varying different model parameters with non-explicit dependencies. a Central magnetic field strength, defined in ( 32 ) b Magnetic field dependence on gas density, defined in ( 32 ). c Shape function dependence on gas density, defined in ( 33 ). d Spectral index of radio emission. e Maximum extent of cluster. f Cluster core radius, defined in ( 34 ). g Exponent of gas density profile, defined in ( 34 ). h Fraction of initial CRp kinetic energy contained in each pion produced in hadronic reactions, defined in ( 42 ). i Pion multiplicity (number of pions produced per hadronic reaction), defined in ( 42 ). j Minimum velocities according to ( 51 ). k Minimum diffusion coefficients according to ( 52 ). l Gamma ray energy flux predicted in the 7.5-10 GeV band, assuming CRe are dominated by secondaries, given by ( 55 ).

$B_0^a$	${\alpha }_B^b$	${\alpha }_{\eta }^c$	$s^d$	$r_{\text{max}}^e$	$r_c^f$	${\beta }^g$	$K_{\pi }^h$	${\xi }^i$	$v_{\text{t,min}}^j$	${\kappa }_{\text{min}}^k$	${\Phi }_{\gamma }^l$(7.5-10 GeV)
($\mu$G)				(Mpc)	(kpc)				(km s-1)	(1030 cm2s-1)	(10-10 GeV cm-2s-1)
1	0.3	0.85	1.35	4	294	0.75	0.17	3	140	65	29
2	0.3	0.85	1.35	4	294	0.75	0.17	3	29	13	5.9
3	0.3	0.85	1.35	4	294	0.75	0.17	3	13	5.7	2.4
6	0.3	0.85	1.35	4	294	0.75	0.17	3	4.2	1.6	0.60
10	0.3	0.85	1.35	4	294	0.75	0.17	3	2.4	0.84	0.27
1	0.5	0.38	1.35	4	294	0.75	0.17	3	520	350	250
2	0.5	0.38	1.35	4	294	0.75	0.17	3	110	70	50
3	0.5	0.38	1.35	4	294	0.75	0.17	3	42	27	19
6	0.5	0.38	1.35	4	294	0.75	0.17	3	10	6.0	4.0
10	0.5	0.38	1.35	4	294	0.75	0.17	3	4.3	2.2	1.3
1	0.7	-0.09	1.35	4	294	0.75	0.17	3	3100	2600	2700
2	0.7	-0.09	1.35	4	294	0.75	0.17	3	610	510	530
3	0.7	-0.09	1.35	4	294	0.75	0.17	3	240	200	210
6	0.7	-0.09	1.35	4	294	0.75	0.17	3	48	39	41
10	0.7	-0.09	1.35	4	294	0.75	0.17	3	16	12	12
3	0.5	0.0	1.35	4	294	0.75	0.17	3	51	41	40
3	0.5	0.2	1.35	4	294	0.75	0.17	3	47	34	28
3	0.5	0.38	1.35	4	294	0.75	0.17	3	42	27	19
3	0.5	0.6	1.35	4	294	0.75	0.17	3	37	20	12
3	0.5	1.0	1.35	4	294	0.75	0.17	3	30	12	4.3
3	0.5	0.48	1.15	4	294	0.75	0.17	3	80	48	7.4
3	0.5	0.43	1.25	4	294	0.75	0.17	3	58	36	12
3	0.5	0.38	1.35	4	294	0.75	0.17	3	42	27	19
3	0.5	0.33	1.45	4	294	0.75	0.17	3	31	21	31
3	0.5	0.28	1.55	4	294	0.75	0.17	3	23	16	48
3	0.5	0.38	1.35	1	294	0.75	0.17	3	14	1.8	1.8
3	0.5	0.38	1.35	2	294	0.75	0.17	3	24	7.1	5.4
3	0.5	0.38	1.35	4	294	0.75	0.17	3	42	27	19
3	0.5	0.38	1.35	8	294	0.75	0.17	3	75	100	76
3	0.5	0.38	1.35	16	294	0.75	0.17	3	130	380	310
3	0.5	0.38	1.35	4	100	0.75	0.17	3	35	25	160
3	0.5	0.38	1.35	4	200	0.75	0.17	3	39	27	41
3	0.5	0.38	1.35	4	294	0.75	0.17	3	42	27	19
3	0.5	0.38	1.35	4	400	0.75	0.17	3	45	28	11
3	0.5	0.38	1.35	4	500	0.75	0.17	3	47	28	7.2
3	0.5	0.22	1.35	4	294	0.55	0.17	3	55	31	10
3	0.5	0.31	1.35	4	294	0.65	0.17	3	51	31	15
3	0.5	0.38	1.35	4	294	0.75	0.17	3	42	27	19
3	0.5	0.43	1.35	4	294	0.85	0.17	3	33	22	21
3	0.5	0.47	1.35	4	294	0.95	0.17	3	25	16	20
3	0.5	0.38	1.35	4	294	0.75	0.50	1	88	57	19
3	0.5	0.38	1.35	4	294	0.75	0.25	2	54	35	19
3	0.5	0.38	1.35	4	294	0.75	0.17	3	42	27	19
3	0.5	0.38	1.35	4	294	0.75	0.13	4	36	23	19
3	0.5	0.38	1.35	4	294	0.75	0.10	5	29	19	19

Equations

$${\text{d}}^6{N}_i(\mathbf{r},\mathbf{p})=f_i(\mathbf{r},\mathbf{p}){\text{d}}^3\mathbf{r}{\text{d}}^3\mathbf{p},$$

$$\frac{\partial f_i(\mathbf{r},\mathbf{p})}{\partial t}+\nabla \cdot (f_i\mathbf{v})+{\nabla }_{\mathbf{p}}\cdot \left(f_i\frac{\text{d}\mathbf{p}}{\text{d}t}\right)=0.$$

$$\begin{gathered} \frac{\partial f_i(r,E_i)}{\partial t}+\frac{1}{r^2}\frac{\partial }{\partial r}(r^2{F}_{cr,i}(r,E_i)) \\ +\frac{\partial }{\partial E_i}\left({\dot{E}}_i(r,E_i)f_i(r,E_i)\right)=s_i(r,E_i), \end{gathered}$$

$${\text{d}}^2{N}_i(r,E_i)=4\pi r^2{f}_i(r,E_i)\text{d}r\text{d}{E}_i.$$

$$\frac{\partial }{\partial E_e}\left({\dot{E}}_e(r,E_e)f_e(r,E_e)\right)=s_e(r,E_e).$$

$$f_e(r,E_e)=\frac{1}{|{\dot{E}}_e|}{\int }_{E_e}^{\infty }{s}_e(r,E_e^{\prime })\text{d}{E}_e^{\prime }.$$

$$f_e(r,{\gamma }_e)=f_e(r,E_e)\frac{\text{d}{E}_e}{\text{d}{\gamma }_e}=m_e{c}^2{f}_e(r,E_e),$$

$$s_e(r,{\gamma }_e)=m_e{c}^2{s}_e(r,E_e).$$

$$f_e(r,{\gamma }_e)=\frac{1}{|{\dot{\gamma }}_e|}{\int }_{{\gamma }_e}^{\infty }{s}_e(r,{\gamma }_e^{\prime })\text{d}{\gamma }_e^{\prime }.$$

$$\frac{1}{r^2}\frac{\partial }{\partial r}(r^2{F}_{cr,p}(r,E_p))=s_p(r,E_p).$$

$$F_{cr,p}=\frac{1}{r^2}{\int }_0^r{r}^{\prime 2}{s}_p(r^{\prime },E_p)\text{d}{r}^{\prime }.$$

$$F_{cr,p}(r,E_p)\equiv f_p(r,E_p)v_t(r,E_p).$$

$$f_p(r,E_p)=\frac{1}{r^2{v}_t(r,E_p)}{\int }_0^r{r}^{\prime 2}{s}_p(r^{\prime },E_p)\text{d}{r}^{\prime }.$$

$$F_{cr,p}(r,E_p)\equiv -\kappa (r,E_p)\frac{\partial f_p(r,E_p)}{\partial r}$$

$$\frac{\partial f_p}{\partial r}=-\frac{1}{r^2\kappa (r,E_p)}{\int }_0^r{r}^{\prime 2}{s}_p(r^{\prime },E_p)\text{d}{r}^{\prime }.$$

$$f_p(r,E_p)={\int }_r^{r_{\text{max}}}\frac{\text{d}{r}^{\prime \prime }}{(r^{\prime \prime }{)}^2\kappa (r^{\prime \prime },E_p)}{\int }_0^{r^{\prime \prime }}{r}^{\prime 2}{s}_p(r^{\prime },E_p)\text{d}{r}^{\prime }.$$

$$\mathbf{u}={\mathbf{v}}_0+⟨\frac{3}{2}(1-{\mu }^2)\frac{{\nu }_{+}-{\nu }_{-}}{{\nu }_{+}+{\nu }_{-}}⟩v_A\mathbf{n}$$

$$\begin{gathered} \frac{\partial {\bar{f}}_i}{\partial t}+\frac{\partial }{\partial p^3}(p^3\mathbf{u})\cdot \nabla {\bar{f}}_i=(\nabla \cdot \mathbf{u})p^3\frac{\partial {\bar{f}}_i}{\partial p^3}+ \\ \nabla \cdot (D_{\parallel }\mathbf{nn}\cdot \nabla {\bar{f}}_i)+\frac{\partial }{\partial p^3}\left(9p^4{D}_{\text{pp}}\frac{\partial {\bar{f}}_i}{\partial p^3}\right) \end{gathered}$$

$$D_{\parallel }=v^2⟨\frac{1-{\mu }^2}{2({\nu }_{+}+{\nu }_{-})}⟩,$$

$$D_{\text{pp}}=4{\gamma }^2{m}^2{v}_A^2⟨\frac{1-{\mu }^2}{2}\frac{{\nu }_{+}{\nu }_{-}}{{\nu }_{+}+{\nu }_{-}}⟩.$$

$${\bar{f}}_i(\mathbf{r},p)=\frac{1}{4\pi }\iint {\text{d}}^2{\Omega }_p\frac{{\text{d}}^{6}N}{{\text{d}}^3\mathbf{r}{\text{d}}^3\mathbf{p}}=\frac{1}{p^2}\frac{{\text{d}}^{4}N}{{\text{d}}^3\mathbf{r}\text{d}p}.$$

$$f_i(\mathbf{r},E)=p^2{\bar{f}}_i(\mathbf{r},p)\frac{\text{d}p}{\text{d}E}\approx \frac{p^2}{c}{\bar{f}}_i(\mathbf{r},p).$$

$$\frac{\partial }{\partial p^3}(p^3\mathbf{u})\cdot \nabla {\bar{f}}_i=\nabla \cdot \left[\frac{{\bar{f}}_i}{3p^2}\frac{\partial }{\partial p}(p^3\mathbf{u})\right]-\frac{{\bar{f}}_i}{3p^2}\frac{\partial }{\partial p}(p^3\nabla \cdot \mathbf{u}).$$

$$(\nabla \cdot \mathbf{u})p^3\frac{\partial {\bar{f}}_i}{\partial p^3}=\frac{1}{p^2}\frac{\partial }{\partial p}\left(\frac{p^3}{3}{\bar{f}}_i\nabla \cdot \mathbf{u}\right)-\frac{{\bar{f}}_i}{3p^2}\frac{\partial }{\partial p}\left(p^3\nabla \cdot \mathbf{u}\right).$$

$$\begin{gathered} \frac{\partial {\bar{f}}_i}{\partial t}+\nabla \cdot \left[\frac{{\bar{f}}_i}{3p^2}\frac{\partial }{\partial p}(p^3\mathbf{u})\right]-\frac{1}{p^2}\frac{\partial }{\partial p}\left(\frac{p^3}{3}{\bar{f}}_i\nabla \cdot \mathbf{u}\right)= \\ \nabla \cdot (D_{\parallel }\mathbf{nn}\cdot \nabla {\bar{f}}_i)+\frac{1}{p^2}\frac{\partial }{\partial p}\left(p^2{D}_{\text{pp}}\frac{\partial {\bar{f}}_i}{\partial p}\right), \end{gathered}$$

$$\begin{gathered} \frac{\partial f_i}{\partial t}+\nabla \cdot \left[\frac{f_i}{3E^2}\frac{\partial }{\partial E}(E^3\mathbf{u})-D_{\parallel }\mathbf{nn}\cdot \nabla f_i\right] \\ -\frac{\partial }{\partial E}\left[\frac{E}{3}{f}_i\nabla \cdot \mathbf{u}+D_{\text{pp}}\frac{\partial f_i}{\partial E}\right]=0. \end{gathered}$$

$$F_{cr,\text{adv}}=\frac{f_i}{3E^2}\frac{\partial }{\partial E}(E^3\mathbf{u})\equiv f_i{\mathbf{v}}_{t,\text{stream}}$$

$$\Rightarrow {\mathbf{v}}_{t,\text{stream}}(\mathbf{r},E)=\frac{1}{3E^2}\frac{\partial }{\partial E}(E^3\mathbf{u}),$$

$$F_{cr,\text{diff}}=-D_{\parallel }\mathbf{nn}\cdot \nabla f_i\equiv -{\mathsf{K}}_{\text{stream}}(\mathbf{r},E)\cdot \nabla f_i$$

$$\Rightarrow {\mathsf{K}}_{\text{stream}}(\mathbf{r},E)=D_{\parallel }\mathbf{nn}.$$

Full text

1 Introduction

Galaxy clusters are the largest ($\sim 10^{14}-10^{15}M_{\odot}$) gravitationally bound objects in the Universe and are host to a plethora of physical processes ranging orders of magnitude in scale length and energy. The energy budget of clusters includes several non-thermal components, such as magnetic fields and cosmic rays (CRs), which are potentially important in cluster dynamics. For instance, wave heating by CRs has been proposed as a heating mechanism to prevent cooling catastrophes in cool core (CC) clusters (Loewenstein et al. (1991); Guo & Oh (2008)).

Some information about the magnetic and CR content of clusters can be determined from radio observations, which in some clusters reveal large scale diffuse synchrotron emission from cosmic ray electrons (CRe). These giant radio halos provide insight into the nature of cluster-wide magnetic fields and CR transport. The Coma radio halo in particular has been the target of several radio observations and studies (see Deiss et al. (1997), Thierbach et al. (2003), Brown & Rudnick (2011), and Brunetti et al. (2012) for just a few examples).

Cosmic ray protons (CRp) can be independently detected from their hadronic reactions - neutral pions produced in high energy hadronic collisions decay into $\gamma$-rays (charged pions also produced in these collisions decay into CR electrons and positrons, referred to as ‘CR secondaries’). Since $\gamma$-radiation can also be produced as inverse Compton emission from high energy electrons up-scattering cosmic microwave background (CMB) photons, any detection of $\gamma$-rays is only an upper bound on the rate of hadronic reactions. However, while $\gamma$-radiation has been seen in several galaxies, no diffuse cluster-wide $\gamma$-ray emission has yet been detected despite deep searches with $\gamma$-ray telescopes such as the Fermi Large Area Telescope (see Huber et al. (2013) and Ackermann et al. (2016) for just two examples). These non-detections put strict upper limits on the CRp content of galaxy clusters.

In this paper we describe a theory of CRs which combines a given radio halo detection with a $\gamma$-ray flux upper limit to derive a minimum CRp transport speed. In simple terms, if CRs are being accelerated in a cluster at a rate consistent with its radio emission, they must escape the cluster on short enough time scales in order to bring the CRp density low enough to explain the lack of $\gamma$-rays. This escape rate depends on the assumed parameters of our model, including magnetic field topology and relative acceleration efficiency of protons to electrons. We will specify to the Coma cluster, but the analysis is general.

It should be noted that ours is far from the first study of cosmic ray transport in galaxy clusters. Enßlin et al. (2011) assumed radial transport at the local sound speed. Wiener et al. (2013) and Wiener et al. (2018) assumed radial transport at a speed determined by equating the rate at which cosmic rays excite the waves that confine them through the streaming instability to the rate at which the waves are damped in the cluster plasma (see Zweibel (2017) for a review). These studies are appropriate when the magnetic field is relatively well ordered and the cosmic rays stream down their density gradient at all times. Other work has, in contrast, suggested that CR diffusion times are of order the age of the universe, and so CRs are effectively confined to their host galaxy clusters (Völk et al. (1996), Berezinsky et al. (1997), Brunetti & Jones (2014)). This can occur if there are enough externally-driven MHD waves that the self-confinement picture fails, or if the magnetic field is tangled enough that even CRs free streaming at the speed of light do not leave the cluster in a Hubble time. Here we do not ascribe to any particular magnetic field topology or level of turbulence, but instead seek to translate observations of radio and $\gamma$-ray emission into constraints on a set of cluster parameters which includes them. Additional observational constraints on the magnetic field in a cluster would then serve, combined with this analysis, to constrain the remaining model parameters, and vice versa.

In §2 we briefly review the theory of CR transport in the regime relevant to galaxy clusters. We establish a framework for discussing the evolution and steady state conditions of CR electrons and protons, and make some limiting assumptions. We broadly consider two transport models, bulk advection or streaming at speed $v_{t}$, and diffusion with diffusivity $\kappa$.

In §3 we review the connection between synchrotron emission and the underlying CRe population. From the implied CRe production rate we determine the CRp population in two limits. In §3.4.1, we assume that the CRe are all primaries (i.e. accelerated by the same processes that produce CRp). In this limit, which is akin to the re-acceleration model of radio halos, the CRp production rate is determined, and we must model CRp transport in order to predict an equilibrium CRp density to compare against $\gamma$-ray observations. In §3.4.2, we assume instead that the CRe are all secondary particles, i.e. the hadronic model of radio halos. In this limit the predicted CRe production rate translates directly into a CRp density and no constraints on CRp transport are made.

In §4 we review the connection between $\gamma$-ray emission and CRp to derive constraints on the CRp population from the upper limits on the $\gamma$-ray flux. We combine the analysis of §3 to derive (parameter-dependent) lower limits on CR transport in the case where CRe are dominated by primaries. This is pictorialized in the simplified schematic of CR processes shown in figure 1b. The opposite case, where CRe are dominated by secondaries, is pictorialized in the simplified schematic shown in figure 1c. In this case, rather than limits on CR transport, we get predictions of $\gamma$-ray emission to compare directly with observations. In §5 we apply the analysis in §§3 and 4 to the Coma cluster. In §6 we check some of the assumptions made in the analysis for consistency. We find that the minimum transport rates derived for CRp in the primary model are so slow that appreciable secondaries would be generated before the CRp escape.

In §7 we summarize and draw conclusions. Our main result is that if the cosmic ray electrons are primary particles, then with standard assumptions about the magnetic field strength and spatial profile in Coma, and Milky Way like assumptions about the ratio of proton to electron cosmic ray sources, the streaming speeds calculated in Wiener et al. (2013) and Wiener et al. (2018) are high enough to prevent proton cosmic rays from generating a detectable $\gamma$-ray signal provided the magnetic field is not tangled below about 100 kpc. But, if the magnetic field is dominated by turbulence on 10s of kpc scales, as suggested by the observations of Bonafede et al. (2010), cosmic rays build up in intensity to a level that should produce a detectable $\gamma$-ray signal for certain choices of our model parameters. The $\gamma$-ray non-detections therefore impose limits on these parameters. One potential use of this limit is to constrain the relative acceleration efficiency of protons to electrons in a cluster. In the opposite limit where the cosmic ray electrons are all secondaries, we find that the radio detections and $\gamma$-ray non-detections require a rather strong magnetic field.

2 Cosmic Ray Transport

In this section we introduce the basic CR transport equations. We first set up general evolution equations that we will use in our analysis later on. These are general in the sense that they encompass a large family of specific evolution models. We follow this with a brief description of one of these specific models, the self-confinement regime of CR transport. See Zweibel (2017), Skilling (1971), longair92, or Kulsrud (2005) for some examples of further review of CR transport equations beyond what is covered here.

2.1 A Generalized Treatment

The population of cosmic rays in a system can be completely described by its phase-space density function $f_{i}(\mathbf{r},\mathbf{p})$, which is defined by

[TABLE]

where d${}^{6}N_{i}$ indicates the number of CRs in the phase-space volume $\textrm{d}^{3}\mathbf{r}\textrm{d}^{3}\mathbf{p}$ centered at $(\mathbf{r},\mathbf{p})$. The subscript $i$ indicates different CR species ($e$ for electrons, $p$ for protons, etc.). The evolution of such a CR population is then most generally described by the Vlasov equation, which expresses the conservation of particles in phase-space:

[TABLE]

(As written, equation (2) does not account for injection of new particles; we add a source term below, in equation(3).)

We simplify equation (2) by imposing the limits inherent in our system. We will approximate our galaxy clusters as spherically symmetric systems, so we can integrate (2) over the two angular coordinates. Furthermore, in the limit of small mean free path due to scattering, CRs are nearly isotropic in momentum space. Equation (2) then reduces to

[TABLE]

where the reduced distribution function is now defined by

[TABLE]

In the above, $\textrm{d}^{2}N_{i}$ represents the number of CRs of species $i$ between cluster radii $r$ and $r+\textrm{d}r$ and having energy between $E$ and $E+\textrm{d}E$. In (3), $F_{cr,i}$ is the spherically-averaged CR number flux per unit energy due to transport, which will be determined by the given CR transport model. Changes in energy are represented by the third term, with $\dot{E_{i}}$ being the rate of energy loss or gain per CR. Processes which produce CRs, such as acceleration or production in hadronic collisions, are incorporated in the source function $s_{i}$, which has dimensions of number per volume per unit energy per unit time.

The main goal of this paper will be to use cluster observations to constrain $F_{cr,i}$, and so we do not restrict the analysis to any specific CR transport model. In §2.2 we consider the self-confinement model of transport for illustrative purposes and to connect the general transport framework to a physical theory, but the analysis presented in §§3 and 4 is intended to be as general as possible.

We next apply our limiting assumptions for CRe and CRp respectively to (3) and derive steady state equations for each. Steady state should be a reasonable approximation for CRe, which have short cooling times of $\approx 100$ Myr for the relevant energies of a few GeV. If the CRe sources in the galaxy cluster do not change faster than this time scale, the CRe density will be held in steady state. Even if the sources fluctuate on faster time scales, the radiative losses would respond to an effective time-smoothing of the CRe source. The resulting CRe density would be a steady state between the radiative cooling and this time-averaged source. The “sources” used below should then be understood to be time-smoothed on 100 Myr scales.

The steady state assumption is less robust for CRp which have much longer cooling times. CRp sources are then primarily offset by transport, which may in principle be very slow. The current CRp density in a region of size $L$ is then determined by the CRp injection history within that region over a CRp crossing time. Assuming the CRp are in steady state would therefore overpredict (underpredict) the real density of CRp if the rate of CRp injection is higher (lower) than the average over this time period. Again, any sources discussed below can be thought of as time-smoothed versions of some real CRp source. Further consideration of time dependent models is beyond the scope of this paper.

As stated above, CRe suffer significant energy losses to synchrotron radiation and inverse Compton (IC) emission. We will therefore approximate advection and diffusion losses as negligible, $\mathbf{F_{cr,e}}=0$, an approximation we will check after the fact in §6. In the steady state, (3) then reduces to

[TABLE]

This has the simple interpretation that, in the steady state, the divergence of any flux of CRe in energy space must be balanced by a CRe source at every point. A CRe distribution function $f_{e}(r,E_{e})$ that suffers energy losses $\dot{E}_{e}<0$ and is replenished by source function $s_{e}(r,E_{e})$ therefore has the steady state solution

[TABLE]

It will be helpful to work in terms of Lorentz factor $\gamma_{e}=E_{e}/m_{e}c^{2}$ instead of $E_{e}$. We define $\gamma_{e}$-dependent functions $f_{e}$ and $s_{e}$ as

[TABLE]

[TABLE]

Equation (5) is then written

[TABLE]

CR protons above about 1 GeV have much longer cooling times than electrons of the same energy (see Enßlin et al. (2011) Figure 2). For protons we therefore approximate the energy losses as negligible, $\dot{E}_{p}=0$ (Adiabatic losses, which are also included in $\dot{E}_{p}$, may be more significant. We discuss these losses for the specific case of CR self-confinement in §2.2 and the consider the effect of ignoring them in our analysis in §6). Equation (3) then reduces, in the steady state, to

[TABLE]

Again this has a simple interpretation that in the steady state, any CRp which leave a given volume must be replenished by the source of CRp in that volume. The solution, assuming the net transport is outward ($F_{cr,p}>0$), is

[TABLE]

To get any further we need to define $F_{cr,p}$, i.e. we need categorize a model of CRp transport. We consider two classes of models here, advection and diffusion. We reiterate that we do not limit analysis to any specific model of CR transport in this work. We use advection and diffusion as two common limiting cases of CR transport that encompass several specific transport models.

We first consider advection, wherein we characterize radial transport by some bulk speed $v_{t}$ which may in principle vary with position and CRp energy:

[TABLE]

As we discuss in §2.2, the self-confinement streaming picture is one example of a CR transport model which has this form. We emphasize that this bulk speed $v_{t}$ is not the speed of any individual CR, which is close to the speed of light, but rather is a measure of the net flow defined by (9).

For a CR flux of this form, the steady state solution is simply

[TABLE]

Evaluating this requires knowing $v_{t}$ and $s_{p}$. This is done in §3.4.1.

We next consider a diffusive model for CR transport, where the CRp flux is

[TABLE]

for some diffusion coefficient $\kappa$. In this regime, the steady state requirement (8) becomes

[TABLE]

We can find $f_{p}(r,E_{p})$ from this by assuming the CRp density at some outermost cluster radius $r_{\textrm{max}}$ is negligible and integrating inward:

[TABLE]

We evaluate this expression in §3.4.1.

2.2 A Specific Treatment

In this section we consider a specific model of CR transport and show how it amounts to one case of the general framework of §2.1. We describe the evolution of CRs in the specific context where transport is mediated by scattering off of Alfvén waves traveling in both directions along the background magnetic field. The self-confinement streaming picture is a special case of this. In the self-confinement model, the majority of Alfvén waves present are generated by CR pressure gradients via the CR streaming instability. CRs will then flow in bulk at a speed which is determined by marginal stability of these waves, and waves will only travel in one direction (see Skilling (1971)).

skilling75a lays out the derivation in detail: starting from the Vlasov equation (2), we first assume the fastest timescale in the problem is the cosmic ray gyromotion in the background magnetic field, and average (2) over gyration angle. We next assume the second fastest timescale is the frequency of scattering from Alfvén waves traveling in both directions along the magnetic field. Denote the scattering rate of CRs off of forward- backward-moving Alfvén waves as $\nu_{\pm}$, which will depend on the CR momentum $p$ and pitch angle cosine $\mu$. skilling75a defines a (momentum-dependent) wave frame

[TABLE]

where $\mathbf{v_{0}}$ is the velocity of the gas, $v_{A}=B/\sqrt{4\pi\rho}$ is the Alfvén speed, and $\mathbf{n}$ is a unit vector pointing in the direction of the background field. The angle brackets indicate an average over $\mu$. Note that in the limit where waves only travel in one direction, as would be expected in the self-confinement regime, this frame reduces to $\mathbf{u}=\mathbf{v_{0}}\pm v_{A}\mathbf{n}$ for all momenta. Strictly, the self confinement model describes transport along magnetic fieldlines, which are not necessarily radial, and predicts its own breakdown in situations where the field aligned components of the thermal gas density gradient and cosmic ray pressure gradient are antiparallel (Skilling, 1971; wiener17a), which is one reason why we consider more general models here.

Shifting the gyration-averaged equation into this frame and averaging over pitch angle, skilling75a arrives at the following equation for $\bar{f}_{i}(\mathbf{r},p)$:

[TABLE]

with spatial and momentum diffusion coefficients defined by

[TABLE]

[TABLE]

Here, $D_{\parallel}$ represents diffusion in space due to pitch-angle scattering by the Alfvén waves, while $D_{\textrm{pp}}$ represents diffusion in energy due to Fermi acceleration.

Note that $\bar{f}_{i}(\mathbf{r},p)$ here is the average of the phase-space density $f_{i}(\mathbf{r},\mathbf{p})$ over all momentum angles:

[TABLE]

This is not the same as the energy distribution function $f_{i}(\mathbf{r},E)$ introduced in equation (3) and use throughout §2.1, but they are related by111The angular spatial coordinates have been restored to $f_{i}$ in this section for ease of comparison to skilling75a.

[TABLE]

With this in mind, we cast the equation of motion (14) specific to the wave scattering model in the form of (3). To do this, note that the second term on the left-hand side of (14) can be written:

[TABLE]

Likewise, the first term on the right-hand side of (14) can be written:

[TABLE]

Equation (14) then becomes, after canceling terms and rearranging:

[TABLE]

or, multiplying both sides by $p^{2}/c$ and further rearranging terms,

[TABLE]

This is now of the same form as (3), but without the source term and with an additional diffusion term in energy space. The two terms under the spatial divergence operator nicely resemble the advective and diffusion forms of transport flux respectively,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

This model of CR transport is one example the generalized treatment of §2.1, where the CR flux has an advective piece and a diffusive piece. The energy gains/losses include a Fermi acceleration term and adiabatic energy changes

[TABLE]

associated with the advective transport. This is the expected adiabatic change for particles of a fluid with adiabatic index $\gamma=4/3$ and is likely to be a loss term due to the decrease in gas density with radius. It is straightforward to add to this specific model additional loss processes, like synchrotron radiation, or source terms.

The comparison is even clearer in the self-confinement limit, where waves only travel in one direction. Then $D_{\textrm{pp}}=0$ and $\mathbf{u}=\mathbf{v_{0}}\pm v_{A}\mathbf{n}$ is independent of $E$, giving:

[TABLE]

which is just equation (3) for the specific case of

[TABLE]

[TABLE]

This reflects the physical description of the self-confinement picture: CRs are carried along with the Alfvén waves, with some amount of diffusion, and suffer adiabatic losses.

3 Properties Determined from Radio Emission

3.1 CRe Profile

As per Rybicki & Lightman (1979), an isotropic, spherically symmetric, power law distribution of electrons

[TABLE]

in the presence of a magnetic field $B(r)$ will emit synchrotron radiation with combined power per unit volume per unit frequency given by222Although we formally integrate the contribution of each electron over all $E_{e}$ to obtain this expression, only a certain range of energies will contribute significantly. Namely, there is an exponential cutoff in the differential power emitted at low $E_{e}$ (the location of this cutoff depends on the frequency $\omega$). This means we don’t have to worry about the fact that our assumed CRe distribution formally contains some electrons that violate the assumption in Rybicki & Lightman (1979) that $\beta\approx 1$. This also means we only need the CRe spectrum to follow a power law in this energy range to be able to use (22).:

[TABLE]

That is, a power law CRe distribution of index $\alpha_{e}$ emits a power law radio spectrum of index $s=(\alpha_{e}-1)/2$. Inverting this relationship, an observed power law radio spectrum of index $s$ implies a power law CRe distribution of index $\alpha_{e}=2s+1$. Let us rewrite (22) in terms of $s$, and also switch to frequency $\nu$ instead of $\omega$ by multiplying by d$\omega$/d$\nu$=2$\pi$.

[TABLE]

[TABLE]

To get the observed radio surface brightness at frequency $\nu$ we would take a line integral of $P_{\textrm{tot}}$ along a line of sight through the cluster. This would involve a non-analytic integral at every point on the sky that we wanted to compare to observations. Let us instead consider the total radio power emitted by the entire cluster. Denote the surface brightness as a function of position on the sky (an observed quantity with units Jy sr*-1*) by $S_{\nu}(\Omega)$. The power per unit volume per frequency bin in radio (an intrinsic quantity with units erg cm*-3* s*-1* Hz*-1*) is $P_{\textrm{tot}}(\nu,r)$ as given above in (23). Then the total observed flux can be expressed in two ways:

[TABLE]

[TABLE]

Here we have assumed the cluster is spherically symmetric, is a distance $D$ away, and has some maximum extent on the sky $\theta_{\textrm{max}}=r_{\textrm{max}}/D$. We have also neglected geometric effects of the size of the cluster by approximating the distance to every point in the cluster as $D$. The left hand side of this equation encases the observations, while the right hand side, via equation (23), contains parts of our cluster model $B(r)$ and $C_{e}(r)$.

Let us plug in (23) for $P_{\textrm{tot}}$ in (24), scaling to the central magnetic field $B_{0}$ and central value of $C_{e}(r)$, denoted $C_{e0}$. We have

[TABLE]

[TABLE]

Thus, for a given model of the magnetic field $B(r)$ and a model for the spatial dependence of $C_{e}(r)$, we can use the observations $S_{\nu}$ to determine the normalization of $C_{e}(r)$.

Let us do this by defining a CRe shape function

[TABLE]

with the inherent normalization $\eta_{e}(0)=1$. Denoting the total integrated radio flux per frequency bin at frequency $\nu$ (the left-hand side of (24)) by $S_{\nu}$ and solving for $C_{e0}$, we arrive at

[TABLE]

where

[TABLE]

is a single constant that encompasses the radio observations and

[TABLE]

is an integral containing information about the cluster model. We motivate some possible forms of $\eta_{e}(r)$ and $B(r)$ in §3.3. As we discuss in §5, our final results are insensitive to our choice of $\eta_{e}(r)$.

3.2 CRe Loss and Source Rate

In the relevant CRe energy range that is probed by the radio observations as described above, energy losses are dominated by synchrotron and inverse Compton (IC) losses. The photon field responsible for the IC losses is itself dominated by cosmic microwave background (CMB) photons. The energy loss rate per electron is then

[TABLE]

where $B_{\textrm{cmb}}\approx 3.24\ \mu$G is the magnetic field with equivalent energy density as the CMB.

Recall from §2 that the steady state equation for an isotropic, spherically-symmetric distribution of CRe, neglecting transport losses, is given by (6):

[TABLE]

For our power law distribution, $f_{e}(r,\gamma_{e})=C_{e0}\eta_{e}(r)\gamma_{e}^{-\alpha_{e}}$. Since $\dot{\gamma}_{e}$ goes as $\gamma_{e}^{2}$, the source function’s energy dependence must be $\gamma_{e}^{1-\alpha_{e}}$. Let us then write the source function in the form $s_{e}(r,\gamma_{e})=G_{e}(r)\gamma_{e}^{1-\alpha_{e}}$. Then

[TABLE]

[TABLE]

Inverting this relationship gives us the source function $G_{e}(r)$ required to give a steady state CRe distribution $C_{e0}\eta_{e}(r)$:

[TABLE]

If we use our formula for $C_{e0}$ based on the radio emission (26), we get

[TABLE]

[TABLE]

Put simply, for a given cluster model, the spatial shape of the source function (assuming the CRe are in steady state) is uniquely determined, and the radio observations determine its normalization.

3.3 CRe Shape Function

Before moving on it will be useful to come up with an educated guess for what the CRe shape function $\eta_{e}(r)$ should be. To do this, consider the more complicated problem of looking at the radio emission as a function of position. This involves a line of sight integral

[TABLE]

with $P_{\textrm{tot}}$ given by equation (23).

Recall that the normalization of $\eta_{e}(r)$ is arbitrary - we are only interested in its spatial dependence. If we just want to find this dependence, we can ignore all factors except those which depend on $r$. We are left with the relationship between the radio brightness $S_{\nu}(\theta)$ and a line integral of the shape function:

[TABLE]

From here on we will model the Coma cluster specifically, but the analysis can be applied more generally. We model the B field of the Coma cluster as one that goes as a power $\alpha_{B}$ of the thermal electron density,

[TABLE]

Let us consider shape functions of the same form:

[TABLE]

going as the electron density to some power $\alpha_{\eta}$. The electron density itself can be modeled with a beta profile:

[TABLE]

with core radius $r_{c}=294$ kpc and $\beta=0.75$ for the specific case of the Coma cluster determined from X-ray surface brightness in Briel et al. (1992). The determination of the shape function thus reduces to determining the power $\alpha_{\eta}$ from

[TABLE]

[TABLE]

The primary dependence on $\theta$ (at least for large $\theta$) appears to be $\theta^{-3\beta(\alpha_{B}(s+1)+\alpha_{\eta})}$. We can compare this to the $\theta$ dependence of $S_{\nu}(\theta)$ at large $\theta$ which is about the same as that of the X-ray surface brightness (very roughly speaking - see Deiss et al. (1997) figure 3), $S_{X}(\theta)\sim\theta^{1-6\beta}$. We end up with

[TABLE]

[TABLE]

For $\beta=0.75$, $\alpha_{B}=0.5$, and $s=1.35$, we get $\alpha_{\eta}=.381$.

Since the extent of the radio halo may depend on frequency and data analysis techniques (Brown & Rudnick (2011) find a more extended halo in the Coma cluster than Deiss et al. (1997) even at the same frequency), we will investigate the effects of changing $\alpha_{\eta}$, using the above formula as a fiducial value. As we explain in §5, while there is some dependence of our results with $\alpha_{\eta}$ (which we vary between 0 and 1), it is subdominant compared to other dependencies. This suggests our results are not sensitive to the assumed shape function $\eta_{e}$.

3.4 Cosmic Ray Proton (CRp) Density

The connection between CRe and cosmic ray protons (CRp) can be complicated, but we can consider two limiting cases. In one limit, CRe secondaries produced by hadronic interactions are negligible in number compared to CRe primaries directly accelerated by the plasma. In the other limit, direct acceleration is negligible, and the CRe are entirely sourced by hadronic reactions. We consider both limits in turn, and check for consistency in §6.

3.4.1 Primary limit

In this limit we simply assume some relative acceleration efficiency of protons to electrons $\zeta\approx 100$,333This value is highly uncertain, so we merely use 100 as a reference value. As discussed in §5.1, the dependence of our results on this parameter will turn out to be linear, so our exact choice for this value is not important. and the CRp source function is just our result from the previous section times $\zeta$:

[TABLE]

[TABLE]

[TABLE]

If we want to find a steady state distribution function $f_{p}(r,E_{p})$, we need to balance the above source function with some loss terms. In §2 we found that the steady state equation for CRp, where we ignore energy losses, is given by (8). We specialized to two classes of transport - advection with some bulk flow $v_{t}(r,E_{p})$ and diffusion with some coefficient $\kappa(r,E_{p})$. For advective transport we derived (10):

[TABLE]

Supposing the functional form of the transport speed $v_{t}$ is known, this, combined with (31) and (36), gives us a way to determine the CRp distribution from the radio emission. Namely,

[TABLE]

[TABLE]

where we have defined another moment of the model profiles

[TABLE]

We see from (37) that the energy dependence of the steady state of $f_{p}$ goes as $E_{p}^{-2s}/v_{t}$, and so depends on the energy dependence of $v_{t}$. For instance, if we assume our transport speed is independent of energy, then we have a power law CRp distribution444We use the tilde notation $\tilde{C_{p}}$ here to highlight the difference in dimensionality between it and $C_{e}$ from the previous section. This arises from defining our CRe distribution function as a function of Lorentz factor $\gamma_{e}$ and our CRp distribution function as a function of energy $E_{p}$. $f_{p}(r,E_{p})=\tilde{C}_{p}(r)E_{p}^{-\alpha_{p}}$ with spectral index $\alpha_{p}=2s=\alpha_{e}-1$ and spatial dependence

[TABLE]

Transport speeds which depend on energy as a power law will have the same form, but with a different exponent $\alpha_{p}$ that is straightforward to determine. More complicated transport speeds require their own treatment on a case by case basis and will not be considered here. Since we only invoke $f_{p}$ here at energies between around 10 and 100 GeV, this analysis covers any transport processes that can be reasonably approximated by a power law over this range.

We next consider diffusion. In §2 we found the steady state equation (12) for CRp undergoing diffusive transport:

[TABLE]

Using what we know for $s_{p}$ and assuming the previously given power law form for $f_{p}$ (which is valid when we assume $\kappa$ has no explicit spatial dependence and is a power law of energy), we get

[TABLE]

[TABLE]

For an energy-independent $\kappa$ this implies, as before, $\alpha_{p}=2s$ and gives us the normalization of $f_{p}$,

[TABLE]

[TABLE]

In the above we have defined a third profile integral

[TABLE]

As with our previous discussion, if $\kappa$ is a non-zero power law in energy, we end up with much the same result, but with a different exponent $\alpha_{p}$.

3.4.2 Secondary limit

In the opposite limit, particle acceleration is negligible and the CRe are entirely sourced by hadronic interactions. Schematically, we start with some distribution of CRp $f_{p}(r,E_{p})$ and some cluster gas density model $n(r)$. From these we determine the pion production rate, and thence the production of CRe via the decay of charged pions.

There are a variety of schemes in the literature for determining pion source functions. The delta-function approximation should work for our purposes - we just approximate that every collision of a single CRp of energy $E_{p}$ with a thermal nucleus results in $\xi$ pions of energy $K_{\pi}T_{p}$ each, where $T_{p}=E_{p}-m_{p}c^{2}$ is the kinetic energy of the incoming proton. In other words, the pion number source function for each collision is

[TABLE]

There are different conventions for choosing the values $K_{\pi}$ and $\xi$ to match experiments. Kelner et al. (2006) choose $K_{\pi}\approx 0.17$ with $\xi=1$ when calculating the neutral pion source only. Pfrommer et al. (2008) use $K_{\pi}=0.25$ and $\xi=2$ for all pion species together. To be unambiguous, we use $\xi$ here to refer to the total pion multiplicity across all pion species and assume that charged and neutral pions are produced approximately in the ratio 2:1.

If we have a CRp distribution $f_{p}(r,E_{p})$, then the rate of collisions (per unit volume) is $R=cn_{N}(r)\sigma_{pp}(E_{p})f_{p}(r,E_{p})$ and so the pion source function is

[TABLE]

[TABLE]

Here, $n_{N}=n_{\textrm{H}}+4n_{\textrm{He}}$ is the target nucleon density and $\sigma_{pp}$ is the (energy-dependent) cross section for proton-proton collisions. The above expression gives the number source function for all pion types, so by the assumed symmetry the charged pion source function is $s_{\pi^{\pm}}=2q_{\pi}/3$.

Let us assume a power law CRp distribution $f_{p}(r,E_{p})=\tilde{C}_{p}(r)E_{p}^{-\alpha_{p}}$. In addition, we will omit the energy dependence of the cross-section $\sigma_{pp}$ with the understanding that it is a very weak function of energy in the relevant energy range around $E_{p}\sim 100$ GeV (see Kelner et al. (2006); Pfrommer et al. (2008) account for this weak dependence with an effective cross section that is a function of the CRp power spectrum, $\sigma_{pp}(\alpha_{p})=32(0.96+e^{4.4-2.4\alpha_{p}})$ mbarn). The charged pion source function is then

[TABLE]

For the next step, we must describe the decay of a charged pion of energy $E_{\pi}$ into electrons and positrons. Following Pfrommer et al. (2008), we use a delta-function approximation and assume that every charged pion of energy $E_{\pi}$ decays into exactly one electron/positron of energy $E_{e}=E_{\pi}/4$. This means

[TABLE]

[TABLE]

[TABLE]

In the last step we have assumed we are looking at high energy CRe with $E_{e}\gg K_{\pi}m_{p}c^{2}/4$. We can also write this in terms of the Lorentz factor $\gamma_{e}$:

[TABLE]

[TABLE]

We know from §3.2 that the CRe source function must be (31) to explain the synchrotron emission. In the limit where CRe are completely sourced by hadronic interactions, we therefore have

[TABLE]

[TABLE]

or

[TABLE]

Recall that in the primary limit, we needed, in addition to the radio observations and a cluster model, a model for the transport (in the form of speed $v_{t}$ or diffusion coefficient $\kappa$) and a relative acceleration efficiency $\zeta$ to determine the spatial distribution of CRp. In the secondary limit, the radio observations and cluster model alone determine the CRp distribution.

4 Properties Determined from $\gamma$-Ray Emission

No diffuse $\gamma$-ray emission has yet been detected from galaxy clusters. At high enough energies, any potential such $\gamma$-rays are expected to be dominated by neutral pion decay. As such, the non-detection of $\gamma$-rays translates to an upper limit on the CRp density. We quantify this below, and follow with a discussion of the meaning of this upper limit in the two limiting cases (primary-dominated and secondary-dominated CRe).

4.1 CRp Density Upper Limit

The source function of $\gamma$-rays coming from neutral pion decay is a straightforward function of the pion source function:

[TABLE]

with $E_{\pi,\textrm{min}}=E_{\gamma}+m_{\pi}^{2}c^{4}/(4E_{\gamma})$. In the previous section we discussed the relationships between the CRp distribution $f_{p}$ and the pion source function $q_{\pi}$, as well as between $q_{\pi}$ and the secondary CRe source function $s_{e}$.

The neutral pion source function is $s_{\pi^{0}}=q_{\pi}/3$. In the previous section we calculated $q_{\pi}$ for a power law CRp distribution $f_{p}(r,E_{p})=\tilde{C}_{p}(r)E_{p}^{-\alpha_{p}}$. Plugging this into the $\gamma$-ray equation above, we have

[TABLE]

[TABLE]

The above integral cannot be done analytically without some further assumptions. If we assume $E_{p}=E_{\pi}/K_{\pi}\gg m_{p}c^{2}$, we can throw away the $m_{p}c^{2}$ in the numerator and the integral can be done as an incomplete Beta function. However, in this limit we also have $E_{\pi}\gg m_{\pi}c^{2}$ so the denominator would also simplify, giving us a simple power law. In this limit we also have $E_{\pi,\textrm{min}}\approx E_{\gamma}$. The integral therefore reduces to

[TABLE]

and so the $\gamma$-ray source function (just from neutral pion decay) is

[TABLE]

and the predicted differential $\gamma$-ray number flux per energy bin at Earth is just the volume integral of this quantity divided by $4\pi D^{2}$. To connect to a potential observation we would then integrate above some energy $E_{\gamma}$:

[TABLE]

[TABLE]

[TABLE]

We can now relate a photon number flux upper limit $F_{\gamma,\textrm{max}}$ to an upper limit on an integral of the CRp distribution:

[TABLE]

[TABLE]

Encompassing the $\gamma$-ray upper limit in the following shorthand

[TABLE]

we can express this constraint as

[TABLE]

This inequality holds for all upper limits on $F_{\gamma}$, so the smallest value of $F_{\gamma,\textrm{max}}(>E_{\gamma})(E_{\gamma}/m_{e}c^{2})^{\alpha_{p}-1}$ across all observations gives us the tightest constraint. Since we have already predicted $\tilde{C}_{p}(r)$ in the primary- and secondary-dominated limits, we can just plug it into the above and interpret the results.

4.1.1 Primary limit

In section 3.4.1 we derived a relationship (39) between the CRp distribution function $\tilde{C}_{p}(r)$ and the synchrotron power $S_{\nu}$ in the limit where CRe are dominated by primaries. By combining this with the above upper limit on $\tilde{C}_{p}(r)$ from $\gamma$-ray observations, again taking $\alpha_{p}=2s$, we obtain:

[TABLE]

[TABLE]

[TABLE]

Ultimately this amounts to a constraint on the transport speed $v_{t}$ for a given model and set of observations. In the very special case of a spatially-independent $v_{t}$ we get a simple lower limit

[TABLE]

We evaluate this limit for different cluster models and discuss its characteristics in section 5. It should be kept in mind that (51) is a radial transport speed and is perhaps most easily interpreted if the magnetic field is relatively coherently oriented. If the field is highly tangled, diffusion is a better model.

Suppose instead we consider diffusive transport characterized by a constant diffusion coefficient $\kappa$ as described by equation (40). Plugging this into the upper limit (49) gives

[TABLE]

[TABLE]

which, finally, can be rearranged into a lower limit on the diffusion coefficient:

[TABLE]

We evaluate this limit for different cluster models and discuss its characteristics in section 5.

4.1.2 Secondary limit

In the limit of secondary CRe dominating the population, we found that the radio observations directly translate into a CRp distribution (45). The $\gamma$-ray upper limit then offers a consistency check - if the inequality above is violated, secondaries cannot dominate the CRe population. Plugging (45) into the upper limit (49), we obtain

[TABLE]

[TABLE]

[TABLE]

We can of course frame this in another way, using the secondary-dominated assumption to make a prediction for the $\gamma$-ray number flux:

[TABLE]

or equivalently, a prediction for the $\gamma$-ray energy flux in some energy band

[TABLE]

to be directly compared with measurements such as from Ackermann et al. (2016). See section 5 for discussion on this limit. Note that since invoking the secondary limit for CRe tells us about the CRp density directly, the above prediction does not depend on any transport parameters, or on $\zeta$, $K_{\pi}$, or $\xi$.

5 Results and Dependencies

We consider here the results of our analysis in the various limiting cases in turn and discuss their dependence on different model parameters.

5.1 Primary Limit, Advective Transport

In the limit of primary CRe dominating over secondaries, we found a lower limit on the advective transport speed given by equation (51). To evaluate this, we model the Coma cluster using the density model $n_{e}$ given in §3.3. The total mass of this model is divergent, so we must also assume a cluster extent $r_{\textrm{max}}$. Based on the X-ray data from Briel et al. (1992) we estimate a maximum radius of $r_{\textrm{max}}=4000$ kpc. This choice is somewhat arbitrary, but as we will see the results do not vary much within the range of reasonable values of $r_{\textrm{max}}$. As a fiducial magnetic field model we take $\alpha_{B}=0.5$ and $B_{0}=3\ \mu$G. We will anchor the radio observations at 1.4 GHz with measured intensity $S_{\textrm{1.4 GHz}}\approx 6.4\times 10^{-24}$ erg cm*-2* s*-1* Hz*-1* (Deiss et al. (1997)), and take the spectral index to be $s=1.35$ (see Thierbach et al. (2003)).

For the $\gamma$-ray flux limits, we impose the energy flux limit in the 7.5 - 10 GeV band from Ackermann et al. (2016) of $\approx 10^{-7}$ MeV cm*-2* s*-1* (the exact value of the limit depends on the spatial model of the emission). We use the highest energy band so as to be most in line with our assumption of high energy. For a power law spectrum of index $2s$, this translates into a photon number flux limit above 10 GeV of

[TABLE]

This will go into the $\mathcal{F}$ quantity in our analysis.

We assume a relative acceleration efficiency of $\zeta=100$, but emphasize that we merely use this as a reference value. This parameter is informed by observations of galactic CRs, but in principle is highly uncertain. The acceleration mechanisms in a galaxy cluster such as Coma may be highly varied, and may include galactic, AGN, and ICM sources. We have purposefully avoided any treatment of specific mechanisms and encompassed all our ignorance of them into $\zeta$555By picking a single value in this way we are essentially assuming the nature of the acceleration is the same at all locations in the cluster. This is hardly certain, but any treatment of acceleration mechanisms beyond this simplification is outside the scope of this work.. Fortunately the dependence of our transport limits $v_{t}$ and $\kappa$ on this parameter are linear, and so any uncertainty in $\zeta$ can be propagated easily. This analysis can also be turned around to impose constraints on $\zeta$ (and consequently on the acceleration mechanisms present) if we fix all other model parameters. We discuss this possibility in §7.

For the above fiducial values, we derive a lower limit on bulk CRp transport of $v_{t}\geq 42$ km s*-1*. The dependence of this limit on some of the model parameters is explicit from equation (51): $v_{t}$ is directly proportional to the proton-to-electron relative acceleration efficiency $\zeta$ and the total radio surface brightness $S_{\nu}$, and inversely proportional to the upper limit on $\gamma$-ray flux $F_{\gamma,\textrm{max}}$. It depends on the normalization of the magnetic field as $B_{0}^{1-s}$ for strong fields ($B_{0}\gg B_{\textrm{cmb}}$) and as $B_{0}^{-1-s}$ for weak fields. If we assume $\xi$ and $K_{\pi}$ are related by the observational constraint that the total fraction of CRp kinetic energy in each collision, $\xi K_{\pi}$, is constant, we see that $v_{t}\propto K_{\pi}^{2s-2}$.

The dependencies on other parameters such as cluster extent $r_{\textrm{max}}$, radio spectral index $s$, and magnetic field-density dependence $\alpha_{B}$ are not explicit, so we try many different models to find empirical trends. The minimum speed $v_{t}$ for various magnetic field models is shown in figure 2. $v_{t}$ seems to be most sensitive to $B_{0}$, roughly following the relation

[TABLE]

A complete table of the results is shown in Table 1. This table shows the effects of changing the model parameters with non-explicit dependence, compared to the fiducial model. We see that changing $B_{0}$ from 1 to 10 decreases the resulting $v_{\textrm{t,min}}$ by a factor of about 100 or more, while changing $\alpha_{B}$ from 0.3 to 0.7 results in an increase of $v_{\textrm{t,min}}$ by about a factor of 10 (see first block of Table 1. In contrast, changing the shape parameter $\alpha_{\eta}$ from 0 (uniform CRe density) to 1 (CRe density proportional to $n_{e}$) decreases $v_{\textrm{t,min}}$ by less than a factor of 2 (see second block of Table 1). This is far less than the dependence on $B$ or on the proton-to-electron relative acceleration efficiency $\zeta$, and suggests that our analysis is not sensitive to our choice of the form of the CRe shape function.

There is only slight dependence of $v_{t}$ with the other parameters in Table 1 except perhaps with the assumed extent of the cluster $r_{\textrm{max}}$ (see fourth block of Table 1). Physically this is because a larger cluster has greater gas content, so CRp densities would have to be lower to avoid a detectable $\gamma$-ray signal. But, as the size of the Coma cluster is almost certainly between 3 and 5 Mpc, the possible variation in $v_{\textrm{t,max}}$ due to uncertainty in this parameters is relatively small.

5.2 Primary Limit, Diffusive Transport

If we describe CRp transport as a diffusive process with diffusion coefficient $\kappa$ that is constant over the (relatively small) observationally constrained range we find the $\gamma$-ray upper limits put a lower limit on $\kappa$ according to (52). For the fiducial values described above, this lower limit comes out to $2.7\times 10^{31}$ cm2/s. We plot the values of $\kappa$ for different models in figure 3. Again the strongest dependence is on $B_{0}$, following the same scaling as $v_{t}$:

[TABLE]

The results of changing the other parameters are again shown in Table 1. As with $v_{\textrm{t,min}}$ the dependencies of $\kappa$ on parameters other than the magnetic field are significantly subdominant.

We can convert these limits on diffusion coefficient $\kappa$ into limits on the mean free path due to scattering seen by individual CRs $l_{\textrm{mfp}}$ by the simple relation $l_{\textrm{mfp}}=3\kappa/c$. For our fiducial model, $l_{\textrm{mfp}}$ must be at least 900 pc for CRp to be able to leave the cluster fast enough to bring $\gamma$-radiation under the detection limit in Coma.

In a model where radial transport is due to scattering by magnetic inhomogeneities at the cosmic ray gyroscale, and the background magnetic fieldlines are nearly radial and straight, the advection speed $v_{t}$ and diffusion coefficient $\kappa$ can be combined to produce a lengthscale $R\equiv\kappa/v_{t}$ which is representative of the size of the system; for the values given here, $R\sim$ 2 Mpc, which is reasonable for Coma (the general expression for this quantity is

[TABLE]

where $\mathcal{I}_{2}$ and $\mathcal{I}_{3}$ are defined in equations (38) and (41)).

5.3 Magnetic Field Tangling

Although our aim in this work is to put limits on the transport of CRs in a general framework, it is useful to consider these limits in the context of specific transport models and see what they imply. For example, we could attribute the transport described above to turbulent motions of the gas in the Coma cluster. If we describe the turbulence by a length scale $L_{\textrm{MHD}}$, where the characteristic speed is $v_{A}$, the eddy motions will induce diffusion of order

[TABLE]

However, for the typical Alfvén speeds of $\sim 100$ km s*-1* this would require $L_{\textrm{MHD}}$ to be as large as the Coma cluster itself. For turbulence this strong, wave damping effects would also be important, and it is likely that streaming motions would outpace any diffusion from turbulence (see Wiener et al. (2013)). Put another way, the diffusivities estimated in Table 1 are much larger than the eddy diffusivity corresponding to large scale cluster turbulence unless the turbulent length scale is of order the global size of the cluster.

Alternatively, motivated by observations which show that cluster magnetic fields have strong random components (Vogt & Enßlin (2005), Bonafede et al. (2010)) we can interpret $\kappa$ as the result of the random walk of CRs along a magnetic field which is tangled on a lengthscale $l_{\textrm{corr}}$ (this works if $l_{\textrm{corr}}$ is intermediate between the size of the system and the cosmic ray gyroradius). In this case, the maximum CR diffusion rate is the one given by field line random walk at the speed of light (Minnie et al. (2009)): $\kappa\sim l_{\textrm{corr}}c/3$. This is also approximately the maximum rate of transport across magnetic fieldlines found in Desiati & Zweibel (2014). Additional scattering up and down the tangled fieldlines, which could be represented by a streaming velocity such as $v_{s}<c$, only reduces the transport rate (Rechester & Rosenbluth (1978)). Therefore, the minimum tangling length $l_{\textrm{corr}}$ consistent with the diffuse $\gamma$-ray upper limits is the same as the mean free path $l_{\textrm{mfp}}$ estimated in the preceding section.

To drive this point home, we have predicted what the streaming speed $v_{s}$ would be as a function of $r$ for the 100 GeV CRp that would be responsible for (undetected) $\gamma$-radiation for a specific self-confinement model. We calculate $v_{s}$ by taking the resulting CRp density profile from our model (40) and plugging it into the formalism in Wiener et al. (2013) and Wiener et al. (2018) for cosmic ray self-confinement with wave damping due to MHD turbulence, characterized by the length scale $L_{\textrm{MHD}}=100$ kpc. According to this model, the bulk streaming speed $v_{s}$ of CRs along magnetic field lines is set by equating the growth rate of the CR streaming instability to the rate of damping by the background medium. This gives us $v_{s}$ as a function of position in the cluster, which can be translated into a minimum allowed field tangling length $l_{\textrm{corr}}=3\kappa/v_{s}$, where $\kappa$ is the minimum diffusion coefficient derived from the observations according to our analysis in this work. We record the maximum (most restrictive) $l_{\textrm{corr}}$ over all values of $r$ for each model and plot the results in figure 4, which also includes the range of length scales in the magnetic turbulence power spectrum inferred from Faraday rotation observations of Coma (Bonafede et al. (2010)).

We see from figure 4 that the magnetic field observed in Coma is tangled too much ($l_{\textrm{corr}}$ is too small) to explain the lack of $\gamma$-radiation under this set of assumptions (primary CRe dominate, CRp are transported via streaming with the assumed level of wave damping). This tension is relieved if the cluster magnetic field is in the upper range of values tested, which reduces the CRe population required to explain the radio emission, or if the relative proton-to-electron acceleration efficiency factor $\zeta$ is smaller than Milky Way values.

5.4 Secondary Limit

In the limit where the CRe are dominated by secondaries from hadronic interactions, we derived a flux prediction (55) which can be used to predict $\gamma$-ray energy fluxes in arbitrary energy bands. We compare with the highest energy band used in Ackermann et al. (2016), which is about 7.5 - 10 GeV. Ackermann et al. (2016) find an upper limit in this bin of (depending on their emission model) around 6-10$\times 10^{-11}$ GeV cm*-2* s*-1*. We compare this to the prediction of our models $\Phi_{\gamma}$ in figure 5. Variations of $\Phi_{\gamma}$ with all the non-explicit parameters are shown in the last column of Table 1. Our fiducial model predicts a $\gamma$-ray flux of $1.9\times 10^{-9}$ GeV cm*-2* s*-1*, far above the Fermi upper limit. A secondary-only model for Coma’s radio halo is thus highly disfavored.

6 Consistency Checks

In this section we check some of the assumptions we have made in various regimes.

6.1 CRe Transport Losses

Throughout this work we have assumed that the CRe losses in the relevant energy range were dominated by synchrotron and IC emission. Namely, we neglected the transport term the CRe evolution equation (second term in equation (3)). We estimate these transport effects here and compare them to the radiative losses.

Suppose we can characterize the CRe transport with a bulk transport speed $v_{t}(r,\gamma_{e})$, as we did for the CRp. Then the number loss (in particles per volume per unit $\gamma_{e}$ per time) due to transport is:

[TABLE]

[TABLE]

This is to be compared to the radiative loss rate (the third term in equation (3)),

[TABLE]

[TABLE]

In the above we have assumed a power law for the CRe distribution, $\dot{\gamma}_{\textrm{loss}}$ refers to the radiative losses per electron (29), and $G_{e}(r)$ refers specifically to the solution from our previous analysis (30,31).

In the simple case of a constant, energy-independent transport speed $v_{t}(r,\gamma_{e})=v_{t}$, the ratio of transport losses to radiative losses for CRe is

[TABLE]

Using our estimate of the CRe shape function $\eta_{e}(r)=n_{e}^{\alpha_{\eta}}(r)=(1+(r/r_{c})^{2})^{-3\alpha_{\eta}\beta/2}$, we get

[TABLE]

and so the transport to radiative loss ratio is

[TABLE]

[TABLE]

[TABLE]

The last term is a fairly weak function of $r$ of order unity, so let us discard it. Picking some scale values, the rest of the expression comes to

[TABLE]

Recall that $r_{c}=297$ kpc for the Coma cluster. The scale field 4.4 $\mu$G is chosen as the effective combined field of a 3 $\mu$G field with $B_{\textrm{cmb}}=3.24\ \mu$G.

We see that our assumption that this quantity is much much less than one is valid except when the transport speed $v_{t}$ is very large, or $r$ is very small. The approximation is better for higher energies.

6.2 CRp Adiabatic Losses

As mentioned in §2.1, energy losses to processes like Coulomb scattering and hadronic collisions are small for CRp above 1 GeV. But adiabatic losses may be significant. It is possible for CRs to be transported without suffering adiabatic losses, such as in a free-streaming model where the magnetic field is completely static. But if CR transport is governed by wave scattering from self-generated turbulence, as in §2.2, we would expect significant adiabatic losses to be present which have not been accounted for in our analysis.

A detailed treatment describing how adiabatic loss times compare with transport times in the most general sense is beyond the scope of this paper. Instead we consider our limiting cases of advection with constant outward velocity $v_{t}$ and diffusion with constant coefficient $\kappa$. If the transport is diffusion-dominated, then there are no associated adiabatic losses and our analysis is correct as is.

Consider instead the advective limit for the specific case of constant $v_{t}$. If we suppose adiabatic losses are of the same order as expected from the specific model in §2.2, then we have, using equation (20) with $\mathbf{u}=v_{t}\hat{r}$,

[TABLE]

The adiabatic loss time $t_{\textrm{adia}}=E/\dot{E}_{\textrm{adia}}$ is then

[TABLE]

While loss times are larger than transport times $r/v_{t}$, they are of the same order. Our derivation of equation (10), which inherently included the assumption that adiabatic losses were negligible, would not be valid in this case (it would still be valid for a transport model with no adiabatic losses).

While we do not calculate the quantitative effect adiabatic losses on our results, the qualitative result of including these losses is clear. The total loss rate of energy in CRp is dictated by the observations. In our analysis we attribute all of this loss to transport. If some energy is lost through other processes, the transport need not be as fast as found in §5 to be consistent with observations. The limits derived in §5 would be relaxed, i.e. the values of $v_{t,\textrm{min}}$ would be lower by some amount.

6.3 Estimate of CRe Secondaries

We have done our analysis in two limiting cases, where the CRe are dominated by primaries and secondaries respectively. We have already seen in §4.1.2 that the upper limits on $\gamma$-ray flux rule out the secondary-dominated limit in Coma unless the magnetic field is fairly high. However, we have yet to check how many secondaries would be produced in the primary-dominated limit. If too many are produced, the primary-dominated limit assumption is inconsistent.

The CRp distribution in the primary-dominated limit inferred from radio observations was determined to be (39). The secondary CRe source function expected from a general CRp distribution was determined to be (44). The expected secondary CRe source function in the primary-dominated limit is then

[TABLE]

[TABLE]

This must be compared with the total CRe source function given by (31). Their ratio is

[TABLE]

This can be found numerically, but we can get a good estimate by recognizing that the quantity in the braces is of order $r$:

[TABLE]

For $s$ close to 1, and a fully ionized medium, $(8/3)(\xi/K_{\pi})(K_{\pi}/4)^{2s}\zeta\mu_{e}$ is approximately 10, and $\sigma_{pp}$ is approximately 50 mbarn. Picking scale values for $n_{0}$ and $r$, we have

[TABLE]

Considering how large we expect $c/v_{t}$ to be, this suggests that our assumption of primary CRe domination is inconsistent. Even at the large end of the range of values of $v_{t}$ we found, around 100 km s*-1*, we may expect some 50% of CRe to be secondaries. For expected relative acceleration efficiencies and CR bulk transport speeds, secondary CRe injected from hadronic processes must make up a significant fraction of the total CRe source.

What does this mean for our earlier analysis? Qualitatively, the presence of secondary CRe suggests that for the same level of radio emission, fewer primary CRe are necessary. That is, we can explain the radio emission with some primary CRe source function that is some fraction of the one we derived in §3.2. This makes it easier to get under the $\gamma$-ray upper limits, implying we can get away with lower transport speeds $v_{t}$. A more quantitative treatment is beyond the scope of this work.

We can alternatively derive a new, more restrictive lower limit on the transport speed $v_{t}$ such that secondaries are negligible. Let’s say we required that secondaries make up less than one percent of all CRe. Then the transport speed must be at least

[TABLE]

Such speeds are not expected and we can reasonably surmise that secondary CRe in the Coma cluster are non-negligible. We touch on the implications of this in the next section.

7 Discussion and Conclusions

Diffuse non-thermal radio continuum emission detected in galaxy clusters shows that cosmic ray electrons are present, but, contrary to expectations, there is no evidence that cosmic ray protons are present. The purpose of this paper was to explore the implications of these observations and find limits for cosmic ray acceleration and propagation in galaxy clusters.

We presented an analysis of CR populations based on radio and $\gamma$-ray observations. Given a cluster thermal plasma density and magnetic field strength model, radio observations constrain the CRe population and, assuming steady state, the CRe source strength. $\gamma$-ray upper limits put upper limits on the CRp density present in the cluster. We then combined these results in two limiting cases. In the case where we assume secondary CRe production is negligible, we arrived at a minimum transport speed necessary to reduce the CRp density below the upper limits. In the opposite case where we assume primary CRe production is negligible, we obtain a parameter-dependent $\gamma$-ray flux prediction to be compared to the upper limits. These equations are general in that they can be applied to any cluster for which the assumptions of the model are reasonably well satisfied.

Using radio (Thierbach et al. (2003)) and $\gamma$-ray (Ackermann et al. (2016)) observations of the Coma cluster we have found that under the assumption that all synchrotron-emitting CRe are directly accelerated, we require CRp to be transported outward at a bulk speed of 10-100 km s*-1* or with diffusion coefficient 1031-1032 cm2 s*-1* in order to explain the lack of $\gamma$-ray detection for our fiducial set of cluster parameters. These speeds are sub-Alfvénic, so a self-confinement model with a magnetic field which is coherent on large scales is consistent with the radio and $\gamma$-ray observations. If the field is tangled on the scales of a few kpc, as suggested by Bonafede et al. (2010), then CRs which free stream at the speed of light will experience diffusive transport with diffusion coefficient comparable to the values found here. So a free streaming model with a tangled field is also consistent with the observations for our choice of parameters.

Conversely, a self-confinement model with a tangled field is not consistent with observations for our chosen parameters, as CRp would not be transported outward fast enough. To bring this model in line with observations requires changing our assumed parameters. In section 5 we presented the dependence of our results on these parameters, shown in detail in Table 1. We can use observations of the magnetic field to constrain other parameters in our model, namely the relative acceleration efficiency $\zeta$. As mentioned, some studies suggest that the confinement time of CRs in galaxy clusters is of order a Hubble time. This is not incompatible with our steady state analysis here, as it does not imply that no CRs ever leave the cluster. However, with reduced outward transport we require CRp to be accelerated less efficiently than in our fiducial model in order to be compatible with both the radio and $\gamma$-ray observations. We find that our model Coma cluster, with field tangled at kpc scales (Bonafede et al. (2010)) and CR transport determined by self-confinement, is consistent with observations if the relative proton to electron acceleration efficiency is less than 15. The exact value of this limit depends on the level of assumed wave damping. Various regions of the (extensive) parameter space can be ruled out in this way.

Additionally, at the comparatively low transport speeds derived here we expect CRe secondaries to make a significant contribution to the population. We also found that unless the magnetic field in the center of the Coma cluster is around 10 $\mu$G or above, much higher than observed, the CRe population cannot be dominated by secondaries or else the $\gamma$-ray flux prediction exceeds upper limits. This suggests neither component, primaries or secondaries, can be consistently ignored. The observations of Coma can only be explained by a hybrid population of primary CRe which are being continuously injected into the cluster and secondary CRe which arise from hadronic interactions. A more detailed treatment is therefore necessary to examine CR transport in the Coma cluster. We leave such a treatment for future work, but speculate that the inclusion of a small secondary component in the primary-dominated analysis would alleviate the limits derived here: the presence of secondaries implies less acceleration of primary CRe is needed to produce the observed radio emission, which in turn implies CRp are produced at a slower rate.

Finally, we reiterate that while we have specified our results to the example of the Coma cluster, the analysis is general and can be used on any individual system reasonably approximated by spherical symmetry.

Acknowledgements: We are happy to acknowledge discussions with Gianfranco Brunetti. Comments by anonymous referees helped improve the content and clarity of the presentation. JW and EGZ acknowledge support by NSF Grant AST-1616037, the WARF Foundation, and the Vilas Trust.