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 of 4103 to 102 and 102 to 1.7102 antiprotons/(𝑚2 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.

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) 𝐹̃𝑝(𝐸̃𝑝)was obtained through the leaky-box model as a solution of the following system of integro-differential equations [6]:𝐹̃𝑝𝐸̃𝑝𝜆esc+𝐹̃𝑝𝐹̃𝑝𝜆inel+𝑑𝑑𝐸̃𝑝𝐹̃𝑝𝐸̃𝑝𝑑𝐸𝐸𝑑𝑥̃𝑝=𝑄2𝐸̃𝑝̃𝑝+𝑄3𝐸̃𝑝̃𝑝,𝑄(1a)2𝐸̃𝑝̃𝑝=𝑁𝐴vog𝑗=𝐻,He,𝑜𝑛𝑗𝐴𝑗𝐸𝑝,𝑡𝑑𝜎𝑗𝑑𝐸̃𝑝𝐸̃𝑝,𝐸̃𝑝,𝐹̃𝑛𝑃𝐸𝑃𝑑𝐸𝑃,𝑄(1b)3𝐸̃𝑝̃𝑝=𝑁𝐴vog𝑗=𝐻,He,𝑜𝑛𝑗𝐴𝑗0𝑑𝜎𝑗𝑑𝐸̃𝑝𝐸̃𝑝,𝐸̃𝑝𝐹̃𝑝𝐸̃𝑝𝑑𝐸̃𝑝.(1c) This considers the production of the secondary antiprotons 𝑄2̃𝑝 by galactic cosmic ray (GCR) proton flux 𝐹𝑝(𝐸𝑝) and subsequent tertiary antiprotons 𝑄3̃𝑝 energy losses 𝑑𝐸/𝑑𝑥 in the interstellar matter, the flux decreases due to escape (𝜆esc) and inelastic interaction (𝜆𝑖𝑛𝑒𝑙).

2.1. Antiproton Production Spectrum

The antiproton production spectrum, that is, the source function 𝑄2̃𝑝 + 𝑄3̃𝑝 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 𝑛𝑗 are 1, 0.1, 8·10-4 cm−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 (𝑥,𝐴) of an incident particle 𝑥 with a target nucleus 𝐴 together with total and inelastic cross sections and multiplicities. The projectile can be a hadron (𝑛,̃𝑛,𝑝,̃𝑝) or a meson (𝜋+,𝜋,𝜋0,𝐾+,𝐾,𝐾0) with kinetic energies from 10 MeV up to 1 TeV. The target can be any nucleus with the atomic mass 𝐴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 energies 𝐸<1GeV are simulated on the base of parameterization of experimental data and πΔ-dynamics, for energies 1GeV𝐸10GeV the improved version of the quark gluon string (QGS) model is used, and for 𝐸>10 GeV the model of independent QGS is used.

Figure 1 shows differential cross sections of the antiprotons 𝑑𝜎𝐻(𝐸𝑝,𝐸̃𝑝)/𝑑𝐸𝑝 and antineutrons 𝑑𝜎𝐻(𝐸𝑝,𝐸̃𝑛)/𝑑𝐸𝑝 for proton projectiles of various energies 𝐸𝑝 in (𝑝,𝐻) reaction (scattering of the data in Figures 14 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).

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 (𝑝,𝐻𝑒) and (𝑝,𝑂) 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 (𝑝,𝐻) cross sections, and multiplicities are shown in Figure 2.

Figure 3 presents examples of spectra of the antiprotons𝑑𝜎𝐻(𝐸̃𝑝,𝐸̃𝑝)/𝑑𝐸̃𝑝, and antineutrons 𝜎𝐻(𝐸̃𝑝,𝐸̃𝑝)/𝑑𝐸̃𝑝, produced by antiproton projectiles of various energies 𝐸̃𝑝 in (̃𝑝,𝐻) reactions. The spectra noticeably differ from the uniform distributions ∝1/projectile energy from [12]. In contrast with (𝑝,𝐻) reaction, which produces practically identical spectra of ̃𝑝 and ̃𝑛 at all the projectile energies those spectra from (̃𝑝,𝐻) 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 approximation 𝐹𝑝(𝐸𝑝)=1.37104𝐸𝑝2.732· (m2sr s GeV)−1 for protons and 𝐹He(𝐸He)=7.06103𝐸He2.699· (m2sr s GeV)−1 for He [13] are used. The production spectra 𝑄2̃𝑝 obtained from 𝐹𝑝(𝐸𝑝) are depicted in Figure 4. The cross section invariance for interaction of the GCR He with interstellar 𝐻 was treated as a convolution of He LIS spectrum with 𝑑𝜎He(𝐸𝑝,𝐸̃𝑝)/𝑑𝐸𝑝, where for He the 𝐸𝑝 is expressed as energy per nucleon.

2.2. Antiproton LIS

In the simulation, the rigidity (𝑅) dependent escape path length for antiprotons in the galaxy [14] is used:𝜆esc𝜆=11.8𝛽forR<4.9GV,esc=11.8𝛽(R/4.9)0.54forR4.9GV,(2) the interaction length 𝜆inel of antiprotons including annihilation is also simulated with the MSDM and the stopping power 𝑑𝐸/𝑑𝑥 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 𝑄3̃𝑝 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% at𝐸̃𝑝=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:𝐹1AU𝐸1AU𝑃21AU=𝐹HB𝐸HB𝑃2HB,𝐸1AU=𝑃𝑐𝑃1n1AU+𝐸1AU𝑃𝑐+𝐸𝑐+𝐸𝑐+Φfor𝑃<𝑃𝑐,𝐸1AU=𝐸HB+Φfor𝑃<𝑃𝑐,𝑅Φ=𝑉HB1AU,𝐸3𝐴21AU=𝑃21AU+𝑚20,𝐸2HB=𝑃2HB+𝑚20.(3)

Here 𝐹HB, 𝐸HB, 𝑃HB, and 𝐹1AU, 𝐸1AU, and 𝑃1AU 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, 𝑚0 is the proton rest mass in GeV. 𝑉 is the average solar wind speed in 103 km/hr. Physical sense of the solution implies a conservation of the distribution function 𝐹/𝑃2 for particle energy decreases from 𝐸HB down to 𝐸1AU in travel from heliosphere 𝑅HB to the Earth at 1 AU.

The heliospheric conditions are described by the “force field” parameter Φ determined by the solar wind speed 𝑉 and the heliospheric boundary distance 𝑅HB. In our simulation we used 𝐴=17, 𝑃𝑐=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 𝛽𝐸HB𝛾 (where 𝛽=𝑃HB/𝐸HB 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 𝐹1AU assuming the interstellar spectrum as 𝐹HB(𝐸HB)=16470𝛽𝐸HB2.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 Figure 5 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 𝑉=103 km/hour and 𝑅HB=70AU.

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]:𝑄2̃𝑝𝐿,𝐸̃𝑝+𝑄3̃𝑝𝐿,𝐸̃𝑝=𝑓𝜆inel𝐿2𝜕1𝜕𝐿𝐿2𝐷𝐿𝐿𝐿,𝐸̃𝑝𝜕𝑓+1𝜕𝐿𝜇𝜕𝜕𝜇𝜇𝑑𝜇𝑓,𝑑𝑡(4) where 𝜇 is the antiproton magnetic moment, 𝐿 is the McIlwain 𝐿-shell parameter in 𝑅Earth,𝑓(𝐿,𝜇)=𝐹(𝐿,𝐸̃𝑝)/𝑃(𝐸̃𝑝)2 represents the phase space distribution function (as in (2)), 𝑑𝜇/𝑑𝑡 describes antiproton energy losses due to Coulomb scattering, 𝜆inel describes inelastic and annihilation interactions as in (1a), the radial diffusion coefficient 𝐷𝐿𝐿 is a sum of those of the magnetic 𝐷𝑀𝐿𝐿 and electric 𝐷𝐸𝐿𝐿 diffusion:𝐷𝐿𝐿(𝐿,𝜇)=𝐷𝑀𝐿𝐿+𝐷𝐸𝐿𝐿=7.109𝐿10+104𝐿10𝐿10+𝜇2,𝑅2Earth𝑠.(5)

The source functions 𝑄2̃𝑝(𝐿,𝐸̃𝑝)and 𝑄3̃𝑝(𝐿,𝐸̃𝑝) are the same as in (1b) and (1c) but depend on 𝐿 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 𝐿-shell 1.2 sink in the underlying atmosphere or propagate to higher 𝐿-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.

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: (𝑝+𝑝,𝑛+̃𝑛+𝑝+𝑝) 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.

𝐿 dependence of the antiproton flux from CRAÑD source is shown in Figure 7. Maxima of the fluxes are located at the 𝐿=1.31.6 [27] that is slightly higher being compared with that from the direct reaction (𝐿=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.

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 𝐿=1.25 in Figure 8. The SAA region where the PAMELA trapped antiproton flux was measured (𝐿=1.151.4) is characterized by strong 𝐿 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].

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 at 𝐿=1.2 regions being compared with the interstellar antiprotons at sub-GeV energies. This agrees well with the last PAMELA experiment [28].