ACCELERATION OF THERMAL PROTONS BY GENERIC PHENOMENOLOGICAL MECHANISMS

Most particle acceleration mechanisms invoked for astrophysical sources must start with the acceleration of low-energy particles of the background magnetized plasma with a thermal or Maxwellian distribution confined in a finite volume with density n, temperature T, and magnetic field B. This initial phase of acceleration can produce a distribution consisting of a quasi-thermal component plus a nonthermal tail. Such distributions can be approximated by the kappa distribution (see, e.g., Pierrad & Lazar 2010) and can be the seeds in the so-called thermal leakage injection model in diffusive shock acceleration (see, e.g., Kang et al. 2014). The main goal of this and an earlier paper (Petrosian & East 2008, hereafter PE08) is to explore the possibility of producing a prominent nonthermal tail by some generic phenomenological acceleration process. PE08 treated the generation of a nonthermal tail in the distribution of electrons neutralized by a cold noninteracting proton population. Here we evaluate the conditions required for the generation of nonthermal proton spectra subject to a similar energizing mechanism using the coupled Fokker–Planck (FP) kinetic equations.

The nonthermal tails in kappa-like electron distributions could also be responsible for nonthermal emission in the hard X-ray regime for hot plasmas ($T\sim {10}^{6}-{10}^{8}$ K) such as in solar flares, as demonstrated in Hamilton & Petrosian (1992), and possibly in the hot intracluster medium (ICM) of some clusters of galaxies. Many clusters show a significant signature of nonthermal activity first observed as radio radiation that is due to synchrotron emission by a population of relativistic (Lorentz factor $\gamma \sim {10}^{4}$) electrons in a $B\sim \mu {\rm{G}}$ field. In addition, earlier observations by several hard X-ray instruments indicated that, in addition to the well known thermal soft X-ray radiation, some clusters show excess radiation above tens of keV that could not be fitted by a single-temperature thermal bremsstrahlung model. For a review of early observations see Durret et al. (2008), Rephaeli et al. (2008), and Ferrari et al. (2008). Initially the excess radiation was assumed to arise from bremsstrahlung of a nonthermal tail of electrons extending to ∼100 keV. However, as stressed by Petrosian (2001) and later rigorously proved by PE08, this scenario would cause excessive heating of the ICM. Subsequently it was suggested by Wolfe & Melia (2006) that a nonthermal tail in proton distributions extending to hundreds of MeV could also be a source of hard X-rays (10–100 keV) via inverse nonthermal bremsstrahlung produced by their scattering of lower velocity thermal electrons. This process was also proposed as a possible source of hard X-rays in solar flares (see review by Emslie & Brown 1985) because of the possible higher X-ray yield compared to electrons. Exploring this possibility is the subject of our paper. As we will show, there are similar difficulties with this scenario as well.

It should be noted that in the ICM hard X-rays could also be produced by the inverse Compton (IC) scattering of cosmic microwave background photons by the radio-producing relativistic electrons. The difficulty with this model is that it requires a lower magnetic field than indicated by Faraday rotation observations (see, e.g., Feretti et al. 1995; Clarke et al. 2001) or that expected from equipartition. Reviews of these emission processes and possible acceleration mechanisms are given in Petrosian et al. (2008) and Petrosian & Bykov (2008).³ However, more recent observations have cast doubt on the reality of the claimed hard X-ray excesses. In particular, observation of Coma by Suzaku (Wik et al. 2009) and NuStar (Wik et al. 2014), and the Bullet cluster (Astaldello et al. 2015) give only upper limits consistent with the claimed relatively high magnetic fields from Faraday rotation measurements. The non-detection of many clusters by Fermi-LAT (Ackermann et al. 2010) supports this view and puts stringent constraints on population of high-energy cosmic rays that may be produced in the shocks arising from mergers during large-scale structure formations (see the recent review by Brunetti & Jones 2014).

The crucial feature in the production of nonthermal electron or proton tails is the interplay between the particle–particle (and particle–external fields) interactions that causes energy loss and momentum diffusion, and the acceleration or energizing mechanism due to interaction of particles with plasma turbulence and/or converging flows (e.g., shocks). For nonrelativistic background plasma $({T}\lt {10}^{9}\;{\rm{K}})$ the energy loss is dominated by Coulomb collisions, which tend to equilibrate the distribution of electrons and protons to a single-temperature Maxwellian. Acceleration rates lower than the Coulomb rates (i.e., with longer timescales) tend to only heat the plasma and higher acceleration rates lead to a runaway distribution. Thus, to produce a significant but not excessive nonthermal population we need acceleration rates comparable to the Coulomb rates, so that in PE08, dealing with acceleration of electrons only, we used acceleration timescales comparable to electron–electron (e–e) Coulomb times. Since the latter are much shorter than electron–proton (e–p) (and proton–proton (p–p)) collision times at all energies, protons remained decoupled and kept their initial distribution. As shown in PE08, production of nonthermal electron tails required special conditions, at least for ICM conditions where one needs to maintain a relatively constant temperature of ${kT}\sim 1-10\;\mathrm{keV}$ over a Hubble time. This is not an issue for short-lived solar flares (Hamilton & Petrosian 1992).

The situation is more complicated for protons because the longer p–p collisional time requires a longer acceleration time. Such a rate for electrons would lead to their heating. Therefore, we assume no acceleration, or a much lower rate of acceleration, of electrons. This can come about for pure Alfvénic turbulence which does not interact with low-energy electrons. (Low-energy electrons interact with the whistler waves of the fast mode branch.) However, this does not mean that electrons can be decoupled from protons. This is because proton–electron (p–e) and p–p collision times can be comparable, so that some of the energy gained by protons can be transferred to electrons. Therefore, we must treat the coupled electron and proton kinetic equations simultaneously.

In addition, the energy dependence of the acceleration rate also plays an important role, with acceleration rates that increase with energy being more efficient in producing nonthermal tails.

In Section 2, we derive and present coefficients for Coulomb interactions and stochastic acceleration (SA) to be used in a FP kinetic equation, and explain our algorithm for solving the coupled FP equations of protons and electrons. In Section 3, we apply the algorithm to the SA of thermal background particles, and examine the time evolution of the proton and electron distributions for various acceleration model parameters. In Section 4, we summarize our results and discuss their implications. Some details of the acceleration and a new formula for the proton inelastic energy loss rate due to pion production are presented in the appendixes.

In order to compute the time evolution of the proton and electron distributions when the two species of particles interact not only internally but also mutually via Coulomb collisions, we need to solve their coupled kinetic equations. We use the FP kinetic equation with the following simplifying assumptions. We assume a magnetized background plasma where particles are tied to magnetic field lines and undergo momentum and pitch angle diffusion. We also assume that the mean free path for scattering is much smaller than the size of the region, which means that the particle momentum distribution is isotropic. We further assume that the system is closed, that is, there is no escape or injection of particles. In this case the pitch angle-averaged, spatially integrated momentum distribution $f(t,p)$ satisfies the simple equation $\partial f/\partial t=(1/{p}^{2})$ $\partial [{p}^{2}{D}_{\mathrm{pp}}(\partial f/\partial p)]/\partial p.$ For convenience we use the energy distribution defined as $N(t,E){dE}=4\pi {p}^{2}f(t,p){dp}$ to obtain the following two commonly used forms of the kinetic equation:

$$\begin{eqnarray}\displaystyle \frac{\partial N}{\partial t} & = & \displaystyle \frac{\partial }{\partial E}\left({D}_{\mathrm{EE}}\displaystyle \frac{\partial N}{\partial E}\right)-\displaystyle \frac{\partial [(A(E)N]}{\partial E}\\ & = & \displaystyle \frac{{\partial }^{2}({D}_{\mathrm{EE}}N)}{\partial {E}^{2}}-\displaystyle \frac{\partial [\tilde{A}(E)N]}{\partial E},\end{eqnarray} \tag{ 1 }$$

where the energy diffusion coefficient and the direct energy gain (or loss) rates are given as

$$\begin{eqnarray}&&{D}_{\mathrm{EE}}={\beta }^{2}{D}_{\mathrm{pp}}\;\mathrm{and}\;A(E)=\zeta {D}_{\mathrm{EE}}/E\;\mathrm{with}\;\zeta =\displaystyle \frac{2{\gamma }^{2}-1}{{\gamma }^{2}+\gamma }.\end{eqnarray} \tag{ 2 }$$

Here $\beta =v/c$ is the particle velocity v in units of the speed of light c and the Lorentz factor $\gamma ={(1-{\beta }^{2})}^{-1/2}.$ As shown in Petrosian & Chen (2014; see also Tsytovich 1977) $\tilde{A}(E)\equiv A(E)+{{dD}}_{\mathrm{EE}}/{dE}$ provides a more accurate representation of energy gain (or Coulomb energy loss) rate than $A(E).$

In what follows we will present results from numerical solutions of the first form of the above FP equation including the effects of SA by turbulence and Coulomb collision. We will assume a background plasma of constant density n and field B (but not temperature T), and a constant energy level and spectrum of turbulence so that the transport coefficients D_EE and A are constant in time. Following PE08 we use the simple form

$$\begin{eqnarray}&&{A}_{\mathrm{SA}}(E)=E/[{\tau }_{0}{(1+{E}_{c}/E)}^{q}],\end{eqnarray} \tag{ 3 }$$

or the characteristic acceleration timescale

$$\begin{eqnarray}&&{\tau }_{{ac}}\equiv E/A(E)={\tau }_{0}{(1+{E}_{c}/E)}^{q},\end{eqnarray} \tag{ 4 }$$

described by the three parameters q, ${\tau }_{0},$ and E_c. The thin black (solid and dotted) curves in Figure 1 show some examples of acceleration times used in the next section.

Zoom In Zoom Out Reset image size

As described above, ${\tau }_{0},$ which is a measure of the acceleration time, should be comparable to the relevant collision times and E_c allows us to introduce a break in the energy dependence of the acceleration rate (making it steeper or flatter at higher energies depending on the value of the index q). The break will be important only for low values of E_c and will have little effect for ${E}_{c}\gg {kT}.$ For SA by turbulence the index q is related to the spectral index q' of turbulence. For relativistic energies and Alfvénic turbulence $q=q^{\prime} -2.$ In the inertial range one expects a Kolmogorov or Kraichnin spectrum with $q^{\prime} =5/3$ or 3/2 so that $q=-1/3$ and $-1/2,$ respectively. But the turbulence spectrum is expected to be steeper in the damping range at smaller scales (large wavevectors) where q may become positive. At lower energies this relation becomes more complicated as interactions with many other modes become important (see, e.g., Pryadko & Petrosian 1997; Petrosian & Liu 2004). In what follows we will use $q=1,0,-1$ to include all possible cases.⁴

2.1. Coulomb Collisions

In the cold-target approximation, in which a test particle of charge Z_ie, mass m_i, and velocity v_i collides with stationary target particles of charge Z_je, mass m_j, and number density n_j, the Coulomb energy diffusion rate is negligible but the Coulomb energy loss rate is

$$\begin{eqnarray}&&{\dot{E}}_{{ij}}^{\mathrm{cold}}\equiv -{\tilde{A}}_{\mathrm{Coul}}=4\pi \mathrm{ln}{\rm{\Lambda }}{Z}_{i}^{2}{Z}_{j}^{2}{e}^{4}{n}_{j}/({m}_{j}{v}_{i}),\end{eqnarray} \tag{ 5 }$$

where $\mathrm{ln}{\rm{\Lambda }}\sim$ 20–40 is the Coulomb logarithm. In what follows we will be interested in electron and proton test particles with Z_i = 1. In most astrophysical situations we are dealing with fully ionized ions so the primary target particles are electrons and protons so Z_j = 1 also. Contributions to Coulomb losses from other ions are negligible except those from alpha particles with a relative number density ${n}_{\alpha }/{n}_{{\rm{p}}}\simeq 0.08$ and ${Z}_{j}=2,$ which may contribute up to 20%. Note that this also means that ${n}_{{\rm{p}}}/{n}_{{\rm{e}}}=0.86$ and ${n}_{\alpha }/{n}_{{\rm{e}}}=0.07.$ If we define ${\tau }_{\mathrm{Coul}}={(4\pi {r}_{0}^{2}{{cn}}_{{\rm{e}}}\mathrm{ln}{\rm{\Lambda }})}^{-1},$ where r₀ is the classical electron radius, then

$$\begin{eqnarray}&&{\dot{E}}_{{ij}}^{\mathrm{cold}}=\left(\displaystyle \frac{{m}_{{\rm{e}}}{c}^{2}}{{\tau }_{\mathrm{Coul}}\beta }\right)\left(\displaystyle \frac{{m}_{{\rm{e}}}}{{m}_{j}}\right)\left(\displaystyle \frac{{Z}_{j}^{2}{n}_{j}}{{n}_{{\rm{e}}}}\right).\end{eqnarray} \tag{ 6 }$$

The cold-target approximation can be used for test particle energies ${E}_{i}\gg {E}_{j}$ of the background plasma. For interaction with a thermal background plasma with temperature T, when the test particle energy E_i approaches ${E}_{j}\sim {kT},$ we must use the hot plasma energy exchange rates, which yield a finite energy diffusion rate and can be an energy gain when ${E}_{i}\lt {kT}.$ In PE08 we used hot plasma rate equations for nonrelativistic electrons given, e.g., by Miller et al. (1996) and Nayakshin & Melia (1998). Here we are interested in interactions of both electrons and protons. This requires generalization of the two functions $G(E,E^{\prime} )$ and $H(E,E^{\prime} )$ of PE08 that describe the loss and diffusion rates, respectively, as follows.

Generalization of Equation (7) in PE08 gives

$$\begin{eqnarray}{G}_{{ij}}({E}_{i},{E}_{j})=\left\{\begin{array}{ll}{-{\beta }_{j}}^{-1}({m}_{{\rm{e}}}/{m}_{i}) & \mathrm{for}\ {\beta }_{i}\lt {\beta }_{j},\\ 0 & \mathrm{for}\ {\beta }_{i}={\beta }_{j},\\ {{\beta }_{i}}^{-1}({m}_{{\rm{e}}}/{m}_{j}) & \mathrm{for}\ {\beta }_{i}\gt {\beta }_{j}.\end{array}\right.\end{eqnarray} \tag{ 7 }$$

The meaning of G_ij is clear; if the test particle is faster than the target particle, it will lose energy to the target particle according to the cold Coulomb loss rate. Conversely, if it is slower than the target particle, it will gain energy from the target particle. From this we can calculate the energy exchange rate (positive for loss, negative for gain) ${\dot{E}}_{{ij}}$ by integrating G_ij over the distribution ${N}_{j}({E}_{j})$ of the target particles:

$$\begin{eqnarray}&&{\dot{E}}_{{ij}}({E}_{i})=\displaystyle \frac{{m}_{{\rm{e}}}{c}^{2}}{{\tau }_{\mathrm{Coul}}}{\displaystyle \int }_{0}^{\infty }{G}_{{ij}}({E}_{i},{E}_{j}){N}_{j}({E}_{j}){{dE}}_{j}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{{m}_{{\rm{e}}}{c}^{2}}{{\tau }_{\mathrm{Coul}}{\beta }_{i}}\left({\displaystyle \int }_{0}^{\frac{{m}_{j}}{{m}_{i}}E}{N}_{j}({E}_{j}){{dE}}_{j}-{\displaystyle \int }_{\frac{{m}_{j}}{{m}_{i}}E}^{\infty }\displaystyle \frac{{\beta }_{i}}{{\beta }_{j}}{N}_{j}({E}_{j}){{dE}}_{j}\right).\end{eqnarray} \tag{ 9 }$$

In a similar manner, we derive the corresponding Coulomb diffusion coefficient D_ij. Generalizing $H(E,E^{\prime} )$ given by Equation (10) of PE08, we define ${H}_{{ij}}({E}_{i},{E}_{j})$ by

$$\begin{eqnarray}{H}_{{ij}}({E}_{i},{E}_{j})=\left\{\begin{array}{ll}{{\beta }_{i}}^{2}/(3{\beta }_{j}) & \mathrm{for}\ {\beta }_{i}\lt {\beta }_{j},\\ {\beta }_{i}/3 & \mathrm{for}\ {\beta }_{i}={\beta }_{j},\\ {{\beta }_{j}}^{2}/(3{\beta }_{i}) & \mathrm{for}\ {\beta }_{i}\gt {\beta }_{j}.\end{array}\right.\end{eqnarray} \tag{ 10 }$$

D_ij is given by integrating H_ij over the target particle distribution:

$$\begin{eqnarray}&&{D}_{{ij}}({E}_{i})=\displaystyle \frac{{({m}_{{\rm{e}}}{c}^{2})}^{2}}{{\tau }_{\mathrm{Coul}}}{\displaystyle \int }_{0}^{\infty }{H}_{{ij}}({E}_{i},{E}_{j}){N}_{j}({E}_{j}){{dE}}_{j}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&=\;\displaystyle \frac{{({m}_{{\rm{e}}}{c}^{2})}^{2}}{{\tau }_{\mathrm{Coul}}{\beta }_{i}}\left({\displaystyle \int }_{0}^{\frac{{m}_{j}}{{m}_{i}}E}\displaystyle \frac{{\beta }_{j}^{2}}{3}{N}_{j}({E}_{j}){{dE}}_{j}+{\displaystyle \int }_{\frac{{m}_{j}}{{m}_{i}}E}^{\infty }\displaystyle \frac{{\beta }_{i}^{3}}{3{\beta }_{j}}{N}_{j}({E}_{j}){{dE}}_{j}\right).\end{eqnarray} \tag{ 12 }$$

Note that, while G_ij is inversely proportional to the mass of the slower particle, H_ij has no explicit dependence on the mass of any particle, and in particular for identical particles ${H}_{{ii}}={\beta }_{i}/3$ is independent of target particle mass and velocity.

When the number of particles in the nonthermal component of ${N}_{j}({E}_{j})$ is much less than that in the thermal component we can approximate the distribution of the target particles by the nonrelativistic Maxwellian form

$$\begin{eqnarray}&&{N}_{T}(E)=\displaystyle \frac{2}{\sqrt{\pi }}{({kT})}^{-3/2}\sqrt{E}{e}^{-E/{kT}},\end{eqnarray} \tag{ 13 }$$

which then yields the hot Coulomb coefficients

$$\begin{eqnarray}{\dot{E}}_{{ij}}^{\mathrm{hot}}(E) & \equiv & -{\tilde{A}}_{\mathrm{Coul}}^{\mathrm{hot}}=\displaystyle \frac{{m}_{{\rm{e}}}{c}^{2}}{{\tau }_{\mathrm{Coul}}\beta }\left(\displaystyle \frac{{m}_{{\rm{e}}}}{{m}_{j}}\right)\left(\displaystyle \frac{{Z}_{j}^{2}{n}_{j}}{{n}_{{\rm{e}}}}\right)\\ & & \times \left[\mathrm{Erf}(\sqrt{x})-2\left(1+\displaystyle \frac{{m}_{j}}{{m}_{i}}\right)\sqrt{\displaystyle \frac{x}{\pi }}{e}^{-x}\right],\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}&&{D}_{{ij}}^{\mathrm{hot}}(E)\displaystyle \frac{{m}_{{\rm{e}}}{c}^{2}}{{\tau }_{\mathrm{Coul}}\beta }\left(\displaystyle \frac{{m}_{{\rm{e}}}}{{m}_{j}}\right)\left(\displaystyle \frac{{Z}_{j}^{2}{n}_{j}}{{n}_{{\rm{e}}}}\right){kT}\left[\mathrm{Erf}(\sqrt{x})-2\sqrt{\displaystyle \frac{x}{\pi }}{e}^{-x}\right],\end{eqnarray} \tag{ 15 }$$

where $x=({E}_{i}/{{kT}}_{j})({m}_{j}/{m}_{i})$ and Erf stands for the error function.

As mentioned above, our numerical code solves the first form of the FP equation, given in Equation (1), which has the energy gain (loss) coefficient $A=\tilde{A}-{{dD}}_{\mathrm{EE}}/{dE}.$ This means that the Coulomb energy loss rate should also be transformed to ${\dot{E}}_{{ij}}^{\mathrm{eff}}={\dot{E}}_{{ij}}^{\mathrm{hot}}+\frac{{{dD}}_{{ij}}^{\mathrm{hot}}}{{dE}}$ or

$$\begin{eqnarray}{\dot{E}}_{{ij}}^{\mathrm{eff}} & = & \left(\displaystyle \frac{{m}_{{\rm{e}}}{c}^{2}}{{\tau }_{\mathrm{Coul}}\beta }\right)\left(\displaystyle \frac{{m}_{{\rm{e}}}}{{m}_{j}}\right)\left(\displaystyle \frac{{Z}_{j}^{2}{n}_{j}}{{n}_{{\rm{e}}}}\right)\\ & & \times \left(\mathrm{Erf}(\sqrt{x})-2\sqrt{\displaystyle \frac{x}{\pi }}{e}^{-x}\right)\left(1-\displaystyle \frac{{{kT}}_{j}}{2E}\right).\end{eqnarray} \tag{ 16 }$$

Note that at ${E}_{i}/{{kT}}_{j}=1/2$ the effective Coulomb loss timescale diverges and for ${E}_{i}/{{kT}}_{j}\lt 1/2$ the rate ${\dot{E}}_{\mathrm{hot}}^{\mathrm{eff}}\lt 0$, which means that the test particle actually gains energy from the target particles. Figure 1 shows the effective Coulomb loss timescales, defined as ${E}_{i}/{\dot{E}}_{{ij}}^{\mathrm{eff}},$ as a function of ${E}_{i}/{{kT}}_{j}$ for the four different pairs of test and target particles (e–e, p–p, e–p, p–e).⁵ There are several features in these equations that deserve some discussion.

• 1.  
  The first feature is that these rates and timescales, although obtained using the nonrelativistic Maxwellian distribution approximation, give the correct (cold target) loss timescales for test particles with ${E}_{i}\gg {{kT}}_{j}$.
• 2.  
  These forms of the Coulomb rates also satisfy the time-independent FP equations for both species of particles, when their distributions are Maxwellian of the same temperature (see Appendix A), and as shown in Appendix B they conserve energy.
• 3.  
  As also can be deduced from the above equations, the shapes of the timescales shown in Figure 1 are invariant and independent of kT but they scale with the temperature of the target particles as ${({{kT}}_{j})}^{3/2}$ (see Appendix C). Thus, although for the purpose of comparison with the acceleration timescale we have used ${kT}=2\;\mathrm{keV},$ the loss times as scaled on the vertical axis are valid for all temperatures.
• 4.  
  Finally the relative values of the timescales play an important role. When electrons are energized the acceleration timescale should be comparable to the shortest e–e collision loss time. In this case the electrons gain energy but share very little of it with protons since the e–p timescale is ${m}_{{\rm{p}}}/{m}_{{\rm{e}}}=1836$ times longer. But as electrons are heated, the e–e timescale increases as ${({{kT}}_{{\rm{e}}})}^{3/2}$ and when ${{kT}}_{{\rm{e}}}\to 150\times {{kT}}_{{\rm{p}}}$ some of the energy goes into protons. This is why PE08, limiting their calculations to lower temperatures, did not need to deal with a coupled kinetic equation. On the other hand, when protons are energized, which is the case we are considering here, the initial phase is similar (p–p time is shorter than p–e) and only protons gain energy. However, as can be seen from Figure 1, once the proton temperature increases by a factor less than ten the p–e interactions become important and some of the energy goes to electrons, which are thermalized quickly because of the much shorter e–e timescale. For this reason, as stated above, we need to solve the more complex coupled FP kinetic equations of protons and electrons even though the acceleration mechanism is energizing only the protons. The new algorithm used for this case is described next.

2.2. Algorithm

In this section we show the results from using the above relations and algorithm to investigate whether SA can generate a significant and discernible nonthermal proton tail without excessively heating the plasma. We assume that protons and electrons are both initially in equilibrium with a Maxwellian distribution of ${{kT}}_{0}=2\;\mathrm{keV},$ and interact both internally and mutually via Coulomb collisions while protons undergo SA by turbulence. As explained in Section 2, the acceleration model is specified by the choice of parameters q, E_c, and ${\tau }_{0}.$ We give E_c in units of ${{kT}}_{0},$ and ${\tau }_{0}$ is expressed in terms of ${\tau }_{\mathrm{Coul}}.$ In addition to the other relevant times, ${\tau }_{\mathrm{pp}}\propto {{kT}}_{{\rm{p}}}^{1.5}$ and ${\tau }_{\mathrm{pe}}\propto {{kT}}_{{\rm{e}}}^{1.5}$ in the plateau region; for the assumed initial temperature ${{kT}}_{0}=2\;\mathrm{keV}$ the plateau is at $\sim {\tau }_{\mathrm{Coul}}$ (see Figure 1).

In what follows we show the time evolution of proton spectra for several values of the acceleration parameters $q,{E}_{c}/{kT}$, and ${\tau }_{0},$ at evenly spaced time steps, beginning with an initial Maxwellian distribution with ${{kT}}_{0}=2\;\mathrm{keV},$ which shifts monotonically to higher temperature and develops an increasing nonthermal component resembling the so-called κ distributions. Electron spectra, which tend to be purely Maxwellian because of the short e–e collision times, are shown only in two interesting cases. These spectra are fit to a thermal component of temperature ${{kT}}_{\mathrm{th}}$ plus an excess that we call the nonthermal component. We give the temperature at the final step, and using the method described in Appendix B of PE08, we also give the fraction of the number of protons in the thermal component ${N}_{\mathrm{th}},$ the ratio of the nonthermal to total energy ${R}_{\mathrm{nonth}},$ and an approximate power-law index δ of the nonthermal component (obtained by a fit to the excess spectrum at energies one order of magnitude above the energy at which it peaks).

As in PE08, we also find that acceleration models with $q\leqslant 0,$ which have stronger acceleration at lower energies where thermalization is fast, lead primarily to heating and a small nonthermal component. Two examples of these are shown in Figure 2 for models with ${E}_{c}=25.5{kT},{\tau }_{0}=2.3{\tau }_{\mathrm{Coul}}$, and $q=-1$ (left) and q = 0 (right). In both models (and others with similar values of q) significant deviations from an isothermal distribution appear after the plasma temperature is increased by several decades, which happens after $t\simeq 10\times {\tau }_{\mathrm{Coul}}\sim {10}^{8}({10}^{-3}{\mathrm{cm}}^{-3}/n)\;\mathrm{yr}.$

Zoom In Zoom Out Reset image size

In what follows we focus on q = 1 and show how the time evolution varies with ${\tau }_{0}$ in Figure 3 and with E_c in in Figure 4. For the first set of models with ${\tau }_{{ac}}$ considerably longer than thermalizing timescales of protons (see Equations (26) and (28)), most of the energy input from the turbulence goes into heating rather than development of a distinct nonthermal tail, producing broad distributions similar to the q = −1 and 0 models. A considerable nonthermal component emerges at late times after the temperature of the thermal component (or the average energy of the distribution) becomes large enough so that the thermalization timescales become considerably longer than ${\tau }_{{ac}}$.⁶ The second set of models with a lower value of E_c exhibit a faster increase of temperature and a more prominent nonthermal component. In each case, the spectrum starts to show a significant nonthermal component only when the proton temperature approaches E_c. This behavior is due to the fact that ${\tau }_{{ac}}\to {\tau }_{0}$ for $E\gg {E}_{c}.$ On the other hand, when ${\tau }_{0}$ (or ${\tau }_{{ac}},$ to be more precise) is comparable to or shorter than ${\tau }_{\mathrm{pe}},$ which is necessary for development of a distinct nonthermal tail, one can avoid excessive heating.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Two examples that have large E_c but ${\tau }_{0}$ comparable to or shorter than ${\tau }_{\mathrm{pe}}$ are shown in Figure 5. The left panel shows the later time evolution of the proton spectrum for the model q = 1, ${\tau }_{0}=5.8{\tau }_{\mathrm{Coul}},$ and ${E}_{c}=125{{kT}}_{0}.$ Although the spectrum at t = $57{\tau }_{\mathrm{Coul}}$ exhibits a promising nonthermal tail and a temperature that is only few times larger than T₀, its shape changes rapidly as a considerable fraction of the protons are accelerated beyond E_c. By the time t = $16{\tau }_{0},$ the spectrum is almost completely nonthermal with a very small thermal component. On the right panel we show the evolution of the proton spectrum for the model q = 1, ${\tau }_{0}=1.2{\tau }_{\mathrm{Coul}},$ ${E}_{c}=255{{kT}}_{0}.$ In this case the spectrum remains nearly thermal until t = $8{\tau }_{0},$ then begins to exhibit an appreciable nonthermal tail at t = $16{\tau }_{0},$ and very quickly becomes strongly nonthermal afterwards. The common feature of the time evolution for these two models is that the proton spectrum stays largely thermal before ${{kT}}_{{\rm{p}}}\sim {E}_{c}$ and then gets very rapidly accelerated once a considerable fraction of protons obtains kinetic energy larger than E_c, leading to a runaway spectrum.

Zoom In Zoom Out Reset image size

These results imply that, in order to get a nonthermal tail by an acceleration model with q = 1, small ${\tau }_{0}$, and large E_c, a great deal of fine-tuning in the duration of acceleration would be necessary. Also, we may see from these results that E_c cannot be too high compared to the desired final temperature of the thermal component because the spectrum will start to develop a considerable nonthermal tail only as it comes close to E_c. The above conclusions will also hold for acceleration models with $q\gt 1,$ whose energy dependence at low energy (i.e., $E\lt {E}_{c}$) is even steeper than in q = 1 models.

As expected from the discussion in Section 2.2, the electron spectrum, in every case considered above, remains nearly thermal with its temperature being usually smaller than that of the proton thermal component during the early times. But once the proton energy exceeds 20 keV p–e interactions become important and electrons are quickly heated to proton temperatures. This is true regardless of how nonthermal the proton spectrum is. In Figure 6 we show the evolution of the electron spectra for the same models as shown in Figure 5. As expected, in some cases and at later times (when ${\tau }_{\mathrm{pe}}\lt {\tau }_{\mathrm{pp}}$) the  temperature of the electrons could be much larger than that of the proton thermal component, especially when the proton spectrum has a prominent nonthermal component. And in some cases the electrons are heated to relativistic temperatures. In such a case the electron bremsstrahlung will exceed the proton (inverse) bremsstrahlung. This is another difficulty with production of hard X-rays by a proton rather than an electron nonthermal tail.

Zoom In Zoom Out Reset image size

In summary, the above results strongly suggest that producing a significant proton nonthermal tail requires a significant fine-tuning of the acceleration parameters, the duration of the acceleration, and/or the initial distributions. As a result, for example, maintaining such a spectrum that could produce hard X-rays in the ICM over a Hubble time would be difficult.

In this paper, we have explored the possibility of development of nonthermal tails in the spectrum of protons starting from a background (nonrelativistic) thermal plasma. We use a generic phenomenological energizing mechanism with a simple three-parameter form for energy diffusion and acceleration rates. As briefly described after Equation (4) in Section 2, these forms can simulate a broad range of acceleration processes. The energizing is opposed by Coulomb energy losses. Such a process may be in operation in the initial phase of most acceleration mechanisms and may be important in collisional plasmas such as solar and stellar flares and in the ICM of galaxy clusters, and could produce nonthermal bremsstrahlung radiation in the hard X-ray regime in collisions with the background electrons.

We derived the Coulomb energy-loss and diffusion coefficients in a general setting of a "hot plasma" where the distribution of target particles and the masses of the test and target particles are arbitrary. We present analytic formulae for a Maxwellian distribution that asymptotically approach the cold-target relations at high energies. We have generalized the algorithm developed by PE08 to solve the coupled e–p kinetic equations. We show that it satisfies the steady-state equilibrium distributions (Appendix A) and the energy conservation (Appendix B). In Appendix C we test our algorithm by demonstrating that in the absence of acceleration the particle distributions relax to Maxwellian within the expected timescales starting from several non-equilibrium initial conditions. In Appendix D we discuss the inelastic energy loss rates and present a new analytic formula for the proton energy loss rate due to pion production.

Using the algorithm, we obtain the time evolution of the proton and electron distributions. Because of the shorter collision time of electrons the production of nonthermal electron tails requires a higher rate of acceleration and shorter acceleration time than those used here for the protons. Thus, in our simulations we include only acceleration of protons, which may be the case in the presence of low-frequency Alfvénic turbulence. The electrons gain energy from protons but thermalize quickly. The resulting proton spectra are decomposed into thermal and nonthermal components according to the fitting methods of PE08. Our results can be summarized as follows:

• 1.  
  For q = −1 and 0, even when the acceleration time ${\tau }_{0}$ is sufficiently small that the energizing rate is comparable to or faster than the thermalization rate, an initially Maxwellian proton spectrum evolves into a hotter one and a broader than thermal distribution with little sign of a distinct nonthermal tail. This is because for $q\leqslant 0$ the energizing rate is smaller at higher energies and not efficient in acceleration into a high-energy tail. This suggests that q must be positive for production of a distinct nonthermal tail and a mild temperature increase. This means we need turbulence with a spectrum considerably steeper than Kolmogorov.
• 2.  
  We show that for q = 1 the temperature increases more slowly as long as most of the particles have kinetic energy below E_c. In addition, a significant nonthermal tail may appear if ${\tau }_{0}$ is comparable to or smaller than the proton thermalization time. For such small values of ${\tau }_{0},$ however, once a considerable fraction of the protons have reached E_c, they rapidly evolve into a completely nonthermal, runaway distribution. This will be more pronounced for even higher values of q, where the acceleration rate increases more rapidly with energy.
• 3.  
  The electron spectrum remains nearly thermal in every case considered in this paper. For cases where the proton distribution is nearly Maxwellian the electron temperature remains somewhat below the proton temperature. But when the proton distribution contains a significant or a larger nonthermal component (and ${\tau }_{\mathrm{pe}}\lt {\tau }_{\mathrm{pp}}$ for most protons) the electron temperature can reach the average energy of the protons and be much higher than that of the proton thermal component.

In summary, the above results show that a great deal of fine-tuning and a steeper than Kolmogorov spectrum of turbulence are necessary to create a significant nonthermal proton tail without excessive heating of protons and/or electrons even when only protons are energized. A corollary of this is that, contrary to some claims (Boldt & Serlemistos 1969; Wolfe & Melia 2006), nonthermal inverse bremsstrahlung cannot be a significant contributor to hard X-ray emissions (see also Emslie & Brown 1985).

We thank Qingrong Chen for his extensive help in developing the numerical code. This work was supported by Stanford University VPUE Major Grant No. 4842.

When two species of particles denoted by i and j are in thermal equilibrium, the time-independent FP equation ($\partial N/\partial t=0$) for each of the two of the species, say i, is written as

$$\begin{eqnarray}0 & = & \displaystyle \frac{{\partial }^{2}}{\partial {E}^{2}}[({D}_{{ii}}^{\mathrm{hot}}(E)+{D}_{{ij}}^{\mathrm{hot}}(E)){N}_{i}(E)]\\ & & +\displaystyle \frac{\partial }{\partial E}[({\dot{E}}_{{ii}}^{\mathrm{hot}}(E)+{\dot{E}}_{{ij}}^{\mathrm{hot}}(E)){N}_{i}(E)],\end{eqnarray} \tag{ 17 }$$

where N_i(E) is the (Maxwellian) distribution of species i and all the Coulomb coefficients are evaluated at the same temperature as that of ${N}_{i}(E).$ Since N_i(E) must be time-independent even in the absence of Coulomb interactions with the other species of particles, i.e.,

$$\begin{eqnarray}&&0=\displaystyle \frac{{\partial }^{2}}{\partial {E}^{2}}[{D}_{{ii}}^{\mathrm{hot}}(E){N}_{i}(E)]+\displaystyle \frac{\partial }{\partial E}[{\dot{E}}_{{ii}}^{\mathrm{hot}}(E){N}_{i}(E)],\end{eqnarray} \tag{ 18 }$$

it follows that, in particular,

$$\begin{eqnarray}&&0=\displaystyle \frac{{\partial }^{2}}{\partial {E}^{2}}[{D}_{{ij}}^{\mathrm{hot}}(E){N}_{i}(E)]+\displaystyle \frac{\partial }{\partial E}[{\dot{E}}_{{ij}}^{\mathrm{hot}}(E){N}_{i}(E)].\end{eqnarray} \tag{ 19 }$$

Given N_i(E) and ${\dot{E}}_{{ij}}^{\mathrm{hot}}$ we can integrate this equation to uniquely determine ${D}_{{ij}}^{\mathrm{hot}}(E)={\displaystyle \int }_{E}^{\infty }{\dot{E^{\prime} }}_{{ij}}^{\mathrm{hot}}(E^{\prime} ){N}_{i}(E^{\prime} ){dE}^{\prime} /{N}_{i}(E).$ The diffusion coefficient uniquely determined in this way agrees exactly with those presented in Section 2, which justifies the assumptions used in their derivations.

Here, we show that, for arbitrary particle distributions and FP coefficients given by Equations (9) and (12), our algorithm described in Section 2 satisfies the total energy conservation. Let E_tot be the total energy of the two species of particles. Then, from their FP equations with a no-flux boundary condition (see Equation (3) of Park & Petrosian 1995), we find

$$\begin{eqnarray}&&\displaystyle \frac{\partial {E}_{\mathrm{tot}}}{\partial t}=-(A+B+C+D),\end{eqnarray} \tag{ 20 }$$

where

$$\begin{eqnarray}&&A={\displaystyle \int }_{0}^{\infty }{\displaystyle \int }_{0}^{\infty }{N}_{i}(E){G}_{{ii}}(E,E^{\prime} ){N}_{i}(E^{\prime} ){dEdE}^{\prime} ,\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&B={\displaystyle \int }_{0}^{\infty }{\displaystyle \int }_{0}^{\infty }{N}_{i}(E){G}_{{ij}}(E,E^{\prime} ){N}_{j}(E^{\prime} ){dEdE}^{\prime} ,\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&C={\displaystyle \int }_{0}^{\infty }{\displaystyle \int }_{0}^{\infty }{N}_{j}(E){G}_{{ji}}(E,E^{\prime} ){N}_{i}(E^{\prime} ){dEdE}^{\prime} ,\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&D={\displaystyle \int }_{0}^{\infty }{\displaystyle \int }_{0}^{\infty }{N}_{j}(E){G}_{{jj}}(E,E^{\prime} ){N}_{j}(E^{\prime} ){dEdE}^{\prime} .\end{eqnarray} \tag{ 24 }$$

A and B separately vanish because of antisymmetry of G_ii and G_jj, and because ${G}_{{ij}}(E,E^{\prime} )=-{G}_{{ji}}(E^{\prime} ,E),$ it follows that $B+C=0.$ This proves the claim.

As a test of our algorithm, we consider the thermalization of the proton and electron energy distributions that are initially not in thermal equilibrium. The total particle numbers of protons and electrons are assumed to be the same and normalized to 1. The proton and electron distributions are expected to converge to Maxwellian distributions of the same temperature over certain thermalization timescales, conserving the total energy at each time. Based on Equation (14), we may define four thermalization timescales as

$$\begin{eqnarray}&&{\tau }_{\mathrm{ee}}({{kT}}_{{\rm{e}}})={\tau }_{\mathrm{Coul}}{\left(\displaystyle \frac{{{kT}}_{{\rm{e}}}}{{m}_{{\rm{e}}}{c}^{2}}\right)}^{3/2}\sim |\displaystyle \frac{{{kT}}_{{\rm{e}}}}{{\dot{E}}_{\mathrm{ee}}^{\mathrm{hot}}({{kT}}_{{\rm{e}}})}|,\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{\tau }_{\mathrm{pp}}({{kT}}_{{\rm{p}}})={\tau }_{\mathrm{Coul}}\sqrt{\displaystyle \frac{{m}_{{\rm{p}}}}{{m}_{{\rm{e}}}}}{\left(\displaystyle \frac{{{kT}}_{{\rm{p}}}}{{m}_{{\rm{e}}}{c}^{2}}\right)}^{3/2}\sim |\displaystyle \frac{{{kT}}_{{\rm{p}}}}{{\dot{E}}_{\mathrm{pp}}^{\mathrm{hot}}({{kT}}_{{\rm{p}}})}|,\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{\tau }_{\mathrm{ep}}({{kT}}_{{\rm{p}}},{{kT}}_{{\rm{e}}})={\tau }_{\mathrm{Coul}}\left(\displaystyle \frac{{m}_{{\rm{p}}}}{{m}_{{\rm{e}}}}\right){\left(\displaystyle \frac{{{kT}}_{{\rm{e}}}}{{m}_{{\rm{e}}}{c}^{2}}\right)}^{3/2}\sim |\displaystyle \frac{{{kT}}_{{\rm{e}}}}{{\dot{E}}_{\mathrm{ep}}^{\mathrm{hot}}({{kT}}_{{\rm{e}}})}|,\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{\tau }_{\mathrm{pe}}({{kT}}_{{\rm{p}}},{{kT}}_{{\rm{e}}})={\tau }_{\mathrm{Coul}}\left(\displaystyle \frac{{m}_{{\rm{p}}}}{{m}_{{\rm{e}}}}\right){\left(\displaystyle \frac{{{kT}}_{{\rm{p}}}}{{m}_{{\rm{e}}}{c}^{2}}\right)}^{3/2}\sim |\displaystyle \frac{{{kT}}_{{\rm{p}}}}{{\dot{E}}_{\mathrm{pe}}^{\mathrm{hot}}({{kT}}_{{\rm{p}}})}|,\end{eqnarray} \tag{ 28 }$$

where kT_p and kT_e are defined as 2/3 of the total energy of the initial distributions.⁷ Interestingly, as can be seen from Figure 1, for ${T}_{{\rm{e}}}\sim {T}_{{\rm{p}}}$ the interspecies interaction times ${\tau }_{\mathrm{ep}}$ and ${\tau }_{\mathrm{pe}}$ are almost equal and are $\sqrt{{m}_{{\rm{p}}}/{m}_{{\rm{e}}}}$ and ${m}_{{\rm{p}}}/{m}_{{\rm{e}}}$ times longer than ${\tau }_{\mathrm{pp}}$ and ${\tau }_{\mathrm{ee}},$ respectively. Thus thermalization is governed by the former timescales.

First, we consider the thermalization when the protons and electrons have different initial temperatures. Figure 7 shows the time evolution of the proton and electron distributions with initial temperatures kT_p = 25 keV and kT_e = 5 keV. As expected, both species remain Maxwellian throughout and their temperatures approach the average of the two initial temperatures over several ${\tau }_{\mathrm{pe}}$'s and ${\tau }_{\mathrm{ep}}$'s, respectively calculated at the initial temperature. We get similar results with initial values reversed: kT_p = 5 keV and kT_e = 25 keV. Next, we consider the thermalization of the protons and electrons when they initially have a narrow log-normal distributions. The normalized log-normal distribution function is given by

$$\begin{eqnarray}&&f(E;{E}_{0},\sigma )=\displaystyle \frac{1}{E\sigma \sqrt{2\pi }}\mathrm{exp}\left(-\displaystyle \frac{\mathrm{ln}E/{E}_{0}}{2{\sigma }^{2}}\right),\end{eqnarray} \tag{ 29 }$$

whose mean is ${E}_{0}{e}^{{\sigma }^{2}/2}.$ Note that, for $\sigma \ll 1,$ ${E}_{0}{e}^{{\sigma }^{2}/2}\approx {E}_{0},$ and the log-normal distribution becomes essentially a Dirac delta function centered at E₀. In this case, the proton and electron distributions should approach the Maxwellian shape over several ${\tau }_{\mathrm{pp}}$'s and ${\tau }_{\mathrm{ee}}$'s respectively. At the same time, their temperatures should also converge to $3{E}_{0}/2$ over several ${\tau }_{\mathrm{pe}}$'s and ${\tau }_{\mathrm{ep}}$'s respectively. In Figure 8 we show the time evolution of the proton and electron distributions, when the protons and electrons are initially in log-normal distributions with $(2/3){E}_{0,{\rm{e}}}$= 5 keV, ${\sigma }_{{\rm{p}}}$ = 0.02 and $(2/3){E}_{0,{\rm{p}}}$ = 25 keV, ${\sigma }_{{\rm{e}}}$ = 0.02, respectively. Not shown here, within tens of ${\tau }_{\mathrm{ee}}\ll {\tau }_{\mathrm{ep}}$ and ${\tau }_{\mathrm{pp}}\ll {\tau }_{\mathrm{ep}},$ respectively, electrons and protons acquire Maxwellian distributions (as can be seen in this figure, this is the case at $t\sim {\tau }_{\mathrm{ep}}$ and $t\sim {\tau }_{\mathrm{pe}},$ respectively). However, it takes several tens of these timescale before they both reach the final equilibrium of ${{kT}}_{{\rm{e}}}={{kT}}_{{\rm{p}}}=15\;\mathrm{keV}.$ Similar results are obtained with initially more energetic protons. Similar results are obtained starting with 25 keV protons and 5 keV electrons.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The inelastic energy loss rates for electrons are due to well known processes of synchrotron, IC (with $\dot{E}\propto u\times {\beta }_{{\rm{e}}}^{2}{\gamma }_{{\rm{e}}}^{2},$ where u is the total energy density of magnetic field and soft photons), and electron–proton bremsstrahlung with

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{ep}}^{\mathrm{br}}=\left(\displaystyle \frac{4\alpha }{3\pi \mathrm{ln}{\rm{\Lambda }}}\right)\displaystyle \frac{{m}_{{\rm{e}}}{c}^{2}}{{\tau }_{\mathrm{Coul}}}{\beta }_{{\rm{e}}}{\gamma }_{{\rm{e}}}\chi ({\gamma }_{{\rm{e}}}),\end{eqnarray} \tag{ 30 }$$

where $\alpha =1/137$ is the fine structure constant and $\chi (\gamma )\simeq a+b\mathrm{ln}\gamma$ is a slowly varying function. At relativistic energies this rate should be increased by nearly an equal amount due to e–e bremsstrahlung. Thus, bremsstrahlung loss becomes dominant compared to elastic e–e scattering at Lorentz factors $\geqslant 23\mathrm{ln}{\rm{\Lambda }}.$ As shown in Figure 9, adding these rate modifies the electron energy timescales shown in Figure 1 at high energies.

Zoom In Zoom Out Reset image size

For protons the main inelastic loss rate is due to hadronic interactions. For $\lt 30\;\mathrm{MeV}$ protons the main loss is due to the proton–ion interactions producing ionic de-excitation lines in the 1–7 MeV range. The loss rate of $\gt 30$ to ∼300 MeV protons is dominated by neutron production. Pion production starts at threshold energy ${E}_{\mathrm{th}}\sim 300\;\mathrm{MeV}$ and dominates all losses (including Coulomb) above a few GeV. Proton–electron interactions can produce X- and gamma-rays via the (inverse) bremsstrahlung process. Since energy loss rate is an invariant quantity it follows that ${\dot{E}}_{\mathrm{pe}}^{\mathrm{br}}$ will be equal to that of e–p bremsstrahlung in the rest frame of the proton, so that this rate is also given by Equation (30) with electron velocity and Lorentz factor replaced with those of protons. The inverse bremsstrahlung loss becomes negligible for proton velocities less than the electron thermal velocity, similar to the Coulomb p–e loss (see Figure 1).

One can get a simple approximate expression for the energy loss rate due to p–p hadronic interactions at relativistic energies but there are no simple formulae for semi-relativistic energies and at energies near the threshold. We combine previously published results on these interaction and derive simple formulae in three energy regimes with an accuracy of a few per cent at the boundaries of these regimes.

D.1.  ${E}_{h}=100\lt E\lt {10}^{8}$

For this very high energy regime, we make use of the work of Kelner et al. (2006) on the energy spectra of the secondary particles produced by p–p interactions. If we denote the number dN_i of secondaries in the kinetic energy interval $({E}_{i},{E}_{i}+{{dE}}_{i})$ produced by a proton of energy E_p (in units of m_pc²) by ${{dN}}_{i}={F}_{i}({E}_{i},{E}_{{\rm{p}}}){{dE}}_{i},$ then the energy loss rate can be written as

$$\begin{eqnarray}&&{\dot{E}}_{{\rm{p}}-{\rm{p}}}=\beta {cn}{\sigma }_{\mathrm{inel}}({E}_{{\rm{p}}})\cdot \left(\displaystyle \sum _{i}{\displaystyle \int }_{0}^{{E}_{{\rm{p}}}}{F}_{i}({E}_{i},{E}_{{\rm{p}}}){E}_{i}{{dE}}_{i}\right),\end{eqnarray} \tag{ 31 }$$

where the index i is summed over all secondaries and ${\sigma }_{\mathrm{inel}}({E}_{{\rm{p}}})$ is the inelastic p–p interaction cross section that encodes information as to how many inelastic p–p interactions occur (not the number of secondaries produced, as some cross sections defined elsewhere do).

There are two choices when counting the secondaries. We can count the relatively stable secondaries, which are gamma-rays, electrons and electron and muon neutrinos (including their antiparticles). Alternatively, we can count the intermediate secondaries that decay into the above-mentioned more stable secondaries, which are pions and eta mesons. The energy loss rates calculated in both ways agree with each other within a few per cent, and can be approximated as

$$\begin{eqnarray}&&{\dot{E}}_{{\rm{p}}-{\rm{p}}}=C\times {({E}_{{\rm{p}}}/{E}_{h})}^{1.09},\\ &&\mathrm{where}\;\;C=1.26\times {10}^{-24}{{cnm}}_{{\rm{p}}}{c}^{2}{\beta }_{{\rm{p}}}.\end{eqnarray} \tag{ 32 }$$

D.2.  ${E}_{l}=2.2\lt E\lt {E}_{h}=100$

For this intermediate energy regime, we adopt Stephens and Badhwar's model (Stephens & Badhwar 1981; see also Dermer 1986), where the pions are the main particles produced and we get equal contributions from each pion to the sum in Equation (31) with

$$\begin{eqnarray}{F}_{\pi }({E}_{\pi },{E}_{{\rm{p}}}) & = & \displaystyle \frac{2\pi \sqrt{{E}_{\pi }^{2}-{m}_{\pi }^{2}}{c}^{4}}{\langle \xi {\sigma }_{\pi }(E)\rangle }\\ & & \times {\displaystyle \int }_{\mathrm{cos}{\theta }_{\mathrm{max}}}^{1}({E}_{\pi }{d}^{3}\sigma /d{{p}_{\pi }}^{3})d\mathrm{cos}\theta ,\end{eqnarray} \tag{ 33 }$$

where $\mathrm{cos}{\theta }_{\mathrm{max}}=({\gamma }_{c}{E}_{\pi }-{E}_{\mathrm{max}}^{*})/({\beta }_{c}{\gamma }_{c}{p}_{\pi }),$ ${E}_{\mathrm{max}}^{*}=[{E}_{\mathrm{tot}}^{2}-(4{m}_{{\rm{p}}}^{2}+{m}_{\pi }^{2}){c}^{4}]$ $/(2{E}_{\mathrm{tot}}),$ and

$$\begin{eqnarray}&&{E}_{\pi }\displaystyle \frac{{d}^{3}\sigma }{d{{p}_{\pi }}^{3}}={Af}(E){(1-\tilde{x})}^{q}\mathrm{exp}[-{{Bp}}_{\perp }],\end{eqnarray} \tag{ 34 }$$

where $q={C}_{1}+{C}_{2}{p}_{\perp }+{c}_{3}{p}_{\perp }^{2},$ with experimentally derived values of A = 140, $B=5.43,$ ${C}_{1}=6.1,$ ${C}_{2}=3.3,$ and ${C}_{3}=0.6.$ In these equations ${\tilde{x}}^{2}={({p}_{\parallel }/{p}_{\mathrm{max}}^{*})}^{2}+4({p}_{\perp }^{2}{c}^{2}+{m}_{\pi }^{2}{c}^{4})/{E}_{\mathrm{tot}}^{2};$ all quantities are measured in the center-of-mass system (c.m.s.). Using these equations we find an energy loss rate ${\dot{E}}_{{\rm{p}}-{\rm{p}}}=C\times [({E}_{{\rm{p}}}-{E}_{l})/{E}_{h}+0.012]$

D.3.  ${E}_{\mathrm{th}}=0.3\lt E\lt {E}_{l}=2.2$

For the lowest energy range down to the threshold kinetic energy we use Stecker's model (Stecker 1970; also see Dermer 1986), where we can approximate the ${\pi }^{0}$ production process by the so-called "isobar-plus-fireball model."⁹ In this model, we will assume that all pions are produced in two ways: via intermediate production and decay of the ${\rm{\Delta }}(1.238)$ non-strange isobar, and from a thermal pion gas created from the remaining available energy in the c.m.s. of the collision where the pions are given an energy distribution very similar to a Maxwell–Boltzmann type distribution. In particular, for $E\lt 2.2,$ we can assume that all pion production occurs through the isobar production mode because it is the most dominant mode in this energy regime (the fireball mode becomes dominant when the total energy of the incident proton is greater than 5 GeV). We further assume that the isobars of mass ${m}_{{\rm{\Delta }}}$ carry momentum either directly forward or directly backward in the c.m.s., and that the low-energy isobar production process produces pions through the two-stage decay

$$\begin{eqnarray*}&&p+p\to p+{\rm{\Delta }}(1.238)\\ &&\downarrow \\ &&p+{\pi }^{0}.\end{eqnarray*}$$

The normalized production spectrum of pions, to be used in Equation (31) in this model, is given by

$$\begin{eqnarray}&&{F}_{\pi }({E}_{\pi },{E}_{{\rm{p}}})={w}_{r}({E}_{{\rm{p}}}){\displaystyle \int }_{{m}_{{\rm{p}}}+{m}_{{\pi }^{0}}}^{{s}^{1/2}-{m}_{{\rm{p}}}}{{dm}}_{{\rm{\Delta }}}B({m}_{{\rm{\Delta }}})f({E}_{\pi },E,{m}_{{\rm{\Delta }}}),\end{eqnarray} \tag{ 35 }$$

where

$$\begin{eqnarray}f({E}_{\pi },E,{m}_{{\rm{\Delta }}}) & = & \displaystyle \frac{1}{2{m}_{{\rm{\Delta }}}}\left\{{(2{\beta }_{{\rm{\Delta }}}^{+}{\gamma }_{{\rm{\Delta }}}^{+}{\beta }_{\pi }^{\prime }{\beta }_{\pi }^{\prime })}^{-1}H\right.\\ & & \times [{\gamma }_{\pi };{\gamma }_{{\rm{\Delta }}}^{+}{\gamma }_{\pi }^{\prime }(1-{\beta }_{{\rm{\Delta }}}^{+}{\beta }_{\pi }^{\prime }),{\gamma }_{{\rm{\Delta }}}^{+}{\gamma }_{\pi }^{\prime }(1+{\beta }_{{\rm{\Delta }}}^{+}{\beta }_{\pi }^{\prime })]\\ & & +{(2{\beta }_{{\rm{\Delta }}}^{-}{\gamma }_{{\rm{\Delta }}}^{-}{\beta }_{\pi }^{\prime }{\beta }_{\pi }^{\prime })}^{-1}H\\ & & \times \left[{\gamma }_{\pi };{\gamma }_{{\rm{\Delta }}}^{-}{\gamma }_{\pi }^{\prime }(1-{\beta }_{{\rm{\Delta }}}^{-}{\beta }_{\pi }^{\prime }),{\gamma }_{{\rm{\Delta }}}^{-}{\gamma }_{\pi }^{\prime }\right.\\ & & \times \left.\left.(1+{\beta }_{{\rm{\Delta }}}^{-}{\beta }_{\pi }^{\prime })\right]\right\},\end{eqnarray} \tag{ 36 }$$

and the normalized isobar mass spectrum is given by the Breit–Wigner distribution $B({m}_{{\rm{\Delta }}})=({\rm{\Gamma }}/\pi )\left((,{1/({m}_{{\rm{\Delta }}}-{m}_{{\rm{\Delta }}}^{0})}^{2}+{{\rm{\Gamma }}}^{2},)\right)$ with the normalization factor ${w}_{r}(E)=\pi [\mathrm{arctan}({s}^{1/2}-{m}_{{\rm{p}}}-{m}_{{\rm{\Delta }}}^{0}/{\rm{\Gamma }})-\mathrm{arctan}$ ${({m}_{{\rm{p}}}+{m}_{{\pi }^{0}}-{m}_{{\rm{\Delta }}}^{0}/{\rm{\Gamma }})]}^{-1}.$ However, since Γ is sufficiently small, we can approximate $B({m}_{{\rm{\Delta }}})$ by the Dirac delta function $\delta ({m}_{{\rm{\Delta }}}-{m}_{{\rm{\Delta }}}^{0})$ and therefore ignore the normalization factor ${w}_{r}(E).$

Here $H[x;a,b]=1$ if $a\leqslant x\leqslant b$ and 0 otherwise, and ${\gamma }_{\pi }$ is the pion Lorentz factor in the lab system. The Lorentz factors of the forward (+) and backward (−) moving isobars are ${\gamma }_{{\rm{\Delta }}}^{\pm }={\gamma }_{c}{\gamma }_{{\rm{\Delta }}}^{*}$ $(1\pm {\beta }_{c}{\beta }_{{\rm{\Delta }}}^{*}),$ where ${\gamma }_{c}={s}^{1/2}/2{m}_{{\rm{p}}}$ is the Lorentz factor of the c.m.s. with respect to the lab system, and ${\gamma }_{{\rm{\Delta }}}^{*}=(s+{m}_{{\rm{\Delta }}}^{2}-{m}_{\pi }^{2})/2\;{s}^{1/2}{m}_{{\rm{\Delta }}}$ is the Lorentz factor of the isobar in the c.m.s. The pion Lorentz factor in the rest frame of the Δ-isobar is ${\gamma }_{\pi }^{\prime }=$ $({m}_{{\rm{\Delta }}}^{2}+{m}_{\pi }^{2}-{m}_{{\rm{p}}}^{2})/2{m}_{{\rm{\Delta }}}{m}_{\pi }.$

Again with the help of these equations and Equation (31) we obtain an energy loss rate of $C\times [5.9\times {10}^{-3}({E}_{{\rm{p}}}-{E}_{\mathrm{th}})$ $+\;2.2\times {10}^{-4}{({E}_{{\rm{p}}}-{E}_{\mathrm{th}})}^{2}].$

Now putting all this together we obtain energy loss rate as a function of proton energy E (in units of m_pc²) as

$$\begin{eqnarray}&&{\dot{E}}_{{\rm{p}}-{\rm{p}}}=\displaystyle \frac{1.26{m}_{{\rm{p}}}{c}^{2}\beta }{{\tau }_{\mathrm{Coul}}\mathrm{ln}{\rm{\Lambda }}}\\ \times \;\left\{\begin{array}{ll}0.011(E-{E}_{\mathrm{th}})+2.2 & \\ \quad \times \;{10}^{-4}{(E-{E}_{\mathrm{th}})}^{2}, & {E}_{\mathrm{th}}\leqslant E\leqslant {E}_{l},\\ (E-{E}_{l})/{E}_{h}+0.022, & {E}_{l}\leqslant E\leqslant {E}_{h},\\ {(E/{E}_{h})}^{1.09}, & {E}_{h}\leqslant E\leqslant {10}^{8},\end{array}\right.\end{eqnarray} \tag{ 37 }$$

where ${E}_{\mathrm{th}}=0.3,$ E_l = 2.2, and E_h = 100. In Figure 9 we show energy loss for both electrons and protons, extending it to higher energies than those shown in Figure 1.