Pulsar interpretation of the lepton spectra measured by AMS-02

AMS-02 recently published its lepton spectra measurement. The results show that the positron fraction no longer increases above $\sim$200 GeV. The aim of this work is to investigate the possibility that the excess of positron fraction is due to pulsars. Nearby known pulsars from ATNF catalogue are considered as a possible primary positron source of the high energy positrons. We find that the pulsars with age $T\simeq (0.45\sim4.5)\times10^{5}$ yr and distance $d<0.5$ kpc can explain the behavior of positron fraction of AMS-02 in the range of high energy. We show that each of the four pulsars --- Geminga, J1741-2054, Monogem and J0942-5552 --- is able to be a single source satisfying all considered physical requirements. We also discuss the possibility that these high energy $e^{\pm}$ are from multiple pulsars. The multiple pulsars contribution predicts a positron fraction with some structures at higher energies.


Introduction
The positron fraction spectrum e + /(e + +e − ) in the cosmic ray (CR) contains two components: secondary e ± produced by nuclei collision and primary e − . It is currently believed that these two components, each of which will produce a diffused power low spectrum, predict a positron fraction which goes down with energy. However, the latest results measured by the Alpha Magnetic Spectrometer (AMS-02) with high accuracy indicate that the positron fraction increases with energy above ∼8 GeV and does not increase with energy above ∼200 GeV [1,2]. This behavior, which is also observed by the payload for antimatter matter exploration and light-nuclei astrophysics (PAMELA) [3,54,53] and the Fermi Large Area Telescope (Fermi-LAT) [4,5], is not compatible with only diffused power low components. The "cutoff" behavior above 200 GeV indicates that potential sources produce the exceed of electron and positron pairs. AMS-02 [1,2] is a state-of-the-art astroparticle detector installed on the International Space Station (ISS). It carries a Transition Radiation Detector (TRD) and a Electromagnetic Calorimeter (ECAL). These two sub-detectors provide independent proton/lepton identification, which will achieve a much larger proton rejection power of AMS-02 compared with PAMELA which has only one Electromagnetic Calorimeter for proton/lepton identification. Compared with Fermi-LAT, AMS-02 has a large magnet which can identify charge sign of the particle. Thus, the contamination of electrons (also called "charge confusion" in [1]) in the positron sample of AMS-02 is much smaller that that of Fermi-LAT. For the reasons given above, there is much less proton or charge confusion contamination in AMS-02 measurement than that in PAMELA or Fermi-LAT. Here, we only interpret AMS-02 result due to the lack of knowledge of the contamination control in PAMELA and Fermi-LAT measurements.
The AMS-02's recent measurements of positron fraction [1], e + flux and e − flux [2] were published. The e − flux contains three components: primary e − , secondary e − and e − from unknown sources. The e + flux contains only two components: e + from secondary production and that from unknown sources. To avoid the unnecessary uncertainty of primary e − , the e + flux seems to be an ideal spectrum to study extra sources. However, there is an acceptance uncertainty from the detector itself in the e + and e − fluxes. This uncertainty in e + flux is strongly correlated with that in e − flux [2], especially at high energies. The positron fraction can avoid this systematic uncertainty [1]. For example, one can clearly see a drop at the last point (350 GeV ∼ 500 GeV) in the positron fraction but cannot tell a drop at the last point (370 GeV ∼ 500 GeV) in the e + flux due to its larger error bars. Therefore, positron fraction is used to study extra sources while e − flux is used to estimate the primary e − which will affect the denominator of e + /(e + +e − ).
Recent studies have proposed some interpretations, such as dark matter annihilation or decay [6,7,8,9,10,11,12], supernova remnants (SNRs) [13,14,15,16,17,18], secondary production in the interstellar medium (ISM) [19] and pulsars [20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36]. Cosmic ray flux data can also be together with other observations (like the dark matter relic density and the direct detection experimental results etc.) to give a combined constraint on dark matter models [37,38]. Besides dark matter scenario, the others can provide astrophysical explanations which do not require the existence of unknown particles. SNRs models, for instance in [17], predict a spectrum increase steadily above 500 GeV. The "model-independent" approach from [19], sets an upper limit of the positron fraction by neglecting radiative losses of electrons and positron but does not indicate any obvious cutoff in the spectrum. Among them, the pulsar interpretation is one of the scenarios which predict a cutoff at a few hundred GeV in the positron fraction spectrum.
A pulsar is widely regarded as a rotating neutron star with a strong magnetosphere, which can accelerate electrons, which will induce an electromagnetic cascade through the emission of curvature radiation [39,40,41,42]. This leads to the production of high energy photons which eventually induces pair production. This process produces the same amount of high energy e + and e − , which can escape from the magnetosphere and propagate to the earth. There is a cutoff energy of the photons produced in a pulsar, which leads to a cutoff in the positron fraction.
In this paper, we use DRAGON [43, 44,45] as a numerical approach to model the propagation environment, to tune the related parameters and to estimate the e ± background. ROOT is used to minimize χ 2 to get the best fit results. We consider six nearby pulsars from ATNF catalogue [46,47] as the possible extra single sources of the high energy positrons. We find only four, which are Geminga, J1741-2054, Monogem and J0942-5552, can survive from all considered physical requirements. We then discuss the possibility that these high energy e ± are from multiple pulsars. The multiple pulsars contribution predicts a positron fraction with some structures at higher energies.
The paper is organized as follows. Section 2 shows the way where e ± background is estimated. In Section 3 we review the properties of pulsars and derive the profile of e ± fluxes produced by a pulsar. We try to use single pulsars in Section 4 and multiple pulsars in Section 5 to interpret positron fraction. The conclusions are drawn in Section 6. In Appendix A, we show that the diffusion energy-loss equation for a burst-like source can be solved without using the spherically symmetric approximation.

Propagation parameters and e ± background
The Galacitc background of the lepton fluxes are considered as three main components, which are primary electrons from CR sources, secondary electrons and positrons from the interactions between the CR and the ISM. The propagation of e ± in the Galaxy obeys the following equation [44,48] where p ≡ | p| is the particle momentum; f i ( x, t, p) is the particle number density of a species i per unit momentum interval; v c is the convection velocity, β ≡ v/c is the ratio of velocity to the speed of light; σ in is the total inelastic cross section onto the ISM gas, whose density is n gas ; σ ji is the production cross section of the species i by the fragmentation of the species j (with j > i); and Q i ( x, t, p) is the source term of species i, which can be thought to be steady Q i = Q i ( x, p) for background CR particles.
The spatial diffusion coefficient D in the cylindrical coordinate system (r, z) may be parameterized as [44,48,49] D(ρ, r, z) = D(ρ, r)e |z|/zt and D(ρ, r) where ρ ≡ pc/(Ze) is defined as the particle magnetic rigidity, z t is the scale height of the diffusion coefficient, and δ is the index of the power-law dependence of the diffusion coefficient on the rigidity. D 0 is the normalization of the diffusion coefficient at the reference rigidity ρ 0 = 4 GV. The function f (r) describes a possible radial dependence of D, and it can be taken to be unity for simplicity. The diffusion coefficient in momentum space D pp is related to the spatial diffusion coefficient D by [11,49,50] where v A is the Alfven velocity, and δ is the power-law index as given in (2). w is the ratio of magnetohydrodynamic wave energy density to the magnetic field energy density, and it is usually taken to be 1. DRAGON [43,44,45] is used to tune the propagation parameters according to the B/C ratio, since it is sensitive to the parameters. The parameters are as follows: the diffusion coefficient D 0 ρ 0 =4 GV = 3.95 × 10 28 cm 2 s −1 , the diffusion coefficient's power-law index δ = 0.4, the scale height of the diffusion coefficient z t = 4.04 kpc, the halo size L = 4 kpc, and the Alfven velocity v A = 10 km/s. In order to match the proton spectrum, three breaks, which are 6.7 GV, 11 GV and 300 GV, are introduced in the injection spectrum of nuclei. The spectral indices below and above the breaks are 2.35, 2.45, 2.55 and 2.47, respectively. The Ferriere model [51] is used as the source distribution for the primary components, e.g. SNRs. The reaccelaration model is from Ref. [52]. To assure that the propagation parameters are correct, we need to compare the model prediction with the Boron-to-Carbon ratio [53] and the proton flux [54] from PAMELA. As shown in Fig. 1, the set of parameters we taken can reproduce the Boron-to-Carbon ratio and the proton flux well. According to this set of parameters, we can obtain the fluxes of the secondary positrons and electrons.
A power-law spectrum with two breaks is introduced to parameterize the injection spectrum of the primary electrons as a function of rigidity,

e ± from a single pulsar
The pulsars are potential sources which could produce extra e ± at high energy [24,25,26,27,33]. Electrons can be accelerated by the strong magnetosphere of the pulsars, and this acceleration produces photons. When those photons annihilate with each other, they can produce e ± pairs. Thus, the e ± energies are related to the pulsar magnetosphere. Assuming the pulsar magnetosphere as a magnetic dipole, this magnetic dipole radiation energy is proportional to the spin down luminosity. Due to this spin down (i.e. slowing of rotation), the rotational frequency of a pulsar Ω ≡ 2π/P (with P being the period) is a function of time as follows [24,26,33] where Ω 0 is the initial spin frequency of the pulsar and τ 0 is a time scale which describes the spin-down luminosity decays. τ 0 cannot be directly obtained from pulsar timing observations, and it is assumed to be [26,33] τ 0 10 4 yr (6) The rotational energy of the pulsar is E(t) = (1/2)IΩ 2 (t). Here I is the moment of inertia, which is related to the mass and the radius of the pulsar and can be regarded as a time independent value. The magnetic dipole radiation energy is equal to the energy loss rate, The total energy loss of a pulsar is [26,27,33] The total energy injection of e ± out of a pulsar should be proportional to the total energy loss where η is the efficiency of the injected e ± energy converted from the magnetic dipole radiation energy. The pulsar characteristic age is defined as [47] T For a mature pulsar with t τ 0 , we have T t. In this condition, eqs. (7), (8) and (9) become The propagation equation for the e ± can be described as [24,33] ∂f where f ( x, t, E) is the number density per unit energy interval of e ± ; D(E) = (v/c)D 0 (E/4 GeV) δ is the diffusion coefficient with the velosity v of the particle, the speed c of light, D 0 and δ the same as the parameters used to calculate the background in Section 2; and b( the rate of energy loss due to inverse Compton scattering and synchrotron [4,26,33]. The source term Q( x, t, E) of a pulsar can be described by a burst-like source with a power-law energy spectrum and an exponential cutoff where Q 0 is the normalization factor related to the total injected energy E out , α is the spectral index, and E cut is the cutoff energy. Many previous works (for example, Refs. [21,33]) use the spherically symmetric approximation to get the solution of the equation (14) with the source (15). However, this approximation is actually unnecessary. In Appendix A, we obtain the exact solution of (14) without using the spherically symmetric approximation, and our result agrees with the previous works. Using the results, i.e. eqs. (54) and (58), in Appendix A, we obtain the electron or positron flux observed at the earth as follows: where d is the distance between the earth and the source, the diffusion distance r dif is given by and the diffusion time t dif is the time the charged particle travels in the ISM before it reaches the earth. We simply assume that t dif T . The maxium energy E max is defined as E max = 1/(b 0 T ). The positron fraction from AMS-02 implies an extra source with an E max (500 ∼ 5000) GeV, which corresponds to a pulsar with an age T (0.45 ∼ 4.5) × 10 5 yr. Note that the distance a particle travels in the ISM should be larger than the distance between the earth and the source, otherwise it cannot reach the earth. Thus, we have r dif > d. Looking into the r dif as a function of diffusion time t dif and lepton energy E, we can get the upper limit of the pulsar distance d, as Fig. 2 shows. Given the diffusion time T

Single pulsar interpretation
We give a few simple examples to explain high energy positron fraction of AMS-02 [1] using single pulsars. The background electrons and positrons are described in Section 2. The primary electron flux is scaled by a normalization factor A prim,e − since we can not constrain the electron flux contribution from SNRs. We take the age T and the distance d from the ATNF catalogue and fit the positron fraction to obtain the free parameters in (16), the spectral index α, and the normalization Q 0 . Q 0 is fixed by the relation [18,27,32] , which approximately yields Q 0 E out for α 2. We set the cutoff energy E cut =5000 GeV, which is large enough, as it does not change the shape of pulsar contribution.
We use six nearby single pulsars to fit the positron fraction. Minuit package in ROOT is used to determine the parameters to minimize χ 2 . The best results of the single pulsars are listed in Table 1. The results are also shown in Fig. 3 Table 1: Parameters of six nearby single pulsars from the best fit results. The χ 2 /ndf from the fits of Geminga, J1741-2054, Monogem, J0942-5552 and J1001-5507 are smaller than 1, which show a good agreement between those single pulsar models and the experiment.
parameters of the best fit results, the positron fraction can be well reproduced by these single pulsar's contributions. Table 1 tells us that the normalization A prim,e − are around 0.5 and the spectral indices α of different pulsars are around 2.0, which is consistent with the radio and γ-ray observations of pulsar [55,56]. We can estimate the injection efficiency η from the pulsar. Take Geminga as an example, the spin-down energy loss rate of Geminga |Ė(T )| = 3.2 × 10 34 erg/s. The total radiation energy of the magnetic dipole can be derived from eq. (12) as E tot (T ) |Ė(T )|T 2 /τ 0 = 1.2 × 10 49 erg. From the fit, we get the injection energy E out /2 = 10 50.2 GeV 2.5 × 10 47 erg. From (13), we get η ∼ 4.2%. This efficiency is consistent with the previous studies by [24] and [33]. We can perform similar studies on the other five pulsars, whose results are listed in Table 2. A smaller η means it is easier for this pulsar to produce the same amount of positrons and electrons. The

Multiple pulsars interpretation
The extra high energy positrons may come from serveral pulsars. We perform similar study for multiple pulsars as we do for a single pulsar. Benefiting from the study in Section 4, we can assume that the spectral indices α and the electron injection efficiencies η for all the pulsars are the same to reduce the number of free parameters. The discussion on η from single pulsar in Section 4 tells us that Geminga, J1741-2054 and Monogem will give a much larger contribution to the high energy positron than J0942-5552. In other words, the η of J0942-5552 in Table 2 is much larger than that of Geminga, which implies that the contribution from J0942-5552 in the multiple pulsars interpretation can be negligible compared with that from Geminga. We choose three from the four "surviving" pulsars in the multiple pulsars discussion. The input parameters are the age T , the distance d and the energy loss rateĖ of each pulsar while the parameters we get from the fit is the normalization factor of primary electron A prim,e − , the spectral index α and the electron injection efficiency η. As shown in Fig. 4 (a), we obtain a good result from the multiple pulsar fit where χ 2 /ndf = 0.79. The parameters we get are A prim,e − = 0.48, α = 1.91 and η = 2.0%. The multiple pulsars interpretaion predicts a positron fraction with a decrease up to 600 GeV and after that a bump up to 2000 GeV, which is possible to be observed with more accumulating AMS-02 data.
Using the parameters from the fit, we can reproduce the electron flux measured by AMS-02 [2] in Fig. 4 (b). It shows that our electron background estimation in Section 2 matches the experimental data especially at high energy. For low energy, we simply take the solar modulation potential as 550 MV here. To reproduce the low energy electron flux more accurately, we need a monthly low energy electron fluxes, which may be published by AMS collaboration to model solar modulation.

Discussion and conclusion
In this work, we investigate the probability that the rise of the positron fraction measured by AMS-02 can be explained by pulsars. The propagation parameters and the injection spectrums of nuclei and electrons are tuned according to the Boron-to-Carbon ratio and the proton flux. It will be better to tune those parameters with Boron-to-Carbon ratio and proton flux measured by AMS-02 since they are in the same data taking period as the lepton fluxes. We find both the single pulsar model and the multiple pulsar model can explain the AMS-02 data very well. Six nearby pulsars are investigated as the single pulsar sources of the high energy positrons and finally four survive from all the conditions. With three mostly contributing pulsars, the multiple pulsars model predicts a positron fraction with a decrease up to 600 GeV and a bump up to 2000 GeV. For the low energy, a simple solar modulation potential potential can not explain the measurement well. Thus, we need the monthly electron fluxes which can describe solar activity during the whole period.
It is shown that the positron excess measured by AMS-02 can be explained by the pulsar scenario. Since the multiple pulsars can explain the experimental data well, it will be difficult to exclude pulsar scenario by isotropy. With accumulating AMS-02 data and future experiments, we can see the positron fraction behavior up to higher energy which will either confirm or reject the multiple pulsars scenario. If we consider other scenarios such as Dark Matter, we have to look into other productions, antiproton for instance, which have no contribution from pulsars.

A Solving the diffusion energy-loss equation without spherically symmetric approximation
Many previous studies (for example, Refs. [21,26,29,30,33]) use the spherically symmetric approximation to solve the diffusion energy-loss equation for a burst-like source. We point out that this approximation is not necessary. In this appendix, we show that the solution can be obtained without using the spherically symmetric approximation. The diffusion energy-loss equation is given by where f ( x, t, E) is the particle number density per unit energy interval, D(E) > 0 is the diffusion coefficient, b(E) ≡ −dE/dt > 0 is the energy loss rate, and Q( x, t, E) is the source term.

A.1 Green function for the diffusion energy-loss equation
The Green function G( x, t, E; x 0 , t 0 , E 0 ) of (18) is defined as The solution of (19) has been given in Ref. [57]. To fix the notation, let us briefly review the derivation. Define Substituting G = φ/b into (19) gives Let us make the variable transformation (t, E) → (t , λ) as follows, The Jacobian matrix of this transformation is easily obtained as whose inverse is Thus, It follows from eq. (24) that det which implies where in the second equality we have used the properties: D(E 0 ) > 0 and b(E 0 ) > 0. Substituting (29) into (27), we obtain As is well known, the Green function G( x − x 0 , λ) of the diffusion equation with λ ≥ 0. The solution of (31) is Comparing (30) with (31), we can read off the solution of φ as follows where in the second equality we have used (22). Thus, we finally get the solution of (19) as follows where the functions τ (E, E 0 ) and λ(E, E 0 ) are defined in eqs. (22) and (23), respectively.

A.2 Solution for a burst-like source
Once we know the Green function, eq. (34), we can write down the solution of eq. (18) for a generic source as follows In particular, for a burst-like source, the source function is proportional to where Q(E) is an arbitrary function of E, x 0 is the position of the source, and t 0 is the instantaneous time when the source bursts. Substituting eqs. (34) and (36) into eq. (35), we obtain Denote the solution of the equation is E = E 0 , then we have where we have used Substituting eq. (39) into eq. (37), we have where the initial energy E 0 is defined as the solution of eq. (38). In other words, if we know t − t 0 and E, we can find E 0 by solving the equation Define the diffusion distance r dif as we can rewrite eq. (41) as with r 2 ≡ ( x − x 0 ) 2 .

A.3 Solution for a burst-like source with power-law spectrum
Consider the case when the function Q(E) in eq. (36) is a power-law function, that is, where Q 0 is a normalization constant, α is the spectral index, and E cut is the cutoff energy. We assume that the diffusion coefficient D(E) and the energy loss rate b(E) take the form where β ≡ v/c is the ratio of velocity to speed of light, the constant D 0 and the index δ can be figured out by the background fitting in Sec. 2, the constant b 0 is given in Sec. 3, and E 0 should be determined by eq. (42). Let us calculate the initial energy E 0 first. Denote the diffusion time t dif ≡ t − t 0 . It follows from eqs. (42) and (47) that which implies from which, we see 1 which gives Thus, we obtain the following pieces in eq. (44): Substituting eqs. (52) and (53) into eq. (44), we obtain which is consistent with previous works (for example, eq. (6) of Ref. [33]). Now let us figure out the diffusion distance r dif . To this end, we substitute eqs. (46) and (47) into eq. (23), and get which can also be written as Comparing the above equation with eq. (46), E max ≡ 1/(b 0 t dif ) and E 0 /E = (1 − E/E max ) −1 , we obtain Thus, the diffusion distance r dif is given by which is consistent with previous works (for example, eq. (10) of Ref. [24] and eq. (7) of Ref. [33]).