Semi-analytical model of cosmic ray electron transport

We present a numerical extension to the analytical propagation model introduced in Hein and Spanier (2008) to describe the leptonic population in the galactic disc. The model is used to derive a possible identification of the components that contribute to the leptonic cosmic ray spectrum, as measured by PAMELA, Fermi and HESS, with an emphasis on secondary e+ - e- production in collisions of cosmic ray particles with ambient interstellar medium (ISM). We find that besides secondaries, an additional source symmetric in e+ and e- production is needed to explain both the PAMELA anomaly and the Fermi bump, assuming a power-law primary electron spectrum. Our model also allows us to derive constraints for some properties of the ISM.


Introduction
Recently the leptonic component of the cosmic ray spectrum has gained new attention. New observations from ATIC, PAMELA and Fermi show a deviation from a power-law in the form of an excess in both the electron and positron spectra. Annihilating dark matter (Allahverdi et al., 2009) and nearby pulsars (Büsching et al., 2008;Grimani, 2009;Blasi and Amato, 2011) (among other hypotheses) have been proposed as possible sources of the excess leptons. Regardless of the source, a propagation model is needed to connect the energy spectrum measured on Earth with the injection spectra.
We present our numerical cosmic ray transport model in application to the high energy electron transport in the ISM. Spatial and momentum diffusion, particle escape, acceleration via Fermi I and continuous energy losses were taken into account and their effects on the steady-state energy spectrum analyzed. In solving the transport equation we em-Correspondence to: A. Ivascenko (aivascenko@astro.uni-wuerzburg.de) ployed quasi-linear transport theory, the diffusion approximation and a separation of the spatial and momentum problem to obtain the leaky-box-equation, which was then solved numerically. The spatial problem was solved analytically in cylindrical and prolate spheroidal coordinates.
The transport model was employed to calculate the spectrum of secondary electrons in our galaxy. We assume the leptons from pioArXivn decay to be the dominating component of the secondary spectrum. The lepton spectrum from pion collisions of highly relativistic cosmic ray protons with thermal protons in the ISM was calculated using the parametrization of pion production in p-p-collisions from Kelner et al. (2006) and used as the injection spectrum for the transport model. With realistic simulation parameters the resulting positron flux in the vicinity of the solar system lies remarkably close to the low-energy PAMELA data. Assuming a generic power-law primary e − spectrum and an additional e + − e − source allows us to fit the datasets from PAMELA, Fermi and HESS very nicely and to put constraints on several transport model parameters and the corresponding ISM properties.

The data
The leptonic cosmic ray spectrum gained new attention when the results of the balloon experiment ATIC were published (Chang et al., 2008). The authors claimed an anomalous excess in the spectrum just below 1 TeV. This started wild speculations about the nature of the source for these leptons, which should be in the relative vicinity of the solar system.
Those speculations were fueled further by new measurements from the PAMELA satellite published in Adriani et al. (2009) claiming a rise in the positron spectrum above 10 GeV. The lepton spectrum from the Fermi satellite (Abdo et al., 2009) could not confirm the ATIC-peak quantitatively but 2 A. Ivascenko and F. Spanier: CR electron transport model showed a deviation from a straight power-law in the same energy region.
The HESS atmospheric Cerenkov telescope extended the spectrum measurement beyond 1 TeV showing a cut-off or a steeper power-law but could also not confirm the ATIC excess (Aharonian et al., 2009).
Recently the Fermi lepton spectrum was extended to lower energies (7 GeV) in Ackermann et al. (2010). Also the PAMELA results were updated with new experimental data and different evaluation techniques in Adriani et al. (2010).
To interpret these measurements, it is necessary to connect the injection spectra at the (possible) sources with the observations in the solar system. To do that we need a transport model that includes all the relevant processes in the ISM (spatial and momentum diffusion, cooling, particle escape, etc.).

Transport model
The derivation of the transport equation follows closely Lerche and Schlickeiser (1988) and Hein and Spanier (2008). Using Quasi Linear Theory and the diffusion approximation the Vlasov equation is transformed into the steady-state diffusion-convection equation for the phase space density f (x,p), which can be written in the form with the spatial operator and the momentum operator . (3) This equation describes the diffusion (coefficients κ and D in space and momentum, respectively), acceleration and cooling (ratesṗ) and escape (timescale τ ) of an injected source distribution S. Following Lerche and Schlickeiser (1988) the equation can be solved by separating the operators and the source term in their spatial and momentum dependencies. Then the resulting steady state particle distribution can be written as an infinite sum: The spatial coefficients A i (x) can be acquired by solving the diffusion equation in the appropriate geometry where ω 2 i take the role of inverse escape times for the momentum modes and u is a convolution variable.
In this work we used the analytical solution for cylindrical geometries as provided by Wang and Schlickeiser (1987) since this offers the best reflection of the symmetries of a spiral galaxy. Analytical solutions in spherical and prolate spheroidal coordinates are also available Hein and Spanier, 2008) and can be used to describe particle transport in elliptical galaxies and galaxy clusters.
The momentum modes R i (p) are solutions of the ODE, commonly referred to as the leaky-box equation: The parameters a 1 and a 2 represent the absolute scales for Fermi I and II processes and depend mainly on the Alfven and shock speed, respectively. As Lerche and Schlickeiser (1988) have shown, an analytical solution for the equation can be found, if the cooling rateṗ scales linearly with the particle energy, which is a good assumption for high energy protons and heavier nuclei. Electrons above a few GeV, however, are cooled primarily by synchrotron radiation and IC scattering, both scaling quadratically in particle energy. Therefore the leaky-box equation was solved using numeric relaxation methods.
To summarize, the model describes a particle distribution that is injected by a homogeneous population of timeindependent sources in the galactic disc (injection region), diffuses throughout the galactic disc and halo (confinement region), cooling by bremsstrahlung, adiabatic deceleration, synchrotron radiation and inverse Compton scattering, and leaves the confinement region at some point. The matter, photon and magnetic field distributions are assumed to be homogeneous over the whole confinement region.
A model with a uniform source distribution is not a good choice to simulate high energy leptons from discrete sources (SNR, PWN. . . ), since the short energy loss timescales make the source distribution critical for the resulting spectra, as discussed in Blasi and Amato (2011). Conversely, production of secondaries from p-p collisions in the galactic halo can be approximated very well as a homogeneous source, since there are no dramatic variations in the CR and ISM densities.
A similar analysis by Moskalenko and Strong (1998) used an analogue model as in this paper, but obviously did not include the new data above 1 GeV. Stephens (2001) has basically the same drawbacks, in addition this work has a more sophisticated secondary production model. Compared with the recent works in this area like Blasi and Amato (2011), we concentrate more on the low energy end of the current PAMELA data and the constraints it can put on propagation parameters if interpreted as secondaries.

Secondary leptons
Fortunately there is a lepton production process that can be approximated with a homogeneous spatial distribution: the secondary leptons from the decay of products of proton collisions. To get a galaxy-averaged flux of secondaries we consider collisions of high energy cosmic ray protons with ambient thermal matter, assuming that both populations are homogeneously distributed with average densities.
The pion production cross-section and the electron source function for a single p-p collision used here are parametrizations of results from the SIBYLL event generator from Kelner et al. (2006). Folding the source function F e Ee Ep ,E e with the cosmic ray spectrum J p (E p ) and the inelastic collision cross-section σ inel (E p ) yields the electron flux: Since the source function is sharply peaked, the CR spectrum can be cut off exponentially at the first knee (≈ 1000 TeV) with negligible deviations to the lepton spectrum below 10 TeV: This parametrisation of the CR spectrum was obtained by fitting the CR data between 10 10 and 10 14 eV in Cronin et al. (1997). The average thermal matter density n H is treated as a free parameter. The resulting electron flux with a slightly harder power-law (spectral index 2.62) is used as the injection function in the transport model. Since the parametrizations in Kelner et al. (2006) only consider p-p interactions, the actual flux of secondaries is underestimated by the contribution from heavier nuclei.

Parameter constraints
The principal shape of the resulting steady-state lepton spectrum is shown in Fig. 1. The cooling, diffusion and escape processes introduce breaks in the injected power-law, so that multiple energy intervals with different spectral indices are formed, each of them defined by a different process. At the first break the injected power-law is steepened by 1/3 as the diffusion/escape time scale becomes shorter than the time scale of bremsstrahlung and adiabatic cooling (linear in energy). The second break appears when the time scale for synchrotron and IC losses becomes shorter that the diffusion time scale. These quadratic losses steepen the spectrum by 1. The absolute flux value is determined mainly by the ISM density, since the injected secondary spectrum has a linear n H dependence (Eq. (7)), and to a lesser extent by escape and cooling losses, which transport particles out of the simulation region. So the model has three major parameters: the linear energy loss time scale that reflects the thermal gas density n H , the quadratic loss time scale that incorporates the energy densities of the EM field (composed of a constant CMB density of ≈ 0.3 eV/cm 3 and a variable magnetic field B) and the escape time scale t esc that corresponds to the size of the confinement region.
Measurements of the positron spectrum, in particular the newest results from PAMELA, allow us to constraint those parameters. A rather conservative statement can be made, if we assume that the secondaries have a negligible contribution to the positron spectrum. As shown in Fig. 1, assuming an average ISM density of ≈ 10 3 m −3 , the confinement region has to be less than the assumed 30 kpc or the diffusion process has to be more efficient, so that the particles leave the galaxy 10 times faster (black dotted line). A higher ISM density would require an even more efficient particle escape, so as not to overshoot the measured positron spectrum. Since the spectral shape wouldn't play any role, the magnetic field can be almost arbitrary.
On the other hand, the PAMELA data can be fitted very well with a superposition of two power-laws with a spectral index of about 3.6 for the low energy part, corresponding remarkably well with the secondary spectrum steepened by synchrotron and IC losses. Assuming the low energy positrons to be secondaries puts much stricter constraints on the parameters, since the second break in the spectrum has to appear at ≈ 1 GeV to match the data. The break energy is determined by the equilibrium of the diffusive/escape losses and the radiative losses. The thick black line in Fig. 1 is the best fit with n H = 1.5 · 10 3 m −3 , B = 3µG and t = t esc .  Fig. 2. A fit of Fermi (Ackermann et al., 2010) and HESS (Aharonian et al., 2009) lepton data (solid red line) and PAMELA (Adriani et al., 2010) positron data (solid green line). Grey areas show the systematic errors. The components that contribute to the spectra are secondary positrons (dashed green), a generic power-law primary spectrum (dashed red) and an additional symmetric e + − e − source (solid and dashed black).
The radiative losses in the galactic halo are dominated by IC scattering on the CMB photons, so that lowering the magnetic field has no effect on the break energy. A slight flattening at the low energy end of the spectrum is needed for the best fit of the new PAMELA data. This corresponds to a halo radius of 30 kpc and a spatial diffusion coefficient κ(E = 1GeV) = 4.5 · 10 24 m 2 s −1 , yielding an average confinement time t esc = 2 · 10 17 s. The absolute flux fixes the third parameter, the ISM density, to 1.5 · 10 3 m −3 . Taking into consideration that about 99% of the confinement region represent the galactic halo with an average gas density of about 10 3 m −3 (McKee and Ostriker, 1977), the parameters seem very realistic.

Lepton excess
In Fig. 2 the latest Fermi lepton measurement, in which the spectrum has been extended down to 7 GeV, is shown. A generic fit of it was used to calculate a positron flux from the PAMELA positron fraction measurement. This data is used because it was taken simultaneously in the same solar cycle making interpretation easier. The HESS high-energy data was shifted down by 15% to better coincide with the Fermi data. This is still within the systematic error margins of HESS, not to mention the systematic error of Fermi being of the same order of magnitude. Once the secondary background has been subtracted from the PAMELA positron flux, the excess can be fitted very well with a simple power-law with a spectral index of 2.3 (Fig. 2, black dashed line). Remarkably, if we assume the source of these "extra" positrons to be symmetric in e + − e − , we can also fit the Fermi bump leaving a primary lepton component of the shape E −3.2 · exp(−E/1.46 TeV) (red dashed line). The (shifted) high-energy HESS data provides a means to fix the cut-off energy of the "extra" leptons to about 920 GeV.

Conclusions
With the help of our semi-analytic transport model we were able to identify the low-energy PAMELA data as secondary positrons from collisions of CR protons with ISM gas using very realistic values for the simulation parameters. An indirect implication from that is that lepton transport is dominated by radiative losses (synchrotron and IC) above 1 GeV.
Assuming a generic primary lepton spectrum of the form E −3.2 ·exp(−E/1.46 TeV), the PAMELA anomaly, the Fermi bump and the high energy HESS data can all be fitted by a single component, that is symmetric in e + − e − with a power-law spectral index of 2.3 and a cut-off energy of 920 GeV. This could either be a lepton population injected into the ISM by nearby source with E −2.3 , so that the propagation time scale is shorter than the radiative loss time scale, or a more distant source with a much harder spectrum E −1.3 which is then steepened by radiative losses. A supernova remnant would be a good candidate in the first case, pulsars of different ages would suite both cases.