# Modulation of cosmic ray anti-protons in the heliosphere: simulations for a
solar cycle

O.P.M. Aslam Shandong Institute of Advanced Technology (SDIAT)
250100 Jinan, Shandong Province
People’s Republic of China  M.S. Potgieter Shandong Institute of Advanced
Technology (SDIAT)
250100 Jinan, Shandong Province
People’s Republic of China  Institute for Experimental and Applied Physics
(IEAP)
Christian-Albrechts-University in Kiel
24118 Kiel, Germany  Xi Luo (罗熙) Shandong Institute of Advanced Technology
(SDIAT)
250100 Jinan, Shandong Province
People’s Republic of China  M.D. Ngobeni Centre for Space Research, North-West
University
2520 Potchefstroom, South Africa

###### Abstract

The precision measurements of galactic cosmic ray protons from PAMELA and AMS
are reproduced using a well-established 3D numerical model for the period July
2006 - November 2019. The resulting modulation parameters are applied to
simulate the modulation for cosmic antiprotons over the same period, which
includes times of minimum modulation before and after 2009, maximum modulation
from 2012 to 2015 including the reversal of the Sun’s magnetic field polarity,
and the approach to new minimum modulation in 2020. Apart from their local
interstellar spectra, the modulation of protons and antiprotons differ only in
their charge-sign and consequent drift pattern. The lowest proton flux was in
February-March 2014, but the lowest simulated antiproton flux is found to be
in March-April 2015. These simulated fluxes are used to predict the proton to
anti-proton ratios as a function of rigidity. The trends in these ratios
contribute to clarify to a large extent the phenomenon of charge-sign
dependence of heliospheric modulation during vastly different phases of the
solar activity cycle. This is reiterated and emphasized by displaying so-
called hysteresis loops. It is also illustrated how the values of the parallel
and perpendicular mean free paths, as well as the drift scale, vary with
rigidity over this extensive period. The drift scale is found to be at its
lowest level during the polarity reversal period, while the lowest level of
the mean free paths are found to be in March-April 2015.

Particle astrophysics (96) — Cosmic rays (329) — Heliosphere (711) — Solar
cycle (1487) — Quiet Sun (1322)

##  1 Introduction

The level and extent of the occurrence of charge-sign dependent modulation in
the heliosphere has always been a puzzling aspect because of the continuously
varying nature and complexity of the Sun’s activity cycles. This has been
exasperated by the unavailability of continuous and long-term observations of
galactic cosmic rays (GCRs) of opposite charge-signs with good precision and
over an extended range of rigidity. This aspect, in particular, has changed
markedly with the deployment of space borne experiments such as PAMELA, the
Payload for Antimatter Matter Exploration and Light-nuclei Astrophysics
(Adriani et al., 2017) in 2006 and AMS, the Alpha Magnetic Spectrometer
(Aguilar et al., 2021a) on the International Space Station in 2011. Before
these missions, charge-sign dependent modulation was observed by the Ulysses
mission (Heber et al., 1999, 2002, 2003), which explored the heliosphere for
the first time at high helio-latitudes, including two solar minima and a solar
maximum period but separate measurements of electrons and positrons could not
be made; see https://sci.esa.int/web/ulysses and the review on GCR
observations by Heber & Potgieter (2006); Heber (2013). The two Voyager
missions could also not distinguish between electrons and positrons but
contributed significantly to our knowledge of the outer regions and boundary
of the heliosphere, from the termination shock (TS) to the heliopause (HP),
and when crossing the latter provided the first direct measurements of what
can be considered the very local interstellar spectra (VLIS’s) of various GCR
species (Stone et al., 2013, 2019; Cummings et al., 2016), also see
https://voyager.jpl.nasa.gov/. However, these reported VLIS’s do not include
those for positrons and antiprotons. For a recent review on the Voyager
observations, and on other related GCR observations, see Rankin et al. (2022);
for summaries of various aspects of charge-sign dependent observations by
PAMELA, see e.g., Boezio et al. (2017); Di Felice et al. (2017); Munini et al.
(2017).

Even for the PAMELA and AMS missions reporting positron spectra and
particularly anti-proton spectra over short time intervals and over a wide
rigidity range remain an effort. Knowledge of the long-term modulation of
these particular GCRs over solar cycles is therefore not at the same level
than for protons and most other GCRs. However, numerical modulation models
have become quite sophisticated and can easily be utilized to study the long-
term modulation of these particles. Such a study for positrons was performed
by Aslam et al. (2019) who discussed the challenge of what to assume for the
positron VLIS and what exactly to use as modulation parameters over the long-
term. Following their approach, we simulated the proton spectra first and then
study the modulation of anti-protons from one solar minimum to the next,
including a period of the reversal of the magnetic field direction of the Sun.
We find insightful results which are useful in making conclusions about the
extent of charge-sign dependence over a complete solar cycle. For this we
utilize the reported galactic proton measurements from PAMELA, combined with
that from AMS, for the period from July 2006 to November 2019 (Adriani et al.,
2013a; Aguilar et al., 2018; Martucci et al., 2018).

In what follows, we describe the specific objectives, the utilized long-term
observations, our modulation modelling approach and assumptions made, the
model itself, the numerical solutions (results) for protons and anti-protons
and the calculated corresponding flux ratios over the mentioned period.

##  2 Objectives and Observations

Our first objective is to simulate and reproduce GCR proton spectra at the
Earth. We utilize Carrington rotation averaged spectra for the period July
2006 - May 2011 (Adriani et al., 2013a) and from January 2010 - February 2014
(Martucci et al., 2018) over the kinetic energy range 80 MeV to 50 GeV
measured by PAMELA. For a detailed description of the PAMELA detector and
mission, see e.g., Adriani et al. (2017). It is noted that the highest proton
spectra ever measured were reported by Adriani et al. (2013a) for the end of
2009. In this context, see also Potgieter et al. (2014); Krainev et al.
(2021).

We further utilize proton spectra from AMS from May 2011 - May 2017 (BRs
2426-2506) as published by Aguilar et al. (2018) and also for May 2011 -
November 2019 (BRs 2426-2540) as published by Aguilar et al. (2021b). For the
later period, up to November 2019, BR averages are calculated from the
reported daily resolution measurements over the rigidity range 1 - 100 GV. For
a detailed description about the AMS experiment, see e.g., Aguilar et al.
(2021a). Combining the AMS and PAMELA proton spectra over the period July 2006
- November 2019 provides us with more than 13 years of proton spectra as basis
for our observational motivated and data driven simulations.

Our second objective is to simulate the corresponding antiproton modulation
for the period mentioned above. In order to accomplish this, all the
modulation parameters found with the numerical model in reproducing the whole
set of proton spectra are used unchanged to simulate the antiproton spectra
for the exact same BRs. The only difference between the modulation of protons
and antiprotons in the heliosphere is their drift patterns, oppositely
directed during the same polarity phase, and, of course, their different
VLIS’s which are discussed in section 4.2.

A third objective is to use the simulated proton and anti-proton spectra to
calculate their ratios at various rigidity values to establish how they change
over time. Detailed time variations of the proton flux at 1.0 - 1.16 GV over
the mentioned periods are explored. This provides insight of how particle
drift and consequently charge-sign dependent modulation evolves over time with
vastly different phases of the solar cycle.

##  3 Solar activity proxies

The continuously varying nature of the Sun’s activity (Usoskin, 2013;
Hathaway, 2015) creates an 11-year time-dependence in the modulation of GCRs,
and with the Sun’s magnetic field changing its direction during each solar
maximum, a 22-year cycle becomes evident in how GCRs are modulated. In this
context, the Sun’s magnetic ’polarity’ is named positive (A $>$ 0) when the
field lines are outward in its northern hemisphere and inward in the southern
hemisphere. When the field lines change to be outward in the southern
hemisphere and inward in the northern hemisphere, the ’polarity’ is named
negative (A $<$ 0). During e.g., the positive phase, positively charged GCRs
are assumed to drift towards the Earth mainly, but not exclusively, through
the polar regions of the heliosphere whereas negatively charged GCRs drift
inwards mostly through the equatorial regions for the heliosphere. They then
encounter the wavy heliospheric current sheet (HCS). See the reviews by e.g.,
Heber (2013); Kóta (2013) and Potgieter (2014, 2017).

During Solar Cycle 24, the change from $A<0$ to $A>0$ occurred from November
2012 - March 2014. This reversal is considered a relatively slow process and
get completed in $\sim$16 months, see Sun et al. (2015). For our modeling
purposes we consider this reversal as a period with no well-defined polarity
in neither the northern nor the southern hemisphere of the Sun.

The time (solar) related activity is introduced to the model by using the tilt
angle, $\alpha$, of the HCS along with the magnitude, $B$, of the heliospheric
magnetic field (HMF) at the Earth as external parameters; both are considered
good proxies for changing solar activity levels.

Figure 1: Lower panel: Magnitude of the HMF observed at the Earth averaged
over Bartels Rotations (BR), adopted from http://omniweb.gsfc.nasa.gov for
January, 2005 to December, 2019 (BRs 2340 to 2542) shown as solid line,
together with $\sim$15-months (17-BRs) running averages, shown as red circles
and $\sim$10-months (11-BRs) running averages, shown as blue circles. Upper
panel: Tilt angle ($\alpha$) of the Heliospheric Current sheet (HCS) as
Carrington Rotation (CR) averages adopted from http://wso.stanford.edu/ for
the same period (CRs 2025 to 2226) shown as a solid line together with
$\sim$15-months (17 CRs) running averages shown as red circles. Shaded parts
indicate the period with no well-defined magnetic field polarity (November
2012 - March 2014), when changing cycles from $A<0$ to $A>0$.

In Figure 1, the variation of $\alpha$ is shown in the top panel from January
2005 - December 2019, averaged over Carrington Rotations from
http://wso.stanford.edu/, along with calculated $\sim$15-months (17 rotations)
running averages as implemented in the numerical model. In the lower panel,
the observed $B$ at the Earth are averaged over Bartels Rotations (BRs), from
http://omniweb.gsfc.nasa.gov, together with calculated $\sim$15 months (17
BRs) running averages and $\sim$10 months (11 BRs) running averages, as
implemented in the numerical model. These averages reflect the assumed
propagation times of these entities from the Earth to the HP. For the HMF, the
$\sim$10-month-averages are used for the solar minimum of July 2006 - October
2010 and again for August 2017 - November 2019, when the averaged $\alpha$ is
below 20∘, while the $\sim$15-month-averages are used for other periods.
Selecting these running averages to be used as solar activity proxies in our
modeling approach is explained and discussed by Potgieter et al. (2015); Aslam
et al. (2019, 2021).

Another solar activity related feature changing with time in the model is the
heliolatitude dependence of the solar wind speed on a global scale. In
accordance with Ulysses observations, it is assumed that it changes from 430
km s-1 in the equatorial plane to 750 km s-1 in the polar regions during solar
minimum activity conditions, see McComas et al. (2002); Heber & Potgieter
(2006). As the Solar Cycle progresses towards maximum, this latitudinal
dependence gradually decreases in the polar regions to become 450 km s-1 when
$\alpha$ is 65∘ or higher.

The width of the inner heliosheath affects the modulated GCR intensity at the
Earth (Vos & Potgieter, 2016). The TS moves inward and outward in response to
changes in the dynamic pressure of the solar wind, and that changes the width
of the inner heliosheath, as reported by Richardson & Wang (2011) based on the
Voyager 2 measurements. This dynamic nature of the TS is incorporated in the
model, in addition to changing $\alpha$ and the HMF magnitude; for similar
approaches, see Langner & Potgieter (2004, 2005); Manuel et al. (2014). In the
model, the position of the TS is changed gradually from 88 au (July 2006) to
80 au (December 2009), gradually increasing back to 88 au by April - May 2014,
and then gradually decreasing to 78 au by mid-2019; see Figure 3 by Vos &
Potgieter (2016). The spherically shaped HP is assumed to be at a radial
distance of 122 au from the Sun.

##  4 Modeling Approach

Once inside the heliosphere GCRs are subjected to an outward convection by the
solar wind, adiabatic deceleration due to the expanding solar wind, diffusion
along and across the turbulent HMF, and gradient, curvature and HCS drifts in
the global HMF. These processes are responsible for altering the differential
intensities of GCRs as a function of energy (rigidity), position in the
heliosphere and time (solar activity). We simulate the transport of GCRs
within the heliosphere by including all of these modulation processes into a
full three-dimensional (3D) numerical model to compute the differential
intensity of GCRs at different radial distances from the Earth up to the HP
and spectra at the Earth over the rigidity range 10 MV - 70 GV. The VLIS’s are
specified as an initial condition (unmodulated spectra) at the HP. The
predicted small modulation beyond the HP, as simulated by e.g., Scherer et al.
(2011); Luo et al. (2015), is neglected.

###  4.1 Modulation model

The model is based on the numerical solution of Parker’s transport equation
(Parker, 1965) which follows the basic transport equation (TPE) for charged
particle motion in large and small scale fluctuating magnetic fields, and
averaged over the pitch and phase angles of propagation particles, with the
assumption that GCRs are approximately isotropic. For the heliosphere this TPE
is described as:

$\frac{\partial f}{\partial
t}=-\vec{V}_{sw}\cdot\nabla\it{f}-\langle\vec{v}_{D}\rangle\cdot\nabla\it{f}+\nabla\cdot(\bf{K}_{s}\cdot\nabla\it{f})+\frac{1}{3}(\nabla\cdot\vec{V}_{sw})\frac{\partial
f}{\partial\ln p}$ (1)

where $f(\vec{r},p,t)$ is the omnidirectional GCR distribution function,
$\it{p}$ is momentum, $\it{t}$ is time, and $\vec{r}$ is the vector position
in 3D. If $\it{P}$ is particle rigidity, the differential intensity is given
by $\it{P}^{2}f$ ($\it{p}^{2}f$). The first two terms shown on the right-hand
side of this equation represent outward convection caused by the expanding
solar wind with velocity ($\vec{V}_{sw}$) and the averaged particle drift
velocity $\langle\vec{v}_{D}\rangle$ (pitch angle averaged guiding center
drift velocity) i.e. particle drift caused by the HMF’s gradients, curvatures
and the wavy HCS. The third term represents spatial diffusion described by the
symmetry diffusion tensor $\bf{K}_{s}$. When the TPE is written in
heliocentric spherical coordinate system, the three coordinates $\it{r}$,
$\theta$, and $\phi$ represent radial distance, polar and azimuth angles,
respectively. According to this system, the Earth is located at 1 au in the
equatorial plane with $\theta=90^{\circ}$.

In order to solve the TPE, six elements in this tensor related to diffusion
are to be specified: $K_{rr}$, $K_{r\phi}$, $K_{\theta\theta}$, $K_{\phi r}$,
and $K_{\phi\phi}$. These five diffusion coefficients in terms of diffusion
parallel and perpendicular to the HMF direction are as follows:

$K_{rr}=K_{\parallel}{\cos}^{2}\psi+K_{\perp r}{\sin}^{2}\psi.$ (2)
$K_{\theta\theta}=K_{\perp\theta}.$ (3) $K_{\phi\phi}=K_{\perp
r}{\cos}^{2}\psi+K_{\parallel}{\sin}^{2}\psi.$ (4) $K_{\phi r}=(K_{\perp
r}-K_{\parallel}){\cos}\psi{\sin}\psi=K_{r\phi}.$ (5)

This means, for example, that the effective radial diffusion coefficient
$K_{rr}$ is a combination of the parallel diffusion coefficient
($K_{\parallel}$) and the radial perpendicular diffusion coefficient
($K_{\perp r}$), with $\psi$ the spiral angle of the background HMF;
$K_{\theta\theta}$ = $K_{\perp\theta}$ is the effective perpendicular
diffusion coefficient in the polar direction; $K_{\phi\phi}$ describes the
effective diffusion in the azimuth direction and $K_{\phi r}$ diffusion in the
$\phi r$-plane.

The fourth term in the right hand side of Equation 1 represents adiabatic
energy changes which depends on the sign of the divergence of $\vec{V}_{sw}$.
If $\nabla\cdot\vec{V}_{sw}>0$, adiabatic energy loss occurs as is the case in
most of the heliosphere, except perhaps inside the heliosheath where we assume
that $\nabla\cdot\vec{V}_{sw}=0$; see also Langner et al. (2006). Adiabatic
energy changes becomes dominant at lower rigidity in the case of modulated
spectra for protons, antiprotons and GCR isotopes; see illustrations of this
process by Moraal & Potgieter (1982); Strauss et al. (2011).

The drift velocity can be written as follows:

$\langle\vec{v}_{D}\rangle=\nabla\times K_{D}\frac{\vec{B}}{B},$ (6)

where $K_{D}$ is the generalized drift coefficient, and $\vec{B}$ is the HMF
vector with magnitude $\it{B}$. When Equation 1 is rewritten not showing the
drift velocity explicitly, then a generalized propagation tensor replaces the
diffusion tensor and four elements must then be added to describe gradient and
curvature drift: $K_{r\theta}$, $K_{\theta r}$, $K_{\theta\phi}$, and
$K_{\phi\theta}$. The three components of the drift velocity (in the radial,
polar and azimuth directions) in terms of these elements of the propagation
tensor are then given by:

$\langle\vec{v}_{D}\rangle_{r}=-\frac{A}{r\sin\theta}\frac{\partial}{\partial\theta}(\sin\theta
K_{\theta r}),\\\
\langle\vec{v}_{D}\rangle_{\theta}=-\frac{A}{r}\left[\frac{1}{\sin\theta}\frac{\partial}{\partial\phi}(K_{\phi\theta})+\frac{\partial}{\partial
r}(rK_{r\theta})\right],\\\
\langle\vec{v}_{D}\rangle_{\phi}=-\frac{A}{r}\frac{\partial}{\partial\theta}(K_{\theta\phi}),$
(7)

where $A$ is a constant that denotes the HMF polarity. In essence, it
determines the drift directions of charged particles:

$A=+1(qA>0);A=-1(qA<0),$ (8)

with $q$ the particle’s charge-sign.

Evidently, the geometry and magnitude of the HMF are important for the
modulation process, particularly in the case of particle drift. We start with
a straight-forward HMF as described by Parker (1958):

$\vec{B}=B_{0}A\Bigg{[}\frac{r_{0}}{r}\Bigg{]}^{2}(\hat{e}_{r}-\tan\psi\hat{e}_{\phi}),$
(9)

where $\hat{e}_{r}$ and $\hat{e}_{\phi}$ are unit vectors in the radial and
azimuth directions; $B_{0}$ is the magnitude of the HMF at the position
$r_{0}$, in our approach assumed to be at the Earth, and

${\tan}\psi=\Omega\frac{(r-r_{\odot})}{V_{sw}}\sin\theta,$ (10)

where $r_{\odot}$ is the solar radius (0.005 au) and $\Omega$ is the average
angular rotation speed of the Sun (2.66 $\times$ 10-6 rad s-1). The magnitude
of this Parkerian HMF is

$B=B_{0}\Bigg{[}\frac{r_{0}}{r}\Bigg{]}^{2}\sqrt{1+(\tan\psi)^{2}.}$ (11)

However, this HMF is modified as proposed by Smith & Bieber (1991) based on
their insight that the magnetic field spirals are relatively more tightly
wound than that predicted by the original Parker theory. They argued that the
differential rotation of the Sun would cause small azimuthal magnetic field
components to develop and that would lead to larger spiral angles at larger
radial distances. They proposed a modification so that the expression for the
HMF spiral angle becomes:

$\tan\psi=\frac{\Omega(r-r_{b})\sin\theta}{V_{sw}(r,\theta)}-\frac{rV_{sw}(r_{b},\theta)}{r_{b}V_{sw}(r,\theta)}\Bigg{(}\frac{B_{T}(r_{b})}{B_{R}(r_{b})}\Bigg{)}$
(12)

where $B_{T}(r_{b})/B_{R}(r_{b})$ is the ratio of the azimuthal to the radial
magnetic field components at a position $r_{b}$ near the solar surface. It is
argued that the ratio $B_{T}(r_{b})/B_{R}(r_{b})$ $\approx$ -0.02 at a
position $r_{b}$ = 20$r_{\odot}$. With $r_{\odot}$ = 0.005 au as the solar
radius, the value $r_{b}$ = 20$r_{\odot}$ and the ratio
$B_{T}(r_{b})/B_{R}(r_{b})$ = -0.02 are constants that determine the HMF
modification. This modification keeps the basic Parkerian geometry (and
$\nabla\cdot\vec{B}=0$) but modifies its magnitude progressively toward the
poles of the heliosphere, which essentially give less effective particle drift
in the polar regions of the heliosphere. The motivation for this modification
from a modulation modelling point of view is discussed by Potgieter (2013);
Potgieter et al. (2015); Potgieter & Vos (2017). An elaborate discussion and
illustration of its relevance and effects on GCR modulation is given by Raath
et al. (2015).

The functional forms of the three diffusion coefficients
($K_{\parallel}$,$K_{\perp r}$,$K_{\theta\theta}$) and the drift coefficient
($K_{D}$) as used in this model are given and discussed in section 5.1.

Figure 2: Very local interstellar spectra (VLIS’s) at the HP (122 au) as a
function of rigidity: for GCR protons as the black solid line from Bisschoff
et al. (2019) and for anti-protons as the blue solid line from Potgieter et
al. (2021). Observations at the Earth are shown as red triangles indicating
PAMELA protons, averaged for July 2006 - January 2010 (Adriani et al., 2013a),
with blue triangles indicating PAMELA antiproton averages for the same period
(Adriani et al., 2013b); blue circles are AMS antiproton averages for May 2011
- May 2015 (Aguilar et al., 2016) and the gray circles are AMS protons for the
same period (Aguilar et al., 2018).

Potgieter et al. (2014, 2015) described the global features of the model used
here in detail e.g. the expression for the full TPE in spherical coordinates;
for the expression of the global solar wind structure, see their equations 12
and 5, respectively. Raath et al. (2016) give the expression for the wavy HCS
as used here. Details are also discussed by Potgieter et al. (2015); Potgieter
& Vos (2017) and the same or a similar model was applied by Aslam et al.
(2019, 2021); Corti et al. (2019); Ngobeni et al. (2020, 2022). Comprehensive
reviews of this modeling approach, and underlying theory for the global
modulation of GCRs, are given by Potgieter (2014, 2017) and Potgieter et al.
(2021). As done for previous applications, we continue to assume a steady-
state so that the model is not well suitable to study GCR transient events as
explained by Aslam et al. (2021); the applied Alternating Direct Implicit
(ADI) numerical method is not stable when a fifth numerical dimension is
incorporated into the numerical scheme (that is, three spatial, one for energy
dependence, and one for time dependence). If all three spatial dimensions are
used, the time dependence must be sacrificed. The only alternative for 3D with
a full time-dependence is to use numerical models based on Stochastic
Differential Equations (SDEs) which is a totally different approach to
modulation studies, see e.g. Kopp et al. (2012); Dunzlaff et al. (2015); Luo
et al. (2019, 2020); Song et al. (2021); Shen et al. (2021), and references
therein. Many 2D time-dependent modeling results, based on an ADI numerical
scheme, have been published over the years e.g., le Roux & Potgieter (1992,
1995); Ferreira et al. (2003); Langner & Potgieter (2004, 2005); Ngobeni &
Potgieter (2012); Manuel et al. (2014).

###  4.2 Very Local Interstellar Spectra

In order to reproduce observed proton spectra, as vindication of the model,
and to simulate subsequently antiproton spectra, we have to specify for each
type (specie) of GCR particles an input spectrum as an initial condition in
the numerical model; it is done in the form of a VLIS at the HP. Bisschoff et
al. (2019, 2021) illustrate the approach of how to construct self-consistent
VLIS’s for GCRs by utilizing available comprehensive galactic propagation
models (GALPROP) along with a 3D modulation model and with relevant validating
observations. At high enough rigidity where modulation may be considered
negligible (Strauss & Potgieter (2014a)), observations from AMS and PAMELA are
used, whereas and at low rigidity for protons, observations from Voyager 1 and
Voyager 2 beyond the HP are used (Stone et al., 2013, 2019; Cummings et al.,
2016). See Strong & Moskalenko (1998) for information about GALPROP and a
description of the numerical computation of the propagation of primary and
secondary nucleons, and Moskalenko et al. (2002, 2005) especially for
numerical computation of antiproton propagation. Studies determining VLIS’s
for modulation modeling purposes, have been reported by e.g., Vos & Potgieter
(2015); Bisschoff & Potgieter (2016); Corti et al. (2016); Boschini et al.
(2019, 2020).

The VLIS for galactic protons obtained in this way, is depicted in Figure 2 in
units of ”particles m-2 sr-1 s-1 GV-1” as a function of rigidity over the
range 0.30 - 70.0 GV. (The Voyager data are not repeated here in these units;
see Figure 6 by Potgieter et al. 2021). The VLIS for antiprotons computed with
the GALPROP code and then adjusted as described by Bisschoff et al. (2021) and
Potgieter et al. (2021), is shown in comparison. Observational spectra for
protons and antiprotons at the Earth are shown for illustrative purposes and
described in the figure’s caption.

Apart from showing the VLIS’s applied in the modelling, this figure also
illustrates the total modulation obtained with the model between the HP (122
au) and Earth (1 au) and how vastly different that is for antiprotons compared
to protons. The peak intensity in the proton VLIS occurs at a much lower
rigidity than for antiprotons, subsequently also for the modulated spectra and
it shows that the additional modulation that protons experienced is markedly
more than for antiprotons based on the illustrative averaged intensity in 2006
- 2009 to that in 2011 - 2015. These differences are explored and reported in
what follows.

##  5 Results

We first simulate and reproduce proton spectra for each Carrington Rotation
(from PAMELA) or Bartels Rotation (from AMS) from July 2006 - November 2019
over a wide range of rigidity. Based on these spectra, the differential
intensity for the rigidity bin 1.0 GV - 1.16 GV is shown over time in Figure 3
and compared with observations for the same rigidity bin. Once the observed
proton spectra are reproduced satisfyingly, the same set of modulation
parameters is used to simulate the antiproton spectra for the same period of
time and rigidity as shown in Figure 3. Corresponding published antiproton
measurements are not available for comparison.

Figure 3: Simulated proton flux (blue solid line) from July 2006 - November
2019 (BRs 2361-2540) along with a standard error of mean (SEM; standard
deviation /$\sqrt{n}$ ) for 1.0 - 1.16 GV in comparison with the observed
proton flux (blue circles) of the same rigidity. For July 200 - May 2011
observations are from PAMELA (Adriani et al., 2013a; Martucci et al., 2018)
and for May 2011 - November 2019 from AMS (Aguilar et al., 2018, 2021b).
Corresponding simulated antiproton flux is shown (red solid line) along with
its SEM (gray shading). No published BRs averaged antiproton observations are
available for comparison. Shaded range indicates the period without a well-
defined HMF polarity (November 2012 - March 2014) as it changes from $A<0$ to
$A>0$.

Inspecting the time trends in Figure 3, it is noticed how both the proton and
antiproton fluxes show a gradual increase from July 2006 onward to reach the
highest level by the end of 2009. The proton flux shows $\sim$75% increase
over this period (based on the averaged value over the whole 2006 - 2019
period), but the antiproton flux shows an increase of only $\sim$17%. This
different trend in the intensity-time profiles is considered to be indicative
of the effect of different drift patterns of these oppositely charged
particle. After December 2009, both fluxes progressively decrease as the solar
cycle progresses towards maximum, with the proton flux reaching its minimum
level by the end of the polarity reversal phase (March - April 2014), but the
antiproton flux reaches a first minimum somewhat just after the reversal ended
and then another minimum in February - March 2015. This is related to the
peculiar behaviour of both the tilt angle and HMF in late 2014 as shown in
Figure 1. The rate of decrease of the antiproton flux from December 2009 to
March 2015 is different compared to protons, again considered to be indicative
of the different drift patterns and consistent to drift model predictions that
the negatively charged GCRs should have a flatter intensity-time profile
during $A<0$ cycles. Before the recovery starts in March - April 2014, the
proton flux then is only $\sim$20% of the flux in December 2009. The
antiproton flux is reduced up to only $\sim$50% of the December 2009 level
before the recovery start. After the HMF reversal, as the solar cycle
progresses toward minimum, both the fluxes increase, with the proton flux
reaching almost the same level in September-October 2019 as in December 2009.
Over the same period the antiproton flux at this rigidity reaches a level of
$\sim$95% of the level of December 2009.

###  5.1 Diffusion and particle drift parameters over a solar cycle

Figure 4: Panel (a): Variation of the drift scale over time at the Earth for
1.0-1.16 GV from July 2006 to October 2019. Panel (b): Variation of the three
MFPs $\lambda_{\parallel}$, $\lambda_{\perp r}$ and $\lambda_{\perp\theta}$
over the same period and rigidity. Both $\lambda_{\perp r}$ and
$\lambda_{\perp\theta}$ are multiplied by $0.5\times 10^{2}$ for a direct
comparison with $\lambda_{\parallel}$ indicating that they have the exact same
time dependence. In these panels the emphasis is on the trends over time from
before the 2009 solar minimum to the next.

In this section, we elaborate on the diffusion and drift coefficients used in
the model to reproduce the proton spectra and to simulate the antiproton flux.
These coefficients contain the essence of what determines the level of solar
modulation over a solar cycle. We emphasize first the variation of the drift
scale, with results shown to vary over time in Figure 4(a) at a chosen
rigidity, and then as a function of rigidity in Figure 5. The parallel and
perpendicular mean free paths are shown in Figure 4(b) at a chosen rigidity
over time, and also as a function of rigidity in Figure 5.

The expression for a generalized drift coefficient, can be written as:

$K_{D}=\frac{\beta P}{3B_{m}}f_{D}=\frac{\beta
P}{3B_{m}}\Bigg{[}\frac{(\omega\tau)^{2}}{1+(\omega\tau)^{2}},\Bigg{]}$ (13)

with $B_{m}$ the magnitude of the modified HMF (see section 4.1) and where
$\beta$ = $v/c$ is the ratio of particle speed to the speed of light, and
$\omega$ is the particle gyro-frequency with $\tau$ the average time between
the scattering of GCR particles in the HMF. In most of the numerical modeling
studies, it is assumed that $\omega\tau\gg 1$ in the heliosphere so that the
drift coefficient takes its simplest form:

$K_{D}=\frac{\beta P}{3B_{m}},$ (14)

known as weak scattering drift. Establishing $\tau$ is complicated, even
controversial; an elaborate turbulence theory is required to understand how
$\omega\tau$ could change as a function of rigidity and space throughout the
entire heliosphere, and especially over the solar cycle; see e.g., Moloto et
al. (2018); Zhao et al. (2018); Engelbrecht & Moloto (2021) and the review by
Engelbrecht et al. (2022). The term $f_{D}$ is called the drift reduction
factor and is determined by how diffusive scattering is described; if $f_{D}$
= 0, then $K_{D}$ and therefore the drift velocity $\langle v_{D}\rangle$
becomes zero so that drift effects vanish from the modulation model to produce
non-drift solutions; if $f_{D}$ approaches 1, drift is at a maximum, and
$K_{D}$ then has the weak scattering value. Theoretically, the largest
possible value for $K_{D}$ is obtained under such conditions and when the
Parker HMF is kept unmodified; for detailed discussions, see e.g., Ngobeni &
Potgieter (2015), and references there-in. If this is the case, drift effects
in GCR modulation are predicted to be very large and dominant as found in
original numerical models by Jokipii & Kopriva (1979, where the $A>0$ and
$A<0$ computed spectra were the opposite of what was observed), Jokipii &
Thomas (1981), Kóta & Jokipii (1983) and Potgieter & Moraal (1985), to mention
only a few. The function $f_{D}$ is also used to adjust the rigidity
dependence of $K_{D}$, which is the most effective direct way of suppressing
drift effects at low rigidity as required by $\it{Ulysses}$ observations of
latitudinal gradients; see reviews by Heber & Marsden (2001); Heber &
Potgieter (2006); Heber (2013). As such, Equation (13) gets a more practical
form (e.g. Ngobeni & Potgieter, 2015; Moloto et al., 2018) and has been widely
used in numerical drift models:

$K_{D}=\frac{\beta P}{3B_{m}}f_{D}=K_{A0}\frac{\beta
P}{3B_{m}}\frac{(P/P_{A0})^{2}}{1+(P/P_{A0})^{2}}.$ (15)

Here, $K_{A0}$ is dimensionless and used to adjust $K_{D}$ in a pragmatic way;
if $K_{A0}$ = 1.0, it is called 100% drift (full weak scattering). In this
study, we keep $K_{A0}$ = 0.90 for the July 2006 - June 2010 period, then
gradually decreases it to a minimum level for the polarity reversal phase from
November 2012 - March 2014. After the completed reversal, it is increased
gradually back and by mid-2016 it is again 0.90, and kept this way until the
end of the simulation period. It is emphasized that this type of variation is
needed to reproduce the proton flux shown in Figure 3. Furthermore, the value
of $P_{A0}$ = 0.75 GV from July 2006 - April 2010, but after November 2010 it
is 0.90 GV. This parameter determines the rigidity value below which $K_{D}$
is reduced in terms of rigidity. Adjusting $K_{A0}$ and $P_{A0}$ along with
${B_{m}}$ over time is the explicit way of drift adjustment; see also
Potgieter (2013, 2014), and Nndanganeni & Potgieter (2016, 2018) who also
discussed the implicit way. The global drift scale, $\lambda_{A}$ =
(${3}/{v})K_{D}$, in units of astronomical units (au), is found to vary from
July 2006 - November 2019 at 1.0 - 1.16 GV, which is selected as the lowest
rigidity observed by AMS. This is plotted in Figure 4(a). The times when
$K_{A0}$ and $P_{A0}$ are kept unchanged are shown so that the time trend in
$\lambda_{A}$ is caused by the changes in HMF strength as shown in Figure 1.
The drift scale as a function of rigidity at the Earth is shown in Figure 5 as
the lower set of graphs, also indicating how it changes with time (solar
activity). The rapid decrease in $\lambda_{A}$ below rigidity $\sim$1 GV is
the direct result of the scaling of $f_{D}$, essentially determined by
$P_{A0}$ in Equation 15; evidently, its value may also change with solar
activity.

Next, we focus on the expressions for the three diffusion coefficients, first,
$K_{\parallel}$ is assumed to be:

$K_{\parallel}=(K_{\parallel})_{0}\beta\Bigg{(}\frac{B_{u}}{B_{m}}\Bigg{)}\Bigg{(}\frac{P}{P_{u}}\Bigg{)}^{c_{1}}\left[\frac{\Bigg{(}\frac{P}{P_{u}}\Bigg{)}^{c_{3}}+\Bigg{(}\frac{P_{k}}{P_{u}}\Bigg{)}^{c_{3}}}{1+\Bigg{(}\frac{P_{k}}{P_{u}}\Bigg{)}^{c_{3}}}\right]^{\frac{c_{2\parallel}-c_{1}}{c_{3}}}$
(16)

where $(K_{\parallel})_{0}$ is a scaling constant in units of $10^{22}$
cm2s-1, with the rest of the equation written to be dimensionless with
$P_{u}=1.0$ GV, and $B_{u}$ = 1.0 nT (in order to preserve the units in cm2
s-1). Here $c_{1}$ is a power index that may change with time;
$c_{2\parallel}$ and $c_{2\perp}$ together with $c_{1}$ determine the slope of
the rigidity dependence, respectively, above and below a rigidity with the
value $P_{k}$ which may also change with time; $c_{3}$ determines the
smoothness of the transition. The rigidity dependence of $K_{\parallel}$ is
thus a combination of two power laws with $P_{k}$ determining the rigidity
where the transition occur, and the value of $c_{1}$ determines the slope of
the power law below $P_{k}$. Its rigidity dependence and how it changes with
time are shown in Figure 5 as the upper set of graphs.

We assume that $K_{\perp r}$ scales spatially similar to Equation (16) but
with a different rigidity dependence at higher rigidity by simply replacing
$c_{2\parallel}$ with $c_{2\perp r}$, which has a different value, to obtain
$\acute{(K)_{\parallel}}$; we then assume

$K_{\perp r}=0.02\acute{(K)_{\parallel}}$ (17)

which is widely used and reasonable assumption, e.g. Giacalone & Jokipii
(1999). On the other hand, $K_{\perp\theta}$, is more complicated, with
consensus that $K_{\perp\theta}$ $>$ $K_{\perp r}$ away from the equatorial
regions as discussed and motivated by Potgieter (2000) and Potgieter et al.
(2014). In this study, $K_{\perp\theta}$ is assumed to scale spatially similar
to Equation (16) but again with a different rigidity dependence at higher
rigidity obtained by replacing $c_{2\parallel}$ with ${c_{2\perp\theta}}$,
which also has a different value than ${c_{2\perp r}}$, to obtain
$\grave{(K)_{\parallel}}$. It is then assumed to be given by

$K_{\perp\theta}=0.02f_{\perp\theta}\grave{(K)_{\parallel}}$ (18)

with

$f_{\perp\theta}=A^{+}\mp A^{-}\tanh[8(\theta_{A}-90^{\circ})\pm\theta_{F}].$
(19)

Here, A${}^{\pm}=(d_{\perp\theta}\pm 1)/2$, $\theta_{F}$ = 35∘, $\theta_{A}$ =
$\theta$ for $\theta\leq 90^{\circ}$ but $\theta_{A}$ = 180∘ \- $\theta$ with
$\theta\geq$ 90∘ and $d_{\perp\theta}$ = 6.0 for protons and antiprotons. This
means that $K_{\perp\theta}$ has a different latitudinal dependence than the
other diffusion coefficients which can be enhanced towards the heliospheric
poles by a factor $d_{\perp\theta}$ with respect to the value of $K_{\perp r}$
in the equatorial region of the heliosphere. We keep $d_{\perp\theta}$ the
same for the whole simulation period. The rigidity dependence and how it
changes with time for these two MFPs are also shown in Figure 5 as the middle
group/band of curves. Note that in the equatorial plane $K_{\perp\theta}$ =
$K_{\perp r}$ but not elsewhere. For a detailed theoretical discussion of
these aspects of the diffusion tensor and theory, see the extensive reviews by
Shalchi (2009, 2020); Engelbrecht et al. (2022). For a motivation and various
applications of this particular modeling approach, see Potgieter (2000),
Ferreira et al. (2003), Potgieter & Langner (2005); Ngobeni & Potgieter
(2012), Nndanganeni & Potgieter (2016), to mention only a few.

Focusing on the time variation (July 2006 - November 2019) for a specified
rigidity (1.0 - 1.16 GV) of $\lambda_{\parallel}$, $\lambda_{\perp r}$ and
$\lambda_{\perp\theta}$ as plotted in Figure 4(b), it is noticed that they
follow the exact same time variation, at this rigidity. However, the interest
is on the trends over the solar cycle, and it is found that they have to
increase gradually from 2006 - 2009 in order to reproduce the observations of
this period. After December 2009, they keep reducing gradually as solar
activity increases to reach the lowest values in February - March 2014, but in
March - April 2015 a second short period of low values occur. This again is
related to the rapid changes displayed for this particular period in Figure 1.
From this time forward, they increase systematically towards the next solar
minimum.

Figure 5: Rigidity dependence of the three MFPs ($\lambda_{\parallel}$,
$\lambda_{\perp r}$, & $\lambda_{\perp\theta}$) and drift scale
($\lambda_{A}$) at Earth for both protons and antiprotons (see Equations 15 to
19). It also shows how they vary from July 2006 - November 2019 with values
for selected BRs indicated with different colors, e.g., July 2006 is when
PAMELA observations began; December 2009 is when the highest proton intensity
was observed; May 2011 is when the AMS observations began; October 2012 is end
of the $A<0$ polarity cycle; March 2014 is when the lowest proton intensity
was observed; May 2016 was when the drift scale recovery is completed.

###  5.2 Flux ratios over a solar cycle

The time variation of the antiproton to proton flux ratio at 1.0 - 1.16 GV for
the simulation period is shown in Figure 6; it is found to be of the order of
10-5. In the same panel, the proton to antiproton flux ratio is also shown.
Noteworthy trends from July 2006 - December 2009 are that the latter ratio
increased to a maximum value, after January 2010 it keeps increasing to reach
a minimum value at the end of the polarity reversal; after the reversal, it
increases rapidly within an year by $\sim$36% with respect to February - March
2014. The changes at the end of 2014 and beginning of 2015 show a large
fluctuation related to the changes in the flux. It then increases
continuously, and more gradual, to reach again a high value when solar minimum
sets in, with the exception of changing rapidly in the second half of 2017.
These fluctuations correspond to what is mentioned above and related to what
is shown in Figures 1 and 3.

Apart from showing flux ratios, another interesting and insightful way to
illustrate how differently the antiproton flux behaves with respect to the
proton flux over the simulated period, is to produce a plot as shown in Figure
7. The obtained loop is popularly called a hysteresis behaviour. The numbers
indicate 3-BR-averages, where ”1” corresponds to July - September 2006, then
e.g., ”16” corresponds to November 2009 - January 2010, ”29-35” to November
2012 - March 2014, which is the polarity reversal phase, and ”60” corresponds
to the September-November 2019, the end of the simulation period. Year numbers
are also indicated to guide the eye, showing how the looped-curve increases
from July 2006 to December 2009 with decreasing solar activity, decreases
along the same slope up to the end of 2011 with increasing solar activity.
Then in 2012 the shape (slope) begins to change, differently in the HMF
reversal phase. Even after the reversal, the shape remains scattered for
$\sim$1 year, then it begins to increase systematically up to the end of the
simulation period but with the slope being somewhat different from 2017
onward, and almost the same at in the previous solar minimum.

Figure 6: Variation of proton to antiproton flux ratio over time, from July
2006 - November 2019 (BRs 2361-2540) for 1.0 - 1.16 GV, shown in blue (left
side). The magnitude is in the order of $10^{5}$. For convenience, the
antiproton to proton flux ratio is shown in red (right side) with the
magnitude in the order of $10^{-5}$.  Figure 7:  Hysteresis-like loop between
antiproton and proton flux for 1.0-1.16 GV from July 2006 - November 2019 (BRs
2361-2540). Numbers (in blue) indicate 3-BR-averages, from ”1” corresponding
to BRs 2361-2363 and ends at ”60” corresponding to BRs 2538-2540.
Corresponding year numbers are marked in red. HMF reversal period is indicated
with a looped grey dashed line.

##  6 Discussions

Concerning solar activity, the recent solar minimum between Solar Cycles 24
and 25 seems comparable with the previous solar minimum between Solar Cycles
23 and 24. The HCS tilt angle $\alpha$ is below 10∘ with less fluctuations
since beginning of 2018, but the HMF magnitude at the Earth is higher than 4
nT and still having large fluctuations until the end of the simulation period.
The lowest value for $\alpha$ recorded during the previous minimum was 4.5∘ in
September-October 2009, with the field magnitude at 3.5 nT in March - April
2009. After this minimum both start to increase rapidly but with different
rates, with $\alpha$ reaching the highest value in January - February, 2013.
Afterwards it shows a very large fluctuation in 2014, and then progressively
decreases with modest fluctuations. The HMF rate on the other hand flattens
off during the polarity reversal period to reach its highest value of 7.9 nT
afterwards, in November - December 2014. Then it decreases progressively with
relative large fluctuations. According to http://wso.stanford.edu/, not shown
in Figure 1, $\alpha$ reaches a minimum of 2.1∘ by March - April 2020 before
starting to increase, with field magnitude reaching a minimum of 3.8 nT, also
in March-April 2020. It is noted that the Oulu Neutron Monitor recorded its
peak count rate in mid-2020, and is observably lower than the counts recorded
in the end-2009, see https://cosmicrays.oulu.fi/. This particular phenomenon
where the flux of GCRs at Neutron Monitor energies (high rigidity) is lower in
$A>0$ polarity cycles than in $A<0$ cycles, in contrast to what happens at low
rigidity as discussed next, is explained in detail by Krainev et al. (2021);
see also references therein.

Concerning the proton flux variations at 1.0 - 1.16 GV, it is noted that the
flux gradually increases since July 2006 with the highest flux during 07
December 2009-03 January 2010 (BR 2407/CR2091). It then decreases as solar
activity level increases and reaches the lowest flux level during 06 March -
02 April, 2014 (BR 2464), around the end of the polarity reversal phase.
Afterwards, the flux starts to recovery and by mid-2019 it reaches almost the
same level than in December 2009. Since protons drift inward mostly through
the heliospheric polar region in an $A>0$ cycle, it is expected that during
the recent minimum, somewhat higher fluxes than during the previous 2009
minimum should be observed at lower rigidity, given that all the heliospheric
conditions are similar; see predictions by Strauss & Potgieter (2014b);
Potgieter & Vos (2017); Krainev et al. (2021) and the recent report of proton
spectra observed between 50-250 MeV in 2020 to be higher than in 2009
(Martucci et al., 2023).

The best-known prediction of drift models with a changing HCS (Jokipii &
Thomas, 1981) is that the proton flux should exhibit a peaked intensity-time
profile during $A<0$ polarity cycles, but flatter profiles during $A>0$
cycles, also evident from Figure 3; see also e.g., Ferreira et al. (2003);
Engelbrecht & Moloto (2021). It implies that the antiproton intensity-time
profiles should be different than for protons, as is found and reported here
but it appears not as spectacularly large as what the original drift or drift
dominated modulation models predicted, e.g., by Jokipii & Kopriva (1979);
Potgieter & Moraal (1985); le Roux & Potgieter (1995); Ferreira & Potgieter
(2004); Tomassetti et al. (2017) and references there-in.

Concerning flux ratios, it is further noted that the $\sim$ 60% increase of
the proton to antiproton ratio by 2009 indicates that even when long periods
of low activity levels and quiet conditions occur, a significant rise in the
antiproton flux does not consequently follow. The simulated antiproton to
proton flux ratio in 2019 is lower than in 2009. The lowest proton to
antiproton flux ratio is found in February-March 2014, although the lowest
levels of the individual fluxes are almost a year apart. The trends in this
flux ratio is qualitatively similar to the positron to electron ratio but with
quantitative differences as simulated by Aslam et al. (2021, 2022); see also
Mikhailov et al. (2020). Concerning the loop-behaviour between the proton and
antiproton fluxes with time shown in Figure 7, it is emphasized that apart
from the different VLIS’s for these GCRs, the only difference between their
modulation is their oppositely directed drift directions under similar
modulation conditions. In this context, taken into consideration what is shown
in Figure 2 and together with the different drift patterns, the loop behaviour
reflects exactly how differently these GCR particles is modulated. The changes
in the slopes are significant from before to after the reversal period as the
new drift directions settle in. After the reversal the loop is not repeating
its ’path’ as before but follows a new one to end up being lower in 2019 than
in 2009 and with a somewhat different slope. This relates directly to the
discussion above about the different behaviour during the two polarity cycles
of oppositely charge particles.

An interesting trend in fluxes is evident during late-2017 when protons
display a large and relatively long decrease, interrupting the recovery to
solar minimum levels, and followed by a systematic recovery. Our simulated
antiproton flux also shows this transient type behaviour but with quantitative
differences which is considered as indicative of charge-sign dependence. This
aspect, specifically what causes this type of transient, needs further
investigation. It is probably an indication of the effect of a merged
interaction region (McDonald, 1998; Wibberenz et al., 1998) far beyond the
Earth; see simulations by e.g., le Roux & Fichtner (1999); Luo et al. (2019).

##  7 Summary and Conclusions

The main focus of this study is to improve our understanding of how
antiprotons are modulated over a solar cycle. Comparing the simulated proton
and antiproton spectra as well as the flux ratios give insight about charge-
sign dependence over such a long period. Our simulations are validated by
using observational proton spectra from PAMELA for July 2006 - May 2011 and
from AMS for May 2011-November 2019, thus covering the solar minimum period
before 2009, the period of maximum solar modulation including the time of the
reversal of the HMF polarity, and then up to solar minimum conditions in 2019.
It has turned out that the recent solar activity minimum around 2020 is
comparable with the minimum between Solar Cycles 23 and 24 in the sense that
the Sun was again rather quiet, also in terms of the tilt $\alpha$ of HCS, the
magnitude of the HMF and the GCR spectra that are observed. See the recent
reports of proton spectra observed at lower rigidity before and during 2020 by
Bartocci et al. (2020); Martucci et al. (2023).

The spectral shape and values of the VLIS’s for various GCRs are an important
aspect of modelling studies. In this context, the proton VLIS is well-
established, perhaps as good as it can be done using models compared to
observations beyond the HP and at the Earth at very high rigidity. The VLIS
for antiprotons is not well-established, although quite reasonably determined
by using GALPROP models and observations at high rigidity at the Earth. It is
shown here that its rather peculiar VLIS plays an important role in the total
amount of antiproton modulation between the HP and the Earth and in how its
modulation develops over a solar cycle with respect to that for protons.
However, comparing our modeling results with what was reported by Webber &
Potgieter (1989); Bieber et al. (1999), it is evident that we have indeed make
good progress, both in terms of what is assumed for the antiproton VLIS and
what is found for the flux ratios over the reversal period of the HMF
direction. Based on the newly determined VLIS for antiprotons, we found that
the total modulation obtained with the model between the HP (122 au) and Earth
(1 au) is vastly different than for protons. The peak intensity in the
antiproton VLIS occurs at a much higher rigidity than for protons,
subsequently also for the modulated spectra. The total modulation that
antiprotons experienced at low rigidity is markedly less than for protons.

A plot between antiproton and proton fluxes over the 2006 - 2019 period shows
a clear hysteresis behaviour indicative of how differently anti-protons are
modulated from protons over a solar cycle. Together with observations these
computational results clarify to a large extent the phenomenon of charge-sign
dependence in heliospheric modulation over a solar cycle. Our simulations show
that even for the quietest heliospheric conditions not significantly more
antiprotons, especially at lower rigidity, is reaching the Earth in sharp
contrast to protons.

Our results demonstrate how the underlying physics in the TPE varies over a
solar cycle displayed by the computed time-variations of the main diffusion
coefficients. In addition we found that the drift scale over time needs to be
adjusted significantly with solar activity, with its lowest level during the
phase of no well-defined HMF polarity. Furthermore, after the reversal period
its recovery to large values associated with low solar activity is a slow
process, it takes at least 2 years. This is consistent to what we reported
when modeling electron and positron spectra over a similar period although the
flux ratio of these particles is quantitatively different compared to protons
and antiprotons.

The authors thank Mirko Boezio, Riccardo Munini and Nadir Marcelli for
discussions of the PAMELA data and making the proton data set electronically
available. The authors also wish to thank the GALPROP developers and their
funding bodies for access to and use of the GALPROPWebRun service. We
acknowledge the use of HCS tilt angles data from the Wilcox Solar Observatory
(http://wso.stanford.edu), and HMF data from NASA’s OMNIWEB data interfaces
(http://omniweb.gsfc.nasa.gov). This work is supported by the Shandong
Institute of Advanced Technology (SDIAT) and NSFC Project U2106201. MDN
acknowledges the SA National Research Foundation (NRF) for partial financial
support under the Joint Science and Technology Research Collaboration between
SA and Russia (Grant no:118915).

## References

  * Adriani et al. (2013a) Adriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al. 2013a, ApJ, 765, 91, doi: 10.1088/0004-637X/765/2/91
  * Adriani et al. (2013b) Adriani, O., Bazilevskaya, G. A., Barbarino, G. C., et al. 2013b, JETP Letters, 96, 621, doi: 10.1134/S002136401222002X
  * Adriani et al. (2017) Adriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al. 2017, Nuovo Cimento Rivista Serie, 40, 473, doi: 10.1393/ncr/i2017-10140-x
  * Aguilar et al. (2016) Aguilar, M., Ali Cavasonza, L., Alpat, B., et al. 2016, Phys. Rev. Lett., 117, 091103, doi: 10.1103/PhysRevLett.117.091103
  * Aguilar et al. (2018) —. 2018, Phys. Rev. Lett., 121, 051101, doi: 10.1103/PhysRevLett.121.051101
  * Aguilar et al. (2021a) Aguilar, M., Ali Cavasonza, L., Ambrosi, G., et al. 2021a, Phys. Rep., 894, 1, doi: 10.1016/j.physrep.2020.09.003
  * Aguilar et al. (2021b) Aguilar, M., Cavasonza, L. A., Ambrosi, G., et al. 2021b, Phys. Rev. Lett., 127, 271102, doi: 10.1103/PhysRevLett.127.271102
  * Aslam et al. (2021) Aslam, O. P. M., Bisschoff, D., Ngobeni, M. D., et al. 2021, ApJ, 909, 215, doi: 10.3847/1538-4357/abdd35
  * Aslam et al. (2019) Aslam, O. P. M., Bisschoff, D., Potgieter, M. S., Boezio, M., & Munini, R. 2019, ApJ, 873, 70, doi: 10.3847/1538-4357/ab05e6
  * Aslam et al. (2022) Aslam, O. P. M., Luo, X., Potgieter, M. S., Ngobeni, M. D., & Song, X. 2022, arXiv e-prints, arXiv:2212.13397.  https://arxiv.org/abs/2212.13397
  * Bartocci et al. (2020) Bartocci, S., Battiston, R., Burger, W. J., et al. 2020, ApJ, 901, 8, doi: 10.3847/1538-4357/abad3e
  * Bieber et al. (1999) Bieber, J. W., Burger, R. A., Engel, R., et al. 1999, Phys. Rev. Lett., 83, 674, doi: 10.1103/PhysRevLett.83.674
  * Bisschoff et al. (2021) Bisschoff, D., Aslam, O. P. M., Ngobeni, M. D., et al. 2021, Phys. Atomic Nuc., 84, 1121, doi: 10.1134/S1063778821130068
  * Bisschoff & Potgieter (2016) Bisschoff, D., & Potgieter, M. S. 2016, Ap&SS, 361, 48, doi: 10.1007/s10509-015-2633-8
  * Bisschoff et al. (2019) Bisschoff, D., Potgieter, M. S., & Aslam, O. P. M. 2019, ApJ, 878, 59, doi: 10.3847/1538-4357/ab1e4a
  * Boezio et al. (2017) Boezio, M., Munini, R., Adriani, O., et al. 2017, in 35th Inter. Cosmic Ray Conf. (ICRC2017), Vol. 301, 1091, doi: 10.22323/1.301.01091
  * Boschini et al. (2019) Boschini, M. J., Della Torre, S., Gervasi, M., La Vacca, G., & Rancoita, P. G. 2019, Adv. Space Res., 64, 2459, doi: 10.1016/j.asr.2019.04.007
  * Boschini et al. (2020) Boschini, M. J., Della Torre, S., Gervasi, M., et al. 2020, ApJ, 889, 167, doi: 10.3847/1538-4357/ab64f1
  * Corti et al. (2016) Corti, C., Bindi, V., Consolandi, C., & Whitman, K. 2016, ApJ, 829, 8, doi: 10.3847/0004-637X/829/1/8
  * Corti et al. (2019) Corti, C., Potgieter, M. S., Bindi, V., et al. 2019, ApJ, 871, 253, doi: 10.3847/1538-4357/aafac4
  * Cummings et al. (2016) Cummings, A. C., Stone, E. C., Heikkila, B. C., et al. 2016, ApJ, 831, 18, doi: 10.3847/0004-637X/831/1/18
  * Di Felice et al. (2017) Di Felice, V., Munini, R., Vos, E. E., & Potgieter, M. S. 2017, ApJ, 834, 89, doi: 10.3847/1538-4357/834/1/89
  * Dunzlaff et al. (2015) Dunzlaff, P., Strauss, R. D., & Potgieter, M. S. 2015, Comp. Phys. Comm., 192, 156, doi: doi:10.1016/j.cpc.2015.03.008
  * Engelbrecht & Moloto (2021) Engelbrecht, N. E., & Moloto, K. D. 2021, ApJ, 908, 167, doi: 10.3847/1538-4357/abd3a5
  * Engelbrecht et al. (2022) Engelbrecht, N. E., Effenberger, F., Florinski, V., et al. 2022, Space Sci. Rev., 218, 33, doi: 10.1007/s11214-022-00896-1
  * Ferreira & Potgieter (2004) Ferreira, S. E. S., & Potgieter, M. S. 2004, Adv. Space Res., 34, 115, doi: 10.1016/j.asr.2003.02.057
  * Ferreira et al. (2003) Ferreira, S. E. S., Potgieter, M. S., Heber, B., & Fichtner, H. 2003, Ann. Geophys., 21, 1359, doi: 10.5194/angeo-21-1359-2003
  * Giacalone & Jokipii (1999) Giacalone, J., & Jokipii, J. R. 1999, ApJ, 520, 204, doi: 10.1086/307452
  * Hathaway (2015) Hathaway, D. H. 2015, Living Rev. Solar Phys., 12, 4, doi: 10.1007/lrsp-2015-4
  * Heber (2013) Heber, B. 2013, Space Sci. Rev., 176, 265, doi: 10.1007/s11214-011-9784-x
  * Heber et al. (2003) Heber, B., Clem, J. M., Müller-Mellin, R., et al. 2003, Geophys. Res. Lett., 30, 8032, doi: 10.1029/2003GL017356
  * Heber & Marsden (2001) Heber, B., & Marsden, R. G. 2001, Space Sci. Rev., 97, 309, doi: 10.1023/A:1011889319933
  * Heber & Potgieter (2006) Heber, B., & Potgieter, M. S. 2006, Space Sci. Rev., 127, 117, doi: 10.1007/s11214-006-9085-y
  * Heber et al. (1999) Heber, B., Ferrando, P., Raviart, A., et al. 1999, Geophys. Res. Lett., 26, 2133, doi: 10.1029/1999GL900482
  * Heber et al. (2002) Heber, B., Wibberenz, G., Potgieter, M. S., et al. 2002, J. Geophys. Res., 107, 1274, doi: 10.1029/2001JA000329
  * Jokipii & Kopriva (1979) Jokipii, J. R., & Kopriva, D. A. 1979, ApJ, 234, 384, doi: 10.1086/157506
  * Jokipii & Thomas (1981) Jokipii, J. R., & Thomas, B. 1981, ApJ, 243, 1115, doi: 10.1086/158675
  * Kopp et al. (2012) Kopp, A., Büsching, I., Strauss, R. D., & Potgieter, M. S. 2012, Comp. Phys. Comm., 183, 53, doi: doi:10.1016/j.cpc.2011.11.014
  * Kóta (2013) Kóta, J. 2013, Space Sci. Rev., 176, 391, doi: 10.1007/s11214-012-9870-8
  * Kóta & Jokipii (1983) Kóta, J., & Jokipii, J. R. 1983, ApJ, 265, 573, doi: 10.1086/160701
  * Krainev et al. (2021) Krainev, M., Kalinin, M., Aslam, O. P. M., Ngobeni, D., & Potgieter, M. S. 2021, Adv. Space Res., 68, 2953, doi: 10.1016/j.asr.2021.05.021
  * Langner & Potgieter (2004) Langner, U. W., & Potgieter, M. S. 2004, J. Geophys. Res., 109, A01103, doi: 10.1029/2003JA010158
  * Langner & Potgieter (2005) —. 2005, Adv. Space Res., 35, 2084, doi: 10.1016/j.asr.2004.12.005
  * Langner et al. (2006) Langner, U. W., Potgieter, M. S., Fichtner, H., & Borrmann, T. 2006, ApJ, 640, 1119, doi: 10.1086/500162
  * le Roux & Fichtner (1999) le Roux, J. A., & Fichtner, H. 1999, Adv. Space Res., 23, 501, doi: 10.1016/S0273-1177(99)00113-1
  * le Roux & Potgieter (1992) le Roux, J. A., & Potgieter, M. S. 1992, ApJ, 397, 686, doi: 10.1086/171824
  * le Roux & Potgieter (1995) —. 1995, ApJ, 442, 847, doi: 10.1086/175487
  * Luo et al. (2019) Luo, X., Potgieter, M. S., Bindi, V., Zhang, M., & Feng, X. 2019, ApJ, 878, 6, doi: 10.3847/1538-4357/ab1b2a
  * Luo et al. (2020) Luo, X., Zhang, M., Feng, X., et al. 2020, ApJ, 899, 90, doi: 10.3847/1538-4357/aba7b5
  * Luo et al. (2015) Luo, X., Zhang, M., Potgieter, M., Feng, X., & Pogorelov, N. V. 2015, ApJ, 808, 82, doi: 10.1088/0004-637X/808/1/82
  * Manuel et al. (2014) Manuel, R., Ferreira, S. E. S., & Potgieter, M. S. 2014, Sol. Phys., 289, 2207, doi: 10.1007/s11207-013-0445-y
  * Martucci et al. (2023) Martucci, M., Bartocci, S., Battiston, R., & et al. 2023, ApJ, in press.
  * Martucci et al. (2018) Martucci, M., Munini, R., Boezio, M., et al. 2018, ApJ, 854, L2, doi: 10.3847/2041-8213/aaa9b2
  * McComas et al. (2002) McComas, D. J., Elliott, H. A., Gosling, J. T., et al. 2002, Geophys. Res. Lett., 29, 1290, doi: 10.1029/2001GL014164
  * McDonald (1998) McDonald, F. B. 1998, Space Sci. Rev., 83, 33, doi: 10.1023/A:1005052908493
  * Mikhailov et al. (2020) Mikhailov, V. V., Aleksandrin, S. Y., Koldobskiy, S. A., et al. 2020, Journal of Physics: Conference Series, 1690, 012004, doi: 10.1088/1742-6596/1690/1/012004
  * Moloto et al. (2018) Moloto, K. D., Engelbrecht, N. E., & Burger, R. A. 2018, ApJ, 859, 107, doi: 10.3847/1538-4357/aac174
  * Moraal & Potgieter (1982) Moraal, H., & Potgieter, M. S. 1982, Ap&SS, 84, 519, doi: 10.1007/BF00651330
  * Moskalenko et al. (2005) Moskalenko, I. V., Strong, A. W., Ormes, J. F., & Mashnik, S. G. 2005, Adv. Space Res., 35, 156, doi: 10.1016/j.asr.2003.08.050
  * Moskalenko et al. (2002) Moskalenko, I. V., Strong, A. W., Ormes, J. F., & Potgieter, M. S. 2002, ApJ, 565, 280, doi: 10.1086/324402
  * Munini et al. (2017) Munini, R., Boezio, M., Potgieter, M., et al. 2017, in 35th International Cosmic Ray Conference (ICRC2017), Vol. 301, 12, doi: 10.22323/1.301.0012
  * Ngobeni & Potgieter (2012) Ngobeni, M. D., & Potgieter, M. S. 2012, Adv. Space Res., 49, 1660, doi: 10.1016/j.asr.2012.02.036
  * Ngobeni & Potgieter (2015) —. 2015, Adv. Space Res., 56, 1525, doi: 10.1016/j.asr.2015.06.034
  * Ngobeni et al. (2022) Ngobeni, M. D., Potgieter, M. S., Aslam, O. P. M., et al. 2022, Adv. Space Res., 69, 2330, doi: 10.1016/j.asr.2021.12.018
  * Ngobeni et al. (2020) Ngobeni, M. D., Aslam, O. P. M., Bisschoff, D., et al. 2020, Ap&SS, 365, 182, doi: 10.1007/s10509-020-03896-1
  * Nndanganeni & Potgieter (2016) Nndanganeni, R. R., & Potgieter, M. S. 2016, Adv. Space Res., 58, 453, doi: 10.1016/j.asr.2016.04.020
  * Nndanganeni & Potgieter (2018) —. 2018, Ap&SS, 363, 156, doi: 10.1007/s10509-018-3377-z
  * Parker (1958) Parker, E. N. 1958, ApJ, 128, 664, doi: 10.1086/146579
  * Parker (1965) —. 1965, Planet. Space Sci., 13, 9, doi: 10.1016/0032-0633(65)90131-5
  * Potgieter (2000) Potgieter, M. S. 2000, J. Geophys. Res., 105, 18295, doi: 10.1029/1999JA000434
  * Potgieter (2013) —. 2013, Living Rev. Solar Phys., 10, 3, doi: 10.12942/lrsp-2013-3
  * Potgieter (2014) —. 2014, Adv. Space Res., 53, 1415, doi: 10.1016/j.asr.2013.04.015
  * Potgieter (2017) —. 2017, Adv. Space Res., 60, 848, doi: 10.1016/j.asr.2016.09.003
  * Potgieter et al. (2021) Potgieter, M. S., Aslam, O. P. M., Bisschoff, D., & Ngobeni, M. D. 2021, Physics, 3, 1190, doi: 10.3390/physics3040076
  * Potgieter & Langner (2005) Potgieter, M. S., & Langner, U. W. 2005, Adv. Space Res., 35, 554, doi: 10.1016/j.asr.2005.01.008
  * Potgieter & Moraal (1985) Potgieter, M. S., & Moraal, H. 1985, ApJ, 294, 425, doi: 10.1086/163309
  * Potgieter & Vos (2017) Potgieter, M. S., & Vos, E. E. 2017, A&A, 601, A23, doi: 10.1051/0004-6361/201629995
  * Potgieter et al. (2014) Potgieter, M. S., Vos, E. E., Boezio, M., et al. 2014, Sol. Phys., 289, 391, doi: 10.1007/s11207-013-0324-6
  * Potgieter et al. (2015) Potgieter, M. S., Vos, E. E., Munini, R., Boezio, M., & Di Felice, V. 2015, ApJ, 810, 141, doi: 10.1088/0004-637X/810/2/141
  * Raath et al. (2016) Raath, J. L., Potgieter, M. S., Strauss, R. D., & Kopp, A. 2016, Adv. Space Res., 57, 1965, doi: 10.1016/j.asr.2016.01.017
  * Raath et al. (2015) Raath, J. L., Strauss, R. D., & Potgieter, M. S. 2015, Ap&SS, 360, 24, doi: 10.1007/s10509-015-2556-4
  * Rankin et al. (2022) Rankin, J. S., Bindi, V., Bykov, A. M., et al. 2022, Space Sci. Rev., 218, 42, doi: 10.1007/s11214-022-00912-4
  * Richardson & Wang (2011) Richardson, J. D., & Wang, C. 2011, ApJ, 734, L21, doi: 10.1088/2041-8205/734/1/L21
  * Scherer et al. (2011) Scherer, K., Fichtner, H., Strauss, R. D., et al. 2011, ApJ, 735, 128, doi: 10.1088/0004-637X/735/2/128
  * Shalchi (2009) Shalchi, A. 2009, Nonlinear Cosmic Ray Diffusion Theories, Vol. 362 (Springer: Berlin, Heidelberg), doi: 10.1007/978-3-642-00309-7
  * Shalchi (2020) —. 2020, Space Sci. Rev., 216, 23, doi: 10.1007/s11214-020-0644-4
  * Shen et al. (2021) Shen, Z., Qin, G., Zuo, P., Wei, F., & Xu, X. 2021, ApJS, 256, 18, doi: 10.3847/1538-4365/ac0a78
  * Smith & Bieber (1991) Smith, C. W., & Bieber, J. W. 1991, ApJ, 370, 435, doi: 10.1086/169830
  * Song et al. (2021) Song, X., Luo, X., Potgieter, M. S., Liu, X., & Geng, Z. 2021, ApJS, 257, 48, doi: 10.3847/1538-4365/ac281c
  * Stone et al. (2019) Stone, E. C., Cummings, A. C., Heikkila, B. C., & Lal, N. 2019, Nature Astron., 3, 1013, doi: 10.1038/s41550-019-0928-3
  * Stone et al. (2013) Stone, E. C., Cummings, A. C., McDonald, F. B., et al. 2013, Science, 341, 150, doi: 10.1126/science.1236408
  * Strauss & Potgieter (2014a) Strauss, R. D., & Potgieter, M. S. 2014a, Adv. Space Res., 53, 1015, doi: 10.1016/j.asr.2014.01.004
  * Strauss & Potgieter (2014b) —. 2014b, Sol. Phys., 289, 3197, doi: 10.1007/s11207-014-0517-7
  * Strauss et al. (2011) Strauss, R. D., Potgieter, M. S., Kopp, A., & Büsching, I. 2011, J. Geophys. Res., 116, A12105, doi: 10.1029/2011JA016831
  * Strong & Moskalenko (1998) Strong, A. W., & Moskalenko, I. V. 1998, ApJ, 509, 212, doi: 10.1086/306470
  * Sun et al. (2015) Sun, X., Hoeksema, J. T., Liu, Y., & Zhao, J. 2015, ApJ, 798, 114, doi: 10.1088/0004-637X/798/2/114
  * Tomassetti et al. (2017) Tomassetti, N., Orcinha, M., Barão, F., & Bertucci, B. 2017, ApJ, 849, L32, doi: 10.3847/2041-8213/aa9373
  * Usoskin (2013) Usoskin, I. G. 2013, Living Rev. Solar Phys., 10, 1, doi: 10.12942/lrsp-2013-1
  * Vos & Potgieter (2015) Vos, E. E., & Potgieter, M. S. 2015, ApJ, 815, 119, doi: 10.1088/0004-637X/815/2/119
  * Vos & Potgieter (2016) —. 2016, Sol. Phys., 291, 2181, doi: 10.1007/s11207-016-0945-7
  * Webber & Potgieter (1989) Webber, W. R., & Potgieter, M. S. 1989, ApJ, 344, 779, doi: 10.1086/167842
  * Wibberenz et al. (1998) Wibberenz, G., Le Roux, J. A., Potgieter, M. S., & Bieber, J. W. 1998, Space Sci. Rev., 83, 309, doi: 10.1023/A:1005083109827
  * Zhao et al. (2018) Zhao, L. L., Adhikari, L., Zank, G. P., Hu, Q., & Feng, X. S. 2018, ApJ, 856, 94, doi: 10.3847/1538-4357/aab362