Modeling of the Near-Earth Low-Energy Antiproton Fluxes

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 of 4⋅10⁻³ to 10⁻² and 10⁻² to 1.7⋅10⁻² antiprotons/(𝑚² 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+𝑑=𝑄₂𝑑𝐸_̃𝑝𝐹^̃𝑝𝐸_̃𝑝𝑑𝑑𝐸𝑥𝐸^̃𝑝𝑄𝑝̃𝐸^̃𝑝+𝑄₃𝑝̃𝐸^̃𝑝,(1a)₂𝑝̃𝐸^̃𝑝=𝑁_𝐴vog_{𝑗=𝐻,He,𝑜}𝑛_𝑗𝐴_𝑗^∞𝐸_𝑝,𝑡ℎ𝑑𝜎_𝑗𝑑𝐸_̃𝑝𝐸_̃𝑝,𝑄𝐸_̃𝑝,𝑛̃𝐹^𝑃𝐸_𝑃𝑑𝐸_𝑃,(1b)₃𝑝̃𝐸^̃𝑝=𝑁_𝐴vog_{𝑗=𝐻,He,𝑜}𝑛_𝑗𝐴_𝑗^∞0𝑑𝜎_𝑗𝑑𝐸_̃𝑝𝐸^̃𝑝,𝐸^′_̃𝑝𝐹^̃𝑝𝐸^̃𝑝𝑑𝐸^̃𝑝.(1c) This considers the production of the secondary antiprotons 𝑄₂𝑝̃ by galactic cosmic ray (GCR) proton flux 𝐹_𝑝(𝐸_𝑝) and subsequent tertiary antiprotons 𝑄₃𝑝̃ 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 𝑄₂𝑝̃ + 𝑄₃𝑝̃ 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^-4cm⁻³ [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 (𝜋⁺,𝜋⁻,𝜋⁰,𝐾⁺,𝐾⁻,𝐾⁰) 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 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 (𝑝,𝐻) reaction. Symbols mark the MSDM results; thick solid lines represent the Tan and Ng [12] approximation.

Figure 2 

Inelastic, total (𝑝,𝐻) cross sections, and multiplicities simulated with the MSDM code compared with the approximation from [12].

Figure 3 

Antiproton and antineutron differentia cross-section of 𝑝̃+𝐻 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 (𝑝,𝐻𝑒) 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.3710⁴𝐸_𝑝−2.732· (m²sr s GeV)⁻¹ for protons and 𝐹_He(𝐸_He)=7.0610³𝐸_He^−2.699· (m²sr s GeV)⁻¹ for He [13] are used. The production spectra 𝑄₂𝑝̃ 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.54forR≥4.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 𝑄₃𝑝̃ 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𝑃²¹_AU=𝐹^HB𝐸_HB𝑃²_HB𝐸_1AU,=𝑃_𝑐1n𝑃^1AU+𝐸_1AU𝑃_𝑐+𝐸_𝑐𝐸_1AU+𝐸_𝑐+Φfor𝑃<𝑃_𝑐,Φ==𝐸𝑉_HB𝑅+Φ_HBfor𝑃<𝑃_𝑐,−1AU𝐸3𝐴_21AU,=𝑃²¹_AU+𝑚²⁰,𝐸²_HB=𝑃²_HB+𝑚²⁰.(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, 𝑚₀ is the proton rest mass in GeV. 𝑉 is the average solar wind speed in 10³ km/hr. Physical sense of the solution implies a conservation of the distribution function 𝐹/𝑃² 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𝛽𝐸_HB^−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 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⁻³ to 10⁻² and 10⁻² to 1.7·10⁻²  antiprotons/m² 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 𝑉=10³ km/hour and 𝑅_HB=70AU.

Figure 5 

Simulated LIS and the modulated spectra compared with experimental observations during the different phase of solar cycles: BESS 95+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]:𝑄₂𝑝̃𝐿,𝐸^̃𝑝+𝑄₃𝑝̃𝐿,𝐸^̃𝑝=𝑓𝜆_inel−𝐿²𝜕𝜕𝐿𝐿1²𝐷_𝐿𝐿𝐿,𝐸_̃𝑝𝜕𝑓+1𝜕𝐿√𝜇𝜕𝜕𝜇√𝜇𝑑𝜇𝑑𝑡𝑓,(4) where 𝜇 is the antiproton magnetic moment, 𝐿 is the McIlwain 𝐿-shell parameter in 𝑅_Earth,𝑓(𝐿,𝜇)=𝐹(𝐿,𝐸_̃𝑝)/𝑃(𝐸_̃𝑝)² 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.10⁻⁹𝐿¹⁰+10⁻⁴𝐿¹⁰𝐿¹⁰+𝜇²,𝑅^2Earth𝑠.(5)

The source functions 𝑄₂𝑝̃(𝐿,𝐸_̃𝑝)and 𝑄₃𝑝(̃𝐿,𝐸^̃𝑝) 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.

Figure 6 

The 𝐿-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: (𝑝+𝑝,𝑛+𝑛̃+𝑝+𝑝) 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 10⁵ up to 10⁹ 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.3−1.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.

Figure 7 

The 𝐿-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 𝐿=1.25 in Figure 8. The SAA region where the PAMELA trapped antiproton flux was measured (𝐿=1.15−1.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].

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