# Modeling of the Near-Earth Low-Energy Antiproton Fluxes

Authors: U. B. Jayanthi; A. A. Gusev
Journal: Advances in Astronomy
 (2011)

---

## Abstract

The local interstellar antiproton spectrum is simulated taking into account antineutron decay, (He,p) interaction, secondary and tertiary antiproton production, and the solar modulation in the “force field” approximation. Inclusive invariant cross-sections were obtained through a Monte Carlo procedure using the Multistage Dynamical Model code simulating various processes of the particle production. The results of the simulations provided flux values of4·10−3 to 10−2
and 10−2 to 1.7·10−2 antiprotons/(m2 s sr GeV) at energies of 0.2 and 1 GeV, respectively, for the solar maximum and minimum epochs. Simulated flux of the trapped antiprotons in the inner magnetosphere due to galactic cosmic ray (GCR) interactions with the atmospheric constituents exceeds the galactic antiproton flux up to several orders. These simulation results considering the assumptions with the attendant limitations are in comprehensive agreement with the experimental data including the PAMELA ones.


---

## Body

## 1. Introduction

The interest of antimatter component in the cosmic radiation ranges from the basics of cosmology [1] to a possible utility of antimatter fuel in future interplanetary missions [2]. The antiprotons are perhaps the most intensively studied constituent of the antimatter accessible for observation. The experimental results of their measurements from [3] to [4] agree well with the hypothesis of their secondary origin [5]. However, it does not exclude primary constituents produced, for example, in evaporation of primordial black holes or dark matter annihilations.Before considering the existence of such exotic sources, it is imperative to evaluate precisely the antiproton production cross sections, uncertainties due to the propagation in the galaxy and heliosphere, and additional sources of antiproton production in the Earth’s environment. The recent experimental data from different experiments at different solar activity phases provided a good opportunity for understanding the modulation process as well as the similarity in the proton and antiproton transport in the galaxy, interplanetary space, and the Earth’s magnetosphere in spite of fluctuations in data and the model approximations.


## 2. Interstellar Antiproton Fluxes

In our simulation the antiproton local interstellar spectrum (LIS)Fp̃(Ep̃)was obtained through the leaky-box model as a solution of the following system of integro-differential equations [6]:(1a)Fp̃(Ep̃)λesc+Fp̃(Fp̃)λinel+ddEp̃(Fp̃(Ep̃)〈dEdx〉(Ep̃))=Q2p̃(Ep̃)+Q3p̃(Ep̃),(1b)Q2p̃(Ep̃)=NAvog∑j=H,He,onjAj∫Ep,th∞dσjdEp̃(Ep̃,Ep̃,ñ)FP(EP)dEP,(1c)Q3p̃(Ep̃)=NAvog∑j=H,He,onjAj∫0∞dσjdEp̃(Ep̃,E′p̃)Fp̃(E′p̃)dE′p̃. This considers the production of the secondary antiprotons Q2p̃        by galactic cosmic ray (GCR) proton flux Fp(Ep) and subsequent tertiary antiprotons Q3p̃ energy losses 〈dE/dx〉 in the interstellar matter, the flux decreases due to escape (λesc) and inelastic interaction (λinel).


### 2.1. Antiproton Production Spectrum

The antiproton production spectrum, that is, the source functionQ2p̃ + Q3p̃        in equations (1a)–(1c) is a sum of the contributions from interactions of the GCR protons, He nuclei, and antiprotons with the interstellar H, He, and O nuclei in the interstellar matter. The corresponding densities nj are 1, 0.1, 8·10-4cm−3 [7]. The antineutron production was also taken into account because they totally decay into antiprotons during confinement in the galaxy. The decaying antineutrons pass practically all their energy to the antiprotons with identical energy spectra.The production of antiprotons and antineutrons was simulated with a Multistage Dynamical Model (MSDM) Monte Carlo code [8]. The code produces energy spectra and angular distributions of the reaction products (x,A) of an incident particle x with a target nucleus A together with total and inelastic cross sections and multiplicities. The projectile can be a hadron (n,ñ,p,p̃) or a meson (π+,π-,π0,K+,K-,K0) with kinetic energies from 10 MeV up to 1 TeV. The target can be any nucleus with the atomic mass A≥1. All the models included in the code were comprehensively tested by adjusting to the experimental data.The MSDM code simulates all the stages of hadron-nucleus and nucleus-nucleus interactions inside the target using an exclusive approach based on the models described in [9]. It considers the cascade and precompound stages of the reaction as well as evaporation/fission, multifragmentation, and Fermi breakup of residual nuclei. At the cascade stage for projectile energies below 1 GeV, when consideration of only nucleons, pions, and D-resonances is sufficient, an original intranuclear cascade model [10] is used fitted to dynamics of formation and absorption of pions [11]. For energies exceeding ≈10 GeV the independent quark-gluon string model is applied.Hadron-hadron interactions for energiesEh<1GeV are simulated on the base of parameterization of experimental data and πΔ-dynamics, for energies 1GeV≤Eh≤10GeV the improved version of the quark gluon string (QGS) model is used, and for Eh>10 GeV the model of independent QGS is used.Figure1 shows differential cross sections of the antiprotons dσH(Ep,Ep̃)/dEp and antineutrons dσH(Ep,Eñ)/dEp for proton projectiles of various energies Ep in (p,H) reaction (scattering of the data in Figures 1–4 is determined by only the statistic accumulated in Monte Carlo simulation. The overall statistics of several hundred millions events provides a necessary convergence of the final result).Figure 1
Antiproton (aP) and antineutron (aN) differential cross sections for(p,H) reaction. Symbols mark the MSDM results; thick solid lines represent the Tan and Ng [12] approximation.Figure 2
Inelastic, total(p,H) cross sections, and multiplicities simulated with the MSDM code compared with the approximation from [12].Figure 3
Antiproton and antineutron differentia cross-section ofp̃+H reaction simulated with the MSDM code.Figure 4
The interstellar secondary antiproton production spectra simulated using the MSDM and [12] cross sections.One can see that the MSDM code provides some excess of the antiprotons over the approximation [12] widely used for simulations in GCR energy range. The spectra produced in (p,He) and (p,O) reactions are identical. The excess of the antineutrons over the antiprotons clearly seen in the figure is explained by a specific mechanism of their generation through the evaporation process favoring neutral particle production in the absence of the electromagnetic interaction. Inelastic, total (p,H) cross sections, and multiplicities are shown in Figure 2.Figure3 presents examples of spectra of the antiprotonsdσH(Ep̃,Ep̃)/dEp̃, and antineutrons σH(Ep̃,Ep̃)/dEp̃, produced by antiproton projectiles of various energies Ep̃ in (p̃,H) reactions. The spectra noticeably differ from the uniform distributions ∝1/projectile energy from [12]. In contrast with (p,H) reaction, which produces practically identical spectra of p̃ and ñ at all the projectile energies those spectra from (p̃,H) reactions are similar only at higher projectile energies and are strongly different at low ones.As a parent LIS spectrum for the simulation a power low approximationFp(Ep)=1.37104Ep-2.732· (m2sr s GeV)−1 for protons and FHe(EHe)=7.06103EHe-2.699· (m2sr s GeV)−1 for He [13] are used. The production spectra Q2p̃ obtained from Fp(Ep) are depicted in Figure 4. The cross section invariance for interaction of the GCR He with interstellar H was treated as a convolution of He LIS spectrum with dσHe(Ep,Ep̃)/dEp, where for He the Ep is expressed as energy per nucleon.


### 2.2. Antiproton LIS

In the simulation, the rigidity(R) dependent escape path length for antiprotons in the galaxy [14] is used:(2)λesc=11.8βforR<4.9GV,λesc=11.8β(R/4.9)0.54forR≥4.9GV,
the interaction length λinel of antiprotons including annihilation is also simulated with the MSDM and the stopping power 〈dE/dx〉 is calculated utilizing standard procedure (e.g., [15]).Equations (1a), (1b), and (1c) is solved through an iteration procedure using the “Mathematica” package. The solution readily converges in the third iteration. As shown in Figure 4 the MSDM cross section provides about two-times larger tertiary output Q3p̃  in the range of 0.3–3 GeV as compared to those from the uniform distribution used in [12]. In the energy range of 0.04–2 GeV the LIS obtained with MSDM slightly exceeds that obtained with the approximation [12] with the maximum deviation of ≤40% atEp̃=0.2 Gev.


### 2.3. Solar Modulation of the LIS

The LIS modulation in the heliosphere is considered on the basis of a transport equation for the spherically symmetric case [16]. The “force field” approximation inherently neglects the heliospheric gradient and curvature drifts but considers the diffusion, convection, and adiabatic deceleration:(3)F1AU(E1AU)P1AU2=FHB(EHB)PHB2,E1AU=Pc1nP1AU+E1AUPc+Ec+Ec+ΦforP<Pc,E1AU=EHB+ΦforP<Pc,Φ=V(RHB-1AU)3A,E1AU2=P1AU2+m02,EHB2=PHB2+m02.HereFHB, EHB, PHB, and F1AU, E1AU, and  P1AU  are the proton flux, proton total energy in GeV, momentum in GV at heliospheric boundary (HB) and at the Earth’s orbit (1 AU), respectively, m0 is the proton rest mass in GeV. V is the average solar wind speed in 103 km/hr. Physical sense of the solution implies a conservation of the distribution function F/P2 for particle energy decreases from EHB down to E1AU in travel from heliosphere RHB to the Earth at 1 AU.The heliospheric conditions are described by the “force field” parameterΦ determined by the solar wind speed V and the heliospheric boundary distance RHB. In our simulation we used A=17, Pc=1.015GV [17], where it was shown that for the energies ≥0.02 GeV  (2) approximates the exact solution of the equation in [16] for the power proton spectrum ∝βEHB-γ (where β=PHB/EHB  is the proton speed) and also consistent with the solar flare proton observations.TheΦ magnitude is determined from the best fit approximation with (2) of the observed proton spectrum F1AU assuming the interstellar spectrum as FHB(EHB)=16470βEHB-2.76protons/m s sr GeV. The fits results are Φmax=0.964 GeV and Φmin=0.368 GeV corresponding to solar maximum and minimum epochs.The results of our simulated spectrum are shown in Figure5 together with the antiproton fluxes obtained in different experiments performed at different solar minimum and maximum periods. The increases in the low-energy fluxes are provided by the higher fluxes of more energetic particles enriching the <1 GeV region due to adiabatic energy losses. The steeper the low-energy branch of the LIS spectrum the more pronounced the above-mentioned increases and [18]. The results of the simulations provided flux values of ≈4·10−3 to 10−2 and 10−2 to 1.7·10−2  antiprotons/m2 s sr GeV at energies of 0.2 and 1 GeV, respectively, corresponding to the solar maximum and minimum epochs. The curve for Φ=1.5 is the lower limit for all the experimental data. It may correspond, for example, to V=103 km/hour and RHB=70AU.Figure 5
Simulated LIS and the modulated spectra compared with experimental observations during the different phase of solar cycles: BESS95+97 [19]; BESS 00 [20], BESS 02 [21], BESS-Polar II [22], and PAMELA [4]. “Force field” parameter Φ in GeV.


## 2.1. Antiproton Production Spectrum

The antiproton production spectrum, that is, the source functionQ2p̃ + Q3p̃        in equations (1a)–(1c) is a sum of the contributions from interactions of the GCR protons, He nuclei, and antiprotons with the interstellar H, He, and O nuclei in the interstellar matter. The corresponding densities nj are 1, 0.1, 8·10-4cm−3 [7]. The antineutron production was also taken into account because they totally decay into antiprotons during confinement in the galaxy. The decaying antineutrons pass practically all their energy to the antiprotons with identical energy spectra.The production of antiprotons and antineutrons was simulated with a Multistage Dynamical Model (MSDM) Monte Carlo code [8]. The code produces energy spectra and angular distributions of the reaction products (x,A) of an incident particle x with a target nucleus A together with total and inelastic cross sections and multiplicities. The projectile can be a hadron (n,ñ,p,p̃) or a meson (π+,π-,π0,K+,K-,K0) with kinetic energies from 10 MeV up to 1 TeV. The target can be any nucleus with the atomic mass A≥1. All the models included in the code were comprehensively tested by adjusting to the experimental data.The MSDM code simulates all the stages of hadron-nucleus and nucleus-nucleus interactions inside the target using an exclusive approach based on the models described in [9]. It considers the cascade and precompound stages of the reaction as well as evaporation/fission, multifragmentation, and Fermi breakup of residual nuclei. At the cascade stage for projectile energies below 1 GeV, when consideration of only nucleons, pions, and D-resonances is sufficient, an original intranuclear cascade model [10] is used fitted to dynamics of formation and absorption of pions [11]. For energies exceeding ≈10 GeV the independent quark-gluon string model is applied.Hadron-hadron interactions for energiesEh<1GeV are simulated on the base of parameterization of experimental data and πΔ-dynamics, for energies 1GeV≤Eh≤10GeV the improved version of the quark gluon string (QGS) model is used, and for Eh>10 GeV the model of independent QGS is used.Figure1 shows differential cross sections of the antiprotons dσH(Ep,Ep̃)/dEp and antineutrons dσH(Ep,Eñ)/dEp for proton projectiles of various energies Ep in (p,H) reaction (scattering of the data in Figures 1–4 is determined by only the statistic accumulated in Monte Carlo simulation. The overall statistics of several hundred millions events provides a necessary convergence of the final result).Figure 1
Antiproton (aP) and antineutron (aN) differential cross sections for(p,H) reaction. Symbols mark the MSDM results; thick solid lines represent the Tan and Ng [12] approximation.Figure 2
Inelastic, total(p,H) cross sections, and multiplicities simulated with the MSDM code compared with the approximation from [12].Figure 3
Antiproton and antineutron differentia cross-section ofp̃+H reaction simulated with the MSDM code.Figure 4
The interstellar secondary antiproton production spectra simulated using the MSDM and [12] cross sections.One can see that the MSDM code provides some excess of the antiprotons over the approximation [12] widely used for simulations in GCR energy range. The spectra produced in (p,He) and (p,O) reactions are identical. The excess of the antineutrons over the antiprotons clearly seen in the figure is explained by a specific mechanism of their generation through the evaporation process favoring neutral particle production in the absence of the electromagnetic interaction. Inelastic, total (p,H) cross sections, and multiplicities are shown in Figure 2.Figure3 presents examples of spectra of the antiprotonsdσH(Ep̃,Ep̃)/dEp̃, and antineutrons σH(Ep̃,Ep̃)/dEp̃, produced by antiproton projectiles of various energies Ep̃ in (p̃,H) reactions. The spectra noticeably differ from the uniform distributions ∝1/projectile energy from [12]. In contrast with (p,H) reaction, which produces practically identical spectra of p̃ and ñ at all the projectile energies those spectra from (p̃,H) reactions are similar only at higher projectile energies and are strongly different at low ones.As a parent LIS spectrum for the simulation a power low approximationFp(Ep)=1.37104Ep-2.732· (m2sr s GeV)−1 for protons and FHe(EHe)=7.06103EHe-2.699· (m2sr s GeV)−1 for He [13] are used. The production spectra Q2p̃ obtained from Fp(Ep) are depicted in Figure 4. The cross section invariance for interaction of the GCR He with interstellar H was treated as a convolution of He LIS spectrum with dσHe(Ep,Ep̃)/dEp, where for He the Ep is expressed as energy per nucleon.


## 2.2. Antiproton LIS

In the simulation, the rigidity(R) dependent escape path length for antiprotons in the galaxy [14] is used:(2)λesc=11.8βforR<4.9GV,λesc=11.8β(R/4.9)0.54forR≥4.9GV,
the interaction length λinel of antiprotons including annihilation is also simulated with the MSDM and the stopping power 〈dE/dx〉 is calculated utilizing standard procedure (e.g., [15]).Equations (1a), (1b), and (1c) is solved through an iteration procedure using the “Mathematica” package. The solution readily converges in the third iteration. As shown in Figure 4 the MSDM cross section provides about two-times larger tertiary output Q3p̃  in the range of 0.3–3 GeV as compared to those from the uniform distribution used in [12]. In the energy range of 0.04–2 GeV the LIS obtained with MSDM slightly exceeds that obtained with the approximation [12] with the maximum deviation of ≤40% atEp̃=0.2 Gev.


## 2.3. Solar Modulation of the LIS

The LIS modulation in the heliosphere is considered on the basis of a transport equation for the spherically symmetric case [16]. The “force field” approximation inherently neglects the heliospheric gradient and curvature drifts but considers the diffusion, convection, and adiabatic deceleration:(3)F1AU(E1AU)P1AU2=FHB(EHB)PHB2,E1AU=Pc1nP1AU+E1AUPc+Ec+Ec+ΦforP<Pc,E1AU=EHB+ΦforP<Pc,Φ=V(RHB-1AU)3A,E1AU2=P1AU2+m02,EHB2=PHB2+m02.HereFHB, EHB, PHB, and F1AU, E1AU, and  P1AU  are the proton flux, proton total energy in GeV, momentum in GV at heliospheric boundary (HB) and at the Earth’s orbit (1 AU), respectively, m0 is the proton rest mass in GeV. V is the average solar wind speed in 103 km/hr. Physical sense of the solution implies a conservation of the distribution function F/P2 for particle energy decreases from EHB down to E1AU in travel from heliosphere RHB to the Earth at 1 AU.The heliospheric conditions are described by the “force field” parameterΦ determined by the solar wind speed V and the heliospheric boundary distance RHB. In our simulation we used A=17, Pc=1.015GV [17], where it was shown that for the energies ≥0.02 GeV  (2) approximates the exact solution of the equation in [16] for the power proton spectrum ∝βEHB-γ (where β=PHB/EHB  is the proton speed) and also consistent with the solar flare proton observations.TheΦ magnitude is determined from the best fit approximation with (2) of the observed proton spectrum F1AU assuming the interstellar spectrum as FHB(EHB)=16470βEHB-2.76protons/m s sr GeV. The fits results are Φmax=0.964 GeV and Φmin=0.368 GeV corresponding to solar maximum and minimum epochs.The results of our simulated spectrum are shown in Figure5 together with the antiproton fluxes obtained in different experiments performed at different solar minimum and maximum periods. The increases in the low-energy fluxes are provided by the higher fluxes of more energetic particles enriching the <1 GeV region due to adiabatic energy losses. The steeper the low-energy branch of the LIS spectrum the more pronounced the above-mentioned increases and [18]. The results of the simulations provided flux values of ≈4·10−3 to 10−2 and 10−2 to 1.7·10−2  antiprotons/m2 s sr GeV at energies of 0.2 and 1 GeV, respectively, corresponding to the solar maximum and minimum epochs. The curve for Φ=1.5 is the lower limit for all the experimental data. It may correspond, for example, to V=103 km/hour and RHB=70AU.Figure 5
Simulated LIS and the modulated spectra compared with experimental observations during the different phase of solar cycles: BESS95+97 [19]; BESS 00 [20], BESS 02 [21], BESS-Polar II [22], and PAMELA [4]. “Force field” parameter Φ in GeV.


## 3. Geomagnetically Trapped Antiproton Fluxes

Another source of low-energy antiprotons in the magnetosphere is interaction of GCR with the Earth’s residual atmosphere within geomagnetic trap. The source is quite the same as for the GCR antiprotons, but specific characteristics of the magnetic trap result in a significantly softer spectrum [23, 24].The flux captured within the plane of the geomagnetic equator can be simulated in frames of the radial diffusion model with an equation [25]:(4)Q2p̃(L,Ep̃)+Q3p̃(L,Ep̃)=fλinel-L2∂∂L(1L2DLL(L,Ep̃)∂f∂L)+1μ∂∂μ(μ〈dμdt〉f),
where μ is the antiproton magnetic moment, L is the McIlwain L-shell parameter in REarth,f(L,μ)=F(L,Ep̃)/P(Ep̃)2 represents the phase space distribution function (as in (2)), 〈dμ/dt〉 describes antiproton energy losses due to Coulomb scattering, λinel describes inelastic and annihilation interactions as in (1a), the radial diffusion coefficient DLL   is a sum of those of the magnetic DLLM and electric DLLE  diffusion:(5)DLL(L,μ)=DLLM+DLLE=7.10-9L10+10-4L10L10+μ2,REarth2s.The source functionsQ2p̃(L,Ep̃)and Q3p̃(L,Ep̃) are the same as in (1b) and (1c) but depend on L because of changing composition of the upper atmosphere [26] and on the cutoff rigidity modulating the parent CR proton and antiproton spectra. During the diffusion antiprotons born in situ at L-shell ≈1.2 sink in the underlying atmosphere or propagate to higher L-shells up to those where its rigidity is equal to that of the antiprotons and escape from the magnetosphere.The Finite Difference Method (FDM) is utilized for the solution of the partial differential equation (4). The results of the simulations for the direct reaction are presented in  Figure 6 for different antiproton energies. In the low-energy region the trapped population exceeds the CR antiproton one by a few orders of magnitude.Figure 6
TheL-dependence of the trapped antiproton flux spectra for various energies of antiprotons. The horizontal segments of the curves correspond to the direct penetration of the cosmic ray flux to the magnetosphere.Equation (4) accounts for only the antiprotons produced in reaction in situ and trapped in the same place. One more source of the antiproton belt stems from decay of antineutrons born in another channel of the pair production reaction: (p+p,n+ñ+p+p) of GCR with the Earth’s atmosphere constituents (this process was first considered in detail by [2] under a project supported with NASA’s grant for search for effective space rocket fuel). In contrast to antiprotons, antineutrons owing to their zero charge travel far from their birthplace. They follow primarily  the parent proton direction, but partially are backscattered in the atmosphere. Antiprotons resulted from decay of these albedo antineutrons acts as a supply for antiparticle radiation belts surrounding the Earth.This is the same CRAND (Cosmic Ray Albedo Neutron Decay) mechanism forming innermost proton radiation belt but with the antineutrons instead of the neutrons CRAÑD Due to that the CRAÑD flux can be calculated through normalization of the empirical CRAND model data by the antiproton-to-proton ratio for the antiproton source function: 1 antineutron for 105 up to 109        protons in dependence on the energy. This provides a few orders excess of this trapped population over the GCR antiproton flux in the low-energy region. The belt is replenished every few years.L dependence of the antiproton flux from CRAÑD source is shown in Figure 7. Maxima of the fluxes are located at the L=1.3-1.6 [27] that is slightly higher being compared with that from the direct reaction (L=1.25 in Figure 6); due to that this population was named as the “external antiproton belt” [23] due to its spatial position relative to the “inner belt” from the “in situ source” described above.Figure 7
TheL-dependence of the trapped antiproton fluxes spectra in the Earth’s magnetosphere from the CRAÑD source for various energies.The flux magnitudes simulated for the PAMELA orbit in the South Atlantic Anomaly (SAA) region are compared with the experimental PAMELA results [28] and simulation [29] for L=1.25 in Figure 8. The SAA region where the PAMELA trapped antiproton flux was measured (L=1.15-1.4) is characterized by strong L and pitch-angle dependence of the trapped flux. Due to that the result of simulation there is very sensitive to small variation of the model parameters used. This may explain the difference between our simulation and that of [29].Figure 8
Comparison of results of simulations with data of PAMELA experiment:ε is the criterion of stable trapping, P is the particle moment and B is geomagnetic field strength module.


## 4. Conclusions

The present simulation of the fluxes of interstellar origin incorporating solar modulation is attempted to explain the recent measurements of antiprotons at solar maximum and minimum. Particularly for the possible excess of the <1 GeV interstellar antiproton observations, initially the simulation considered the tertiary and antineutron decay antiprotons of the LIS source. The interaction cross sections by the MSDM Monte Carlo code provided a slightly larger antiproton flux in the energy range of 0.1–1 GeV compared to the [12] approximation. Then the “force field” solution for the solar modulation with rigidity dependence in compliance with the LIS and the 1 AU spectra showed satisfactory agreement between the simulations and the balloon results at the solar maximum and minimum periods.The simulation of geomagnetically trapped antiproton fluxes generated by GCR in the atmosphere showed a few orders of magnitude higher flux atL=1.2 regions being compared with the interstellar antiprotons at sub-GeV energies. This agrees well with the last PAMELA experiment [28].



---

*Source: 471094-2012-02-13.xml*