Acceleration of Thermal Protons By Generic Phenomenological Mechanisms
Abstract
We investigate heating and acceleration of protons from a thermal gas with a generic diffusion and acceleration model, and subject to Coulomb scattering and energy loss, as was carried out in Petrosian & East (2008) for electrons. As protons gain energy their loss to electrons becomes important. Thus, we need to solve the coupled protonelectron kinetic equation. We numerically solve the coupled FokkerPlank equations and computes the time evolution of the spectra of both particles. We show that this can lead to a quasithermal component plus a high energy nonthermal tail. We determine the evolution of nonthermal tail and the quasithermal component. The results may be used to explore the possibility of inverse bremsstrahlung radiation as a source of hard Xray emissions from hot sources such as solar flares, accretion disk coronas and the intracluster medium of galaxy clusters. We find that emergence of nonthermal protons is accompanied by excessive heating of the entire plasma, unless the turbulence needed for scattering and acceleration is steeper than Kolmogorov and the acceleration parameters, the duration of the acceleration, and/or the initial distributions are significantly finetuned. These results severely constraint the feasibility of nonthermal inverse bremsstrahlung process producing hard Xray emissions. However the nonthermal tail may be the seed particles for further reacceleration to relativistic energies, say by a shock. In the Appendix we present some tests of the integrity of the algorithm used and present a new formula for the energy loss rate due to inelastic protonproton interactions.
1 Introduction
Most particle acceleration mechanisms invoked for astrophysical sources must start with acceleration of low energy particles of the background magnetized plasma with a thermal or Maxwellian distribution confined in a finite volume with density , temperature and magnetic field . This initial phase of acceleration can produce a distribution consisting of a quasithermal 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 socalled thermal leakage injection model in diffusive shock acceleration (see e.g. Haysung et al. 2014). The main goal of this and an earlier paper [Petrosian & East 2008 (PE08)] is to explore the possibility of producing a prominent nonthermal tail by some generic phenomenological acceleration process. PE08 treated the generation of 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 similar energizing mechanism using the coupled FokkerPlanck (FP) kinetic equations.
The nonthermal tails in kappalike electron distributions could also be responsible for nonthermal emission in the hard Xray regime for hot plasmas () 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 signatures of nonthermal activity first observed as radio radiation that is due to synchrotron emission by a population of relativistic (Lorentz factor ) electrons in field. In addition, earlier observations by several hard Xray instruments indicated that, in addition to the well known thermal soft Xray radiation, some clusters show excess radiation above 10’s 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 keV. However, as stressed by Petrosian (2001) and later rigorously proved by EP08, this scenario would cause excessive heating of ICM. Subsequently it was suggested by Wolfe & Melia (2006) that a nonthermal tail in proton distributions extending to 100’s of MeV could also be a source of hard Xrays (10100 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 Xrays in solar flares (see review by Emslie & Brown 1985) because of the possible higher Xray yield compared to electrons. Exploring this possibility is 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 Xrays could also be produced by the inverse Compton (IC) scattering of CMB 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. 1996; 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).^{1}^{1}1All the review articles cited above can be found in the February 2008, Volume 134 issue of the Space Science Reviews, entitled Clusters of Galaxies: Beyond the Thermal View, by Kaastra et al. (2008). However, more recent observations have cast doubt on the reality of the claimed hard Xray 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 nondetection of many clusters by FermiLAT (Ackermann et al. 2011 and 2014; Huber et al. 2011) 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 recent review by Brunetti & Jones 2014).
The crucial feature in production of nonthermal electron or proton tails is the interplay between the particleparticle (and particleexternal 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 K the energy loss is dominated by Coulomb collisions, which tend to equilibrate electrons and protons distribution 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 distributions. Thus, for production of significant but not excessive nonthermal population we need acceleration rates comparable to the Coulomb rates, so that in EP08, dealing with acceleration of electrons only, we used acceleration timescales comparable to electronelectron (ee) Coulomb times. Since the latter are much shorter than that of electronproton (ep) (and protonproton (pp)) 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 temperatures of 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 pp collisional time requires a longer acceleration time. Such a rate for electrons would lead to their heating. Therefore, we assume no, or much lower rate of acceleration of electrons. This can come about for pure Alfvénic turbulence which do 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 protonelectron (pe) and pp collision times can be comparable, so that some of the energy gained by protons can be transfered 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 increasing with energy being more efficient in producing nonthermal tails.
In §2, we derive and present coefficients for Coulomb interactions and stochastic acceleration to be used in a FokkerPlank (FP) kinetic equation, and explain our algorithm for solving the coupled FP equations of protons and electrons. In §3, we apply the algorithm to the stochastic acceleration of thermal background particles, and examine the time evolution of the proton and electron distributions for various acceleration model parameters. In §4, we summarize our results and discuss their implications. Some details of the acceleration and a new formula for proton inelastic energy loss rate due to pion production are presented in the Appendixes.
2 FokkerPlanck Equations and its Coefficients
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 using 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 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 angleaveraged, spatiallyintegrated momentum distribution satisfies the simple equation . For convenience we use the energy distribution defined as to obtain the following commonly used two forms of the kinetic equation:
(1) 
where the energy diffusion coefficient and the direct energy gain (or loss) rates are given as
(2) 
Here is the particle velocity in unit of the speed of light and the Lorentz factor . As shown in Petrosian & Chen (2015) (see also Tsytovich, 1977) provides a more accurate representation of energy gain (or Coulomb energy loss) rate than .
In what follows we will present results from numerical solutions of the first form of the above FP equation including the effects of stochastic acceleration (SA) by turbulence and Coulomb collision. We will assume a background plasma of constant density and field (but not temperature ), and a constant energy level and spectrum of turbulence so that the transport coefficients and are constant in time. Following PE08 we use the simple form
(3) 
or the characteristic acceleration time scale
(4) 
described by the three parameters , , and . The thin black (solid and dotted) curves in Figure 1 show some examples of acceleration times used in the next section.
As described above, , which is a measure of the acceleration time, should be comparable to the relevant collision times and allows 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 ). The break will be important only for low values of and will have little effect for . For SA by turbulence the index is related to the spectral index of turbulence. For relativistic energies and Alfvénic turbulence . In the inertial range one expects a Kolmogorov or Kraichnin spectrum with or 3/2 so that and , respectively. But the turbulence spectrum is expected to be steeper in the damping range at smaller scales (large wave vectors) where 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 and Petrosian & Liu 2004). In what follows we will use to include all possible cases.^{2}^{2}2Although we limit our calculations to the SA by turbulence it should be noted that our results are more general and will be similar to that expected from acceleration by a shock. For example, the energy dependence of acceleration rate by a shock, which depends on the pitch angle scattering, rather than momentum diffusion coefficient (Petrosian 2012), will have similar energy dependence (see Fig. 1 of Petrosian & Chen 2014).
2.1 Coulomb Collisions
In the cold target approximation, in which a test particle of charge , mass , and velocity collides with stationary target particles of charge , mass and number density , the Coulomb energy diffusion rate is negligible but the Coulomb energy loss rate
(5) 
where is the Coulomb logarithm. In what follows we will be interested in electron and proton test particles with . In most astrophysical situations we are dealing with fully ionized ions so the primary target particles are electrons and protons so also. Contributions to Coulomb losses from other ions are negligible except from alpha particles with a relative number density and , which may contribute up to 20 %. Note that this also means that and . If we define , where is the classical electron radius, then
(6) 
The cold target approximation can be used for test particle energies of the background plasma. For interaction with a thermal background plasma with temperature , when the test particle energy approaches , we must use the hot plasma energy exchange rates which yields a finite energy diffusion rate and can be an energy gain when . 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 on interactions of both electrons and protons. This requires generalization of the two functions and of PE08 that describe the loss and diffusion rates, respectively, as follows:
Generalization of Equation (7) in EP08 gives
(7) 
The meaning of is clear; if the test particle is faster than the target particle, it will lose energy to the target particle according to 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) by integrating over the distribution of the target particles:
(8)  
(9) 
In a similar manner, we derive the corresponding Coulomb diffusion coefficient . Generalizing given by Equation (10) of PE08, we define by
(10) 
is given by integrating over the target particle distribution:
(11)  
(12) 
Note that, while is inversely proportional to the mass of the slower particle has no explicit dependence on the masses of any particles and in particular for identical particles is independent of target particle mas or velocity.
When the number of particles in the nonthermal component of is much less than that of the thermal component we can approximate the distribution of the target particles by the nonrelativistic Maxwellian form
(13) 
which then yields the hot Coulomb coefficients
(14) 
and
(15) 
where 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 . This means that the Coulomb energy loss rate should also be transformed to or
(16) 
Note that at the effective Coulomb loss timescale diverges and for the rate which means that the test particle actually gains energy from the target particles. Figure 1 shows the effective Coulomb loss timescales, defined as , as a function of for the four different pairs of the test and target particles (ee, pp, ep, pe).^{3}^{3}3Inclusion of the effects of particles with does not change and but increases and by 1.16. Adding losses to particles reduces the two latter times by 8% (to 1.07) at relativistic energies and by 10 to 30% around . These effects are included in Figure 1. There are several features in these equation that deserve some discussion.

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 .

These forms of the Coulomb rates also satisfy the timeindependent 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.

As also can be deduced from the above equations the shapes of the timescales shown in Figure 1 are invariant and independent of but they scale with the temperature of the target particles as (see Appendix C). Thus, although for the purpose of comparison with the acceleration timescale we have used keV, the loss times as scaled in the vertical axis are valid for all temperatures.

Finally the relative values of the timescales play an important role. When electrons are energized the acceleration timescale should be comparable to the shortest ee collision loss time. In this case the electrons gain energy but share very little of it with protons since ep timescale is times longer. But as electrons are heated the ee timescale increases as and when some of the energy goes into protons. This is why EP08, 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 (pp time is shorter than pe) 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 pe interactions become important and some of the energy goes to electrons which are thermalized quickly because of much shorter ee 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
The algorithm we adapt iteratively solves the coupled FP equations of protons and electrons interacting via Coulomb collisions with only proton subject to SA. The coupling arises from the fact that the energy lost to electrons by protons gives the heating rate of electrons and vice versa. At each time iteration step the algorithm uses the Coulomb coefficients calculated based on particle distributions obtained in the previous time step. PE08 fitted the distributions in each step to a Maxwellian with “effective temperatures” and used the hot plasma Equations (14) and (15) to calculate the Coulomb collision coefficients. This clearly is a reasonable approximation when the distributions are nearly isothermal and saves considerable computing time. However, for particle distributions with strong nonthermal tails that are of interest here this may not be a good approximation. Moreover, when two species are involved it becomes more complicated to determine their respective effective temperatures. For these reasons, we use the more accurate but the more time consuming algorithm whereby at each time step we calculate the Coulomb collision coefficients by integrating the more exact Equations (9) and (12) over the particle distributions determined from solution of the coupled FP equation at the previous time step. Several tests of the integrity of this algorithm are presented in Appendix C.
3 Time Evolution of Proton and Electron Spectra
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 Maxwellian distribution of , and interact both internally and mutually via Coulomb collisions while protons undergo SA by turbulence. As explained in §2, the acceleration model is specified by the choice of parameters , , and . We give in units of , and is expressed in terms of . In addition the other relevant times and in the plateau region; for the assumed initial temperature keV the plateau is at (see Fig. 1).
In what follows we show time evolution of proton spectra, for several values of the acceleration parameters and , at evenly spaced time steps, beginning with initial Maxwellian distributions with keV, which shifts monotonically to higher temperature and develops an increasing nonthermal component resembling the socalled distributions. Electron spectra, which tend to be purely Maxwellian because of the short ee collision times, are shown only in two interesting cases. These spectra are fit to a thermal component of temperature plus an excess which we call the nonthermal component. We give the the temperature at the final step, and using the method described in the Appendix B of PE08, we also give the fraction of the number of protons in the thermal component , the ratio of the nonthermal to total energy , and an approximate powerlaw 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 , 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 and (left) and (right). In both models (and others with similar values of ) significant deviations from an isothermal distribution appears after the plasma temperature is increased by several decades, which happens after yr.
In what follows we focus on and show how the time evolution varies with in Figure 3 and with in in Figure 4. For the first set of models with considerably longer than thermalizing timescales of protons [see Equations (C2) and (C4)], 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 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 .^{4}^{4}4The fact that thermalization timescales are increasing functions of energy also accounts for the latetime emergence of the considerable nonthermal component in the models. The second set of models with a lower value of exhibit 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 . This behavior is due to the fact that for . On the other hand, when (or , to be more precise) is comparable to or shorter than , which is necessary for development of a distinct nonthermal tail, one can avoid excessive heating.
Two examples that have large but comparable to or shorter than are shown in Figure 5. The left panel shows the later time evolution of the proton spectrum for the model , , and . Although the spectrum at = exhibits a promising nonthermal tail and temperature that is only few times larger than , its shape changes rapidly as a considerable fraction of the protons are accelerated beyond . By the time = , the spectrum is completely nonthermal with a very small thermal component. On the right panel we show the evolution of the proton spectrum for the model , , . In this case the spectrum remains nearly thermal until = , then begins to exhibit an appreciable nonthermal tail at = , 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 and then gets very rapidly accelerated once a considerable fraction of protons obtains kinetic energy larger than , leading to a runaway spectrum.
These results imply that, in order to get a nonthermal tail by an acceleration model with , small and large , a great deal of finetuning in the duration of acceleration would be necessary. Also, we may see from these results that 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 . The above conclusions will also hold for acceleration models with , whose energy dependence at lowenergy (i.e. ) is even steeper than models.
As expected from the discussion in §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 pe 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 for the same models shown in Figure 5. As expected, in some cases and later times (when ) the electron temperature 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 Xrays by a proton rather than an electron nonthermal tail.
In summery, 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 Xrays in the ICM over a Hubble time would be difficult.
4 Summary and Discussion
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 Eq. (4) in §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 like in solar and stellar flares and in the ICM of galaxy clusters and could produce nonthermal bremsstrahlung radiation in the hard Xray regime in collisions with the background electrons.
We derived the Coulomb energyloss 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 formulas for Maxwellian distribution that asymptotically approach the cold target relations at high energies. We have generalized the algorithm developed by PE08 to solve the coupled electronproton 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 absence of acceleration the particle distributions relax to Maxwellian within the expected timescales starting from several nonequilibrium initial conditions. In Appendix D we discuss the inelastic energy loss rates and present a new analytic formula for proton energy loss rate due to pion production.
Using the algorithm, we obtain the time evolution of the proton and electron distributions. Because of shorter collision time of electrons the production of nonthermal electron tails require 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:

For , even when acceleration time is sufficiently small so that the energizing rate is comparable to or faster than the thermalization rate, an initially Maxwellian proton spectrum evolves into a hotter one and broader than thermal distribution with little sign of a distinct nonthermal tail. This is because for the energizing rate is smaller at higher energies and not efficient in acceleration into a highenergy tail. This suggests that must be positive for production of a distinct nonthermal tail and a mild temperature increase. This means we need turbulence with spectrum considerably steeper than Kolmogorov.

We show that for the temperature increases more slowly as long as most of the particles have kinetic energy below . In addition, a significant nonthermal tail may appear if is comparable to or smaller than the proton thermalization time. For such small values of , however, once a considerable fraction of the protons have reached , they rapidly evolves into a completely nonthermal, runaway distribution. This will be more pronounced for even higher values of , where acceleration rate increases more rapidly with energy.

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 that of the proton. But when the proton distribution contains a significant or a larger nonthermal component (and 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 shows that a great deal of finetuning and a steeper than Kolmogorov spectrum of turbulence is necessary for creating 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 (Wolfe & Melia 2006: Boldt & Serlemitsos 1969), nonthermal inverse bremsstrahlung cannot be a significant contributor to hard Xray emissions (see also Emslie & Brown 1985).
Acknowledgements
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.
Appendix A Proof of Steady State Maxwellian
When two species of particles denoted by and are in thermal equilibrium, the timeindependent FP equation ( ) for each of the two of the species, say , is written as
(A1) 
where is the (Maxwellian) distribution of the species and all the Coulomb coefficients are evaluated at the same temperature as that of . Since the must be timeindependent even in the absence of Coulomb interactions with the other species of particles, i.e.,
(A2) 
it follows that, in particular,
(A3) 
Given and we can integrate this equation to uniquely determine . The diffusion coefficient uniquely determined in this way agrees exactly with the ones presented in §2, which justifies the assumptions used in their derivations.
Appendix B Energy Conservation
Here, we show that, for arbitrary particle distributions and FP coefficients are given by Equation (9) and (12), our algorithm described in Section 2 satisfies the total energy conservation. Let be the total energy of the two species of particles. Then, from their FP equations with noflux boundary condition (see Equation (3) of Park & Petrosian 1995), we find
(B1) 
where
(B2) 
(B3) 
(B4) 
(B5) 
A and B separately vanish because of antisymmetry of and , and because , it follows that . This proves the claim.
Appendix C Thermalization Tests
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 the 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
(C1)  
(C2)  
(C3)  
(C4) 
where and are defined as of the total energy of the initial distributions.^{5}^{5}5In the last expressions numerical factors of order unity are omitted and for and , appropriate asymptotic forms of the function are used. Interestingly as can be seen from Figure 1, for the interspecies interaction times and are almost equal and are and times longer than and , 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 = 25 keV and = 5 keV. As expected, both species remain Maxwellian throughout and their temperatures approach the average of the two initial temperatures over several ’s and ’s, respectively calculated at the initial temperature. We get similar results with initial values reversed; = 5 keV and = 25 keV.^{6}^{6}6Figures showing this results, and several others not shown in this paper, can be found in Byungwoo Kang’s (2013) Senior Thesis at Stanford. Next, we consider the thermalization of the protons and electrons when they initially have a narrow lognormal distributions. The normalized lognormal distribution function is given by
(C5) 
whose mean is . Note that, for , , and the lognormal distribution becomes essentially a Dirac delta function centered at . In this case, the proton and electron distributions should approach the Maxwellian shape over several ’s and ’s respectively. At the same time, their temperatures should also converge to over several ’s and ’s respectively. In Figure 8 we show the time evolution of the proton and electron distributions, when the protons and electrons are initially in lognormal distributions with = 5 keV, = 0.02, and = 25 keV, = 0.02, respectively. Not shown here, within 10’s of and , respectively, electrons and protons acquire Maxwellian distributions (as can be seen in this figure this is the case at and ,respectively. However it takes several tens of these timescale before they both reach the final equilibrium of keV. Similar results are obtained with initially more energetic protons. Similar results are obtained starting with 25 keV protons and 5 keV electrons.
Appendix D Inelastic Energy Loss Rates and Times
The inelastic energy loss rates for electrons are due to well known processes of synchrotron, inverse Compton (with , where is the total energy density of magnetic field and soft photons) and electron proton bremsstrahlung with
(D1) 
where is the fine structure constant and is a slowly varying function. At relativistic energies this rate should be increased by nearly an equal amount due to electronelectron bremsstrahlung. Thus, bremsstrahlung loss becomes dominant compared to elastic scattering at Lorentz factors . As shown in Figure 9, adding these rate modifies the electron energy time scales shown in Figure 1 at high energies.
For protons the main inelastic loss rate is due to hadronic interactions. For MeV protons the main loss is due to the protonion interactions producing ionic deexcitation lines in the 17 MeV range. Loss rate of to MeV protons is dominated by neutron production. Pion production starts at threshold energy MeV and dominates all losses (including Coulomb) above few GeV. Protonelectron interactions can produce X and gammarays via (inverse) bremsstrahlung process. Since energy loss rate is an invariant quantity it follows that will be equal to that of ep bremsstrahlung in the rest frame of the proton, so that this rate is also given by Eq. (D1) with electron velocity and Lorentz factor replaced with that of protons. The inverse bremsstrahlung loss becomes negligible for proton velocities less than the electron thermal velocity, similar to the Coulomb loss (see Figure 1).
One can get a simple approximate expression for the energy loss rate due to hadronic interactions at relativistic energies but there are no simple furmulae for semirelativistic energies and at energies near the threshold. We combine previously published results on these interaction and derive a simple formulas in three energy regimes with few percent accuracy at the boundaries of these regimes.
d.1
For this very high energy regime, we make use of the work of Kelner et al. (2006) on the energy spectra of the secondaries particles produced by interactions. If we denote the number of secondaries in the kinetic energy interval produced by a proton of energy (in unit of ) by , then the energy loss rate can be written as
(D2) 
where the index is summed over all secondaries and is the inelastic interaction cross section that encodes information as to how many inelastic interactions occur (not the number of the secondaries produced, as some cross sections defined elsewhere do).
There are two choices as to counting the secondaries. We can count the relatively stable secondaries, which are gammarays, 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 rate calculated in both ways agree with each other within few percent, and can be approximated as
(D3) 
d.2
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 (D2) with
(D4) 
where , , and
(D5) 
where , with experimentally derived values of , , , , and . In these equations ; all quantities measured in the c.m.s. Using these equations we find an energy loss rate
d.3
For the lowest energy range down to the threshold kinetic energy we use the Stecker’s model (Stecker 1970; also see Dermer 1986), where we can approximate the production process by the socalled ’isobarplusfireball model’.^{7}^{7}7Stephens and Badhwar (1981) argue that their representation works quite well in this energy range too, despite the fact that the energy range of the experimental data used to determine these constants was around a few ten GeV. However, Dermer (1986) shows that, near the threshold energy, Stecker’s model fits better with experimental pion energy spectra than StephensBadhwar model. In this model, we will assume that all pions are produced in two ways: via intermediate production and decay of the nonstrange 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 MaxwellBoltzmann type distribution. In particular, for , 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 carry momentum either directly forward or directly backward in the c.m.s, and that the lowenergy isobar production process produces pions through the twostage decay
The normalized production spectrum of pions, to be used in Equation (D2), in this model is given by
(D6) 
where
(D7) 
and the normalized isobar mass spectrum is given by the BreitWigner distribution with the normalization factor . However, since is sufficiently small, we can approximate by the Diracdelta function and therefore ignore the normalization factor .
Here if and 0 otherwise, and is the pion Lorentz factor in the lab system (LS). The Lorentz factors of the forward (+) and backward () moving isobars are , where is the Lorentz factor of the center of mass (CM) with respect to the LS, and is the Lorentz factor of the isobar in the CM. The pion Lorentz factor in the rest frame of the isobar is .
Again with the help of these equations and Equation (D2) we obtain an energy loss rate of .
Now putting all this together we obtain energy loss rate as a function of proton energy (in units of ) as