TeV cosmic ray electrons from millisecond pulsars

Recent {\gamma}-ray observations suggest that the {\gamma}-ray millisecond pulsar (MSP) population is separated into two sub-classes with respect to the pair multiplicity. Here, we calculate the cosmic ray electron/positron spectra from MSPs. Based on the assumption of the equipartition in the pulsar wind region the typical energy of electrons/positrons ejected by a MSP with the pair multiplicity of order unity is \sim 50 TeV. In this case, we find that a large peak at 10 - 50 TeV energy range would be observed in the cosmic ray electron/positron spectrum. Even if the fraction of pair starved MSPs is 10%, the large peak would be detectable in the future observations. We also calculate the contribution from MSPs with high pair multiplicity to the electron/positron spectrum. We suggest that if the multiplicity of dominant MSP population is \sim 10^3, electrons/positrons from them may contribute to the observed excess from the background electron/positron flux and positron fraction.


INTRODUCTION
The Fermi Gamma-Ray Space Telescope has detected γ-ray pulsed emissions from more than twenty millisecond pulsars (MSPs) (Abdo et al. 2011), which have a rotation angular frequency Ω ∼ 10 3 s −1 and a stellar surface magnetic field Bs ∼ 10 8.5 G. The detection of the GeV emissions from a pulsar magnetosphere means that electrons and positrons are accelerated to more than ∼ TeV by the electric field parallel to the magnetic field, which arises in a depleted region of the Goldreich-Julian (GJ) charge density (Goldreich & Julian 1969). The γ-ray light curve is an important tool for probing the particle acceleration process in the pulsar magnetosphere. Therefore, the γ-ray emission region has been explored by comparing theoretical models such as polar cap (Daugherty & Harding 1996), outer gap (Cheng, Ho & Ruderman 1986) and slot gap models (Muslimov & Harding 2004a) with the observed light curve (e.g., Venter, Harding & Guillemot 2009;Romani & Watters 2010;Kisaka & Kojima 2011). Venter, Harding & Guillemot (2009) fitted the pulse profiles of the Fermi detected MSPs with the geometries of γ-ray emission region predicted by different theoretical models. They found that the pulse profiles of six of eight MSPs could be fitted by the geometries of either the outer gap or the slot gap model, as was the case of canonical pulsars. They interpreted that copious pairs are produced in the magnetosphere of these MSPs. However, Venter, Harding & Guillemot (2009) also found that the pulse profiles of remaining two MSPs show the unusual behavior in the γ-ray light curves and could not be fitted by the geometry of either the outer gap or the slot gap models. They proposed that these unusual light curves could be fitted by the pair starved polar cap model (Muslimov & Harding 2004b), in which the multiplicity of the pairs is not high enough to completely screen the electric field above the polar cap, and the particles are continuously accelerated up to high altitude over the entire open field line region. Thus, from the model fitting of the γ-ray light curves, Venter, Harding & Guillemot (2009) suggested that the γ-ray MSP population is separated into two sub-classes.
The important fact is that radio pulsed emission is also detected from all currently detected γ-ray MSPs and remarkably similar to that from canonical pulsars. The pulsar radio emission is a highly coherent process because the brightness temperature is extremely high. In the theoretical models of the radio emission mechanisms, some authors have believed the conditions that there are a highly relativistic primary beam with the large Lorentz factor (∼ 10 7 ) and the number density nearly equal to the GJ density, and the secondary electron/positron plasma with relatively small bulk streaming Lorentz factor (∼ 10 -10 3 ) and the large pair multiplicity (∼ 10 3 -10 5 ) in the radio emitting region (e.g., Melrose 1995;Lyutikov, Blandford & Machabeli 1999;Gedalin, Gruman & Melrose 2002). However, the existence of pair starved MSPs suggests that the radio emission mechanisms should be insensitive to the particle number density down to sub-GJ number density. The pulsar radio emission mechanism is still poorly understood, so that the observationally-based constraints are valuable (Melrose 1995). Therefore, another verification for the extent of the MSP multiplicity, especially the existence of the pair starved MSPs is important for the pulsar radio emission mechanisms.
Recently, HESS has discovered a new TeV source, which is located in the close vicinity of the globular cluster Terzan 5 (Abramowski et al. 2011). Several globular clusters, including Terzan 5 also emit GeV γ-ray (Abdo et al. 2009;Kong, Hui & Cheng 2010;Abdo et al. 2010;Tam et al. 2011), which may plausibly be due to a number of MSPs residing in these clusters (Harding, Usov & Muslimov 2005;Venter, de Jager & Clapson 2009). Thus, inverse Compton scattering by the high-energy particles ejected from MSPs are proposed for the origin of the observed TeV emission (Bednarek & Sitarek 2007;Venter, de Jager & Clapson 2009). The high-energy electron/positron spectrum ejected from MSPs would be a useful probe for the multiplicity of the MSPs. However, only from the TeV spectra, we cannot distinguish two models (Bednarek & Sitarek 2007;Venter, de Jager & Clapson 2009), which assume different pair multiplicities (Abramowski et al. 2011). Another way to investigate the electron/positron spectrum ejected from MSPs is its direct measurement. Since high-energy electrons/positrons can propagate only about a few kpc due to the energy losses by the synchrotron and the inverse Compton emission, the direct detection of the electrons/positrons ejected from MSPs in the globular clusters is unlikely. However, for the following reasons, we may detect those from nearby MSPs.
MSPs have much lower spin-down luminosity than canonical pulsars.  investigated the possible contribution of the nearby MSP, PSR J0437-4715 to the cosmic ray electron/positron spectrum. They concluded that unlike canonical pulsars such as Geminga pulsar, the contribution from a MSP to the observed electron/positron flux is negligible. However, since the lifetime of MSPs is much longer than canonical pulsars (> 10 10 yr), there should be much more nearby active MSPs. Furthermore, hereafter KIK11) pointed out that since white dwarf pulsars have long lifetime and continue to inject the electrons/positrons after the nebulae stop expanding, the adiabatic energy losses of electrons/positrons in the pulsar wind nebula region are negligible. Also the synchrotron cooling of electrons/positrons is so small and the high-energy electrons/positrons can escape the nebulae without losing much energy. Although they consider the case of the white dwarfs, their results are also applicable to MSPs. Therefore MSPs could potentially contribute to the observed high energy cosmic ray electrons/positrons and will be detectable by the next generation experiments, such as CALET (Torii et al. 2008a) and CTA (CTA Consortium 2010).
In this paper, we investigate the contribution of electrons/positrons ejected from the MSPs to the observed cosmic ray spectrum. In section 2, we apply KIK11 model to the case of MSPs. We estimate the typical energy of electrons/positrons from the MSPs and show that during the propagation in pulsar wind nebulae, the adiabatic losses and radiative cooling of electrons/positrons are not so large. We also describe the propagation in the interstellar medium (ISM). In section 3, we calculate the energy spectrum of cosmic ray electrons/positrons from the MSPs and show the possibility that the electrons/positrons from these MSPs are detectable for the future observations.

Acceleration and cooling
In order to estimate the energy of electrons/positrons available in the wind region and their adiabatic and radiative cooling in the shocked region, we adopt the model of KIK11. For a pulsar wind nebula formed by a MSP, we consider the conditions that the relativistic wind blasts off from the light cylinder ∼ R lc = c/Ω where c is the speed of light, and two shock fronts are formed between the supersonic pulsar wind and ISM. Since the energy from a MSP is continuously transported to the wind, the shock fronts are expanding until the pressure of the shocked region P sh becomes equal to that of ISM p. Although KIK11 considered the case of white dwarf pulsars, the situations are similar to the case of MSPs because they have a long lifetime and the supernova shock front have already decayed. We assume that the effects of binary companion are negligible, because the fraction of solid angle occupied by companion is small (<1%). We also neglect radiative loss due to curvature radiation within light cylinder (∼ 10%). From now on we set fiducial parameters of the MSP's surface magnetic field strength, angular frequency and radius as B0 = 10 8.5 G, Ω = 10 3 s −1 and R = 10 6 cm, respectively.
We assume the energy equipartition between particles and magnetic field, εeN = B 2 /8π, and the conservation of the particle number flux, 4πr 2 cN ∼ constant, in the MSP wind region. Here, N is the number density of electrons/positrons. The number density can be described as where κ is the multiplicity of electrons/positrons, N lc and B lc are the number density and the magnetic field at the light cylinder, respectively. We assume the magnetic field configuration as pure dipole (B ∝ r −3 ) within light cylinder. Outside the light cylinder, we assume the conservation of the energy flux of the magnetic field, Br ∼ constant. Using these assumptions, the typical energy of electrons/positrons εe can be described as where ψmax is the electric potential difference across the open magnetic field lines described as ( These values are similar to those in the case of white dwarf pulsars (KIK11). Note that the typical energy depends on the pair multiplicity. Next, we estimate the adiabatic and the radiative cooling of electrons/positrons in the shocked region. The outer shock of the pulsar wind nebula finally decays when the pressure of the shocked region P sh becomes equal to that of the ISM p. If the outer shock decaying time is shorter than the lifetime of MSP, the adiabatic cooling is negligible. In order to estimate the outer shock decaying timescale, we solve the equation of motion and the energy conservation law at the outer shock front. The equation of motion can be described as where Rout is the radius of the outer shock front, P sh is the pressure of the shocked region and ρ is the density of ISM. The energy equation is d dt where L sd is the spin-down luminosity of MSP, In the derivation of eq.(5), we assume that in the shocked region the internal energy of particles is 3P sh /2 because the energy of particles is the relativistic regime. Using the typical value for the density of ISM ρ ∼ 10 −24 g cm −3 , the pressure in ISM is p ∼ 10 −13 dyn cm −2 . Solving above equations, the outer shock decays at about where T is the temperature of ISM. The lifetime of a MSP τ can be estimated as where Erot is the rotation energy described as where M is the mass of a MSP. For the fiducial parameters of MSP We find that the outer shock decays at a very early stage of the lifetime of MSP and does not expand after t > t dec . Therefore, in the similar to the case of white dwarf pulsars (KIK11), the adiabatic cooling shall give minor contributions to the cooling process of the high-energy electrons/positrons. Also for the radiative cooling by the synchrotron radiation in the shocked region Rin < r < Rout, we follow the discussion in KIK11 based on the diffusion in the shocked region. We take the Bohm limit, where the fluctuation of the magnetic field is comparable to the coherent magnetic field strength. In this limit, the timescale t diff for the electron/positron trapping in the shocked region is given by where D sh = crg/3 is the diffusion coefficient under the Bohm limit, d is the size of the shocked region and rg = εe/eB is the Larmor radius of electron/positron with energy εe. To estimate the diffusion timescale, we need to know the size and the magnetic field strength of the shocked region. Here we consider the time t > t dec , so that the size of the shocked region is an order of the radius of forward shock front at t = t dec , d ∼ Rout(t = t dec ) ∼ 3 × 10 19 cm. For the radius of the inner shock front at t > t dec , we use the balance condition between the momentum transferred by wind and the pressure of ISM For the fiducial parameters, Rin(t > t dec ) ∼ 2 × 10 17 cm. The strength of the magnetic field at the inner radius Bin can be estimated as Bin ∼ 2 × 10 −6 G, which is almost the same as that of ISM. Then, the diffusion timescale is The synchrotron energy loss of a particle with energy εe is described as where σT is the Thomson scattering cross section, β is the particle velocity normalized by the speed of light and me is the mass of electron/positron. Then, The typical energy loss of the electrons/positrons ∆εe with energy εe can be estimated as This means that the high-energy electrons/positrons injected into the shocked region lose roughly 30% of the energy by the synchrotron radiation before diffusing out into ISM. Therefore, as in the case of white dwarf pulsars (KIK11), we can conclude that the radiative energy loss of electrons/positrons in the pulsar wind nebula is not so large. The above expressions for the estimate of the energy losses are only applicable to the case that the velocity of a MSP is subsonic in ISM. The observed velocity of MSPs is less than that of canonical pulsars in average sense . However, some MSPs have the large velocity and a few MSPs actually forms bow shock nebulae (Stappers et al. 2003;Hui & Becker 2006). In this case, the size of the bow shock is described as (e.g., Wilkin 1996) where V is the velocity of a MSP. Due to the assumption of the energy equipartition, the strength of the magnetic field is The ratio of the Larmor radius of electrons/positrons to the bow shock radius is rg R bow ∼ 0.5κ −1 .
The fact that rg/R bow is close to unity supports that electrons/positrons may escape from the bow shock region. Therefore, in the case of κ ∼ 1, high-energy electrons/positrons can escape with an efficiency of order unity (Bandiera 2008;Blasi & Amato 2010). Even if we consider the case of κ ≫ 1, the synchrotron loss can be estimated by using eqs. (11), (14) and (16) Therefore, we can conclude again that the radiative energy loss of electrons/positrons in the pulsar wind nebula is not so large.

Diffusion in Interstellar medium
The observed electron/positron spectrum after the propagation in ISM is obtained by solving the diffusion equation where f (t, r, εe) is the energy distribution function of electrons/positrons, D(εe) = D0(1+εe/3GeV) δ is the diffusion coefficient, P (εe) is the cooling function of the electrons/positrons which takes into account synchrotron emissions and inverse Compton scatterings during the propagation, and Q(t, εe, r) is the injection term. Here we adopt D0 = 5.8 × 10 28 cm 2 s −1 , δ = 1/3, which is consistent with the boron-to-carbon ratio according to the latest GALPROP code. Atoyan et al. (1995) showed a solution in the case of an instantaneous injection from a single point-like source, i.e. Q(t, εe, r) ≈ Q0(εe)δ(t − ti)δ(r). Then the observed spectrum G(t, r, εe;t) would be G(t, r, εe;t) = Q0(εe,0)P (εe,0) π 3/2 P (εe)d diff (εe, εe,0) 3 exp − where εe,0 is the energy of electrons/positrons at the timet(< t) and which are cooled down to εe at the time t, and d diff is the diffusion length given by The cooling function P (εe) can be described as where utot(εγ )dεγ is the energy density of interstellar radiation fields with the photon energy between εγ and εγ + dεγ (including cosmic microwave background, starlight, and dust emission; Porter et al. 2008), and B is the interstellar magnetic field with we here set as 1µG. Here the function fKN(x) is the correction factor to include the Klein-Nishina (KN) effect, which approaches to unity when x is much smaller than unity. The mathematical expression of fKN can be found in Moderski et al. (2005). As shown in the last section, in general MSPs have a very long lifetime (τ ∼ 5 × 10 10 yr) which is comparable with the cosmic age. In such a case that a point-like source with a finite duration, taking into account the time-dependence of an injection rate (Q0(εe) → Q0(εe,t)), the spectrum can be calculated by integrating G(t, r, εe; τ ) for τ : where ti is the time when the particle injection from a source has started. Here we assume that the electron/positron injection spectrum can be described as mono-energetic distribution or power-law distribution where α is the intrinsic power-law index of an electron/positron spectrum and εcut is the maximum electron/positron energy from a source. Now we can calculate the average electron/positron spectrum by considering the birth rate of MSPs as follows: where t0 is the cosmic age (i.e. the present time), and R is the local pulsar birth rate (yr −1 kpc −2 ). Here εe,i is the energy at the time ti of CR electrons/positrons which are cooled down to εe at t.

RESULTS AND DISCUSSIONS
First, we calculate the cosmic ray electron/positron spectra from the pair starved MSPs. We set the pair multiplicity κ = 1, the lifetime τ = 5 × 10 10 yr, the total energy Erot = 10 52 erg, the local birth rate R = 3 × 10 −9 yr −1 kpc −2 and the fraction of the lost energy due to synchrotron emission 30% for each MSP. We assume that each MSP has the same value of the parameters (B0 = 10 8.5 G, Ω = 10 3 s −1 , R = 10 6 cm), because most MSPs have the almost same spin-down luminosity. For the injection distribution function, we assume mono-energetic distribution eq.(25) with the energy εe = 50 TeV. Even if we consider the power-law distribution eq.(26), the energy range of the distribution is small because the allowed range of the cutoff energy should be only 50 -80 TeV due to the observed constraint by KASKADE/GRAPES/CASA-MIA (Kistler & Yüksel 2009). Thus the distribution should be nearly mono-energetic distribution. Note that our local MSP birth rate is based on the MSP local surface density of 38± 16 pulsars kpc −2 for 430 MHz luminosity above 1 mJy kpc 2 (Lorimer 2008). Actually, now 23 MSPs are detected within 1kpc (ATNF catalog, Manchester et al. 2005). We can only detect MSPs that have the radio beam directed toward us and the radio flux larger than the threshold of detectors. However, the cosmic ray electrons/positrons ejected from MSPs are distributed isotropically because of the effect of Galactic magnetic field, so that a large number of MSPs will contribute to the observed electron/positron spectrum. Therefore, our local MSP birth rate corresponds to the lower limit for the current radio observations. In figure 1, the electron/positron flux from multiple pair starved MSPs is shown. Thick solid, dashed and dash-dotted lines show the total electron/positron spectra if the fraction of the pair starved MSPs is 100%, 25% and 10%, respectively. Thin solid line shows the contribution of electrons/positrons from the pair starved MSPs when the fraction is 100%. The background flux is shown as a thin dashed line. We adopt the background model of an exponentially cutoff power law with an index of -3.0 and a cutoff at 1.5 TeV, which is similar to that shown in Aharonian et al. (2008) and reproduces the data in ∼10 GeV-1 TeV well. It is very interesting that there is a large peak at 10-50 TeV energy range. The existence of this peak cannot be ruled out from the current observations. The high-energy component is more enhanced for the longduration sources (Atoyan, Aharonian & Völk 1995;Kawanaka, Ioka & Nojiri 2010). This is because the longer the duration of injection is, the larger fraction of fresh electrons/positrons are expected to reach the Earth without losing their energy so much during the propagation. MSPs are the continuously injecting sources with the duration as long as the cosmic age, so that the spectrum from them has nearly the same shape with the injection spectrum with the soft energy tail component. This is the difference from other sources such as young pulsars, whose typical duration is only ∼ 10 4 -10 5 yrs. Fitting result of Venter, Harding & Guillemot (2009) showed that the fraction of the pair starved MSPs is 25%, although they have only eight samples. Even if the fraction is 10%, the flux is ∼ 20m −2 s −1 sr −1 GeV 2 at 10 TeV. In this case, we can detect the electron/positron flux with near future missions such as CALET (we assume the geometrical factor times the observation time ∼ 5 yrs as ∼ 220 m 2 sr days) because the predicted electron/positron flux is sufficiently large. It was considered that the number of astrophysical sources contributing to the above several TeV energy range is quite small according to the birth rate of the supernovae and the canonical pulsars in the vicinity of the Earth (Kobayashi et al. 2004;Kawanaka, Ioka & Nojiri 2010;Kawanaka et al. 2011). However, we find that it is possible for multiple pair starved MSPs to contribute to the 10 TeV energy range in the electron/positron spectrum. Therefore, if the anisotropy of the observed electrons/positrons are weak in the 10 TeV energy range, we suggest that the pair starved MSPs may contribute to the spectrum significantly.
Next, we also investigate the contribution of MSPs with high pair multiplicity to the observed cosmic ray electrons/ positrons. In this model, we assume that the injection function of these MSPs is power-law distribution with index α = 2, the cutoff energy εe,cut = 1TeV and the minimum energy εe,min = 1GeV. In this case, the pair multiplicity is ∼ 2000. Other parameters are the same values as in the case of the pair starved MSPs. We assume that the fraction of MSPs with high multiplicity is 100%. The results are shown in figure 2 for the electron/positron spectrum and in figure 3 for the positron fraction. Both figures show that electrons/positrons ejected from MSPs with high multiplicity partially contribute to the excess from background flux observed by PAMELA, HESS and Fermi. In this energy range, the other sources such as canonical pulsars (Shen 1970;Atoyan, Aharonian & Völk 1995;Chi, Cheng & Young 1996;Zhang & Cheng 2001;Kobayashi et al. 2004;Grimani 2007;Profumo 2012;Hooper, Blasi & Serpico 2009;Yüksel, Kistler & Stanev 2009;Malyshev, Cholis & Gelfand 2009;Grasso et al. 2009;Kawanaka, Ioka & Nojiri 2010;Heyl, Gill & Hernquist 2010) would also contribute to the observed electron/positron spectrum. Note that even if other sources are dominant for the observed excess, the total spectrum added to the contribution of MSPs with high multiplicity does not significantly exceed the observed electron/positron spectrum.

SUMMARY
In this paper, we show the possibility that the cosmic ray electrons/positrons from MSPs would significantly contribute to the observed spectrum. Although MSPs have relatively low spin-down luminosity and low birth rate, the lifetime is so long that there are many active MSPs in the vicinity of the Earth. Furthermore, such a long lifetime source continuously injects electrons/positrons after the nebula ceases expanding, so the adiabatic energy losses in a pulsar wind nebula region are negligible. The synchrotron cooling in the nebula is also small, so the high-energy electrons/positrons can escape the nebula without losing much energy. We calculate the diffusive propagation of high-energy electrons/positrons in the ISM taking into account the cooling via synchrotron emissions and inverse Compton scatterings, and predict their spectrum observed at the Earth.
In the case of the MSPs with multiplicity that is unity, the typical energy of the electrons/positrons produced should be ∼ 50 TeV based on the assumption of the equipartition in the MSP wind region. Since the long duration of injection make a hard spectrum, the peak is enhanced at 10 -50 TeV energy range. Even if the fraction of pair starved MSPs is as small as 10%, this peak would be detectable in the future missions such as CALET and CTA. Although a single young source can make the similar spectral feature in this energy range, in the case of pair starved MSPs the anisotropy of the electron/positron flux would be weaker because a number of sources contribute to it. If this peak is detected, that will be a great impact for on the studies of pulsar radio emission mechanisms because the existence of pair starved MSPs suggests that the radio emission mechanisms should be insensitive to the particle number density down to sub-GJ number density. The detection also suggests that the current outer gap model should be modified because Wang & Hirotani (2011) suggested that most MSPs locate above the pair death line of the outer gap and the multiplicity is larger than unity.
We also calculate the electron/positron spectrum from MSPs with high pair multiplicity. We suggest that if multiplicity of these MSPs is the order of ∼ 10 3 , electrons/positrons from them partially contribute to the observed excess of the total spectrum and the positron fraction.  Kistler & Yüksel 2009). We also show the total spectra from pair starved MSPs with the different fractions: 25% (thick dashed line) and 10% (dot-dashed line). We assume that the lifetime τ = 5 × 10 10 yr, the total energy Erot = 10 52 erg, the local birth rate R = 3 × 10 −9 yr −1 kpc −2 and the fraction of the energy loss is 30%.  Figure 2, and the background (dotted line), compared with the PAMELA data as red circles (Adriani et al. 2009). Note that the solar modulation is important below ∼ 10 GeV.