Modeling of Heliospheric Modulation of Cosmic-Ray Positrons in a Very Quiet Heliosphere

The heliospheric conditions were quite unusual during the solar minimum period 2006–2009, between solar cycles 23 and 24. The minimum was unusually long and deep, the heliospheric magnetic field (HMF) was much weaker, and the solar wind was slowest and least dense compared to previous solar cycles. However, the tilt angle (α) of the heliospheric current sheet (HCS) was more warped (less flat) and had not decreased rapidly, as compared to previous minima, but reached a minimum value at the end of 2009 (Aslam & Badruddin 2012; Potgieter et al. 2015). The very low and almost constant intensity of the solar polar field (about half of the previous two minima) during this unusual minimum activity period was also of interest from the point of view of cosmic-ray (CR) modulation (e.g., Wang et al. 2009; Jian et al. 2011). A record high level of CR intensity since the beginning of the era of neutron-monitor (NM) observations was observed by these ground-based detectors, and the highest-ever galactic proton spectrum at Earth was observed by PAMELA during the end of 2009; see the discussion by Heber et al. (2009), Mewaldt et al. (2010), Potgieter et al. (2014), Aslam & Badruddin (2015), and other relevant references therein. This unusual increase in CR intensity also indicates that this minimum had an extraordinary effect on the properties of the magnetic structure shielding the Earth to allow such an increase in CRs (White et al. 2011).

Through the level of solar activity, velocity of the solar wind, tilt angle of the HCS, and magnitude of the HMF, the Sun controls the heliospheric structure and the modulation of CRs (as an example, see McDonald et al. 2010; Vos & Potgieter 2015, 2016). The Sun's magnetic field polarity reverses during every solar maximum activity phase; if the field lines are directed outward from the Sun's northern hemisphere, it is defined as positive polarity (A > 0), and if the field lines are outward from the Sun's southern hemisphere, it is defined as negative polarity (A < 0). The HCS separates these two oppositely directed magnetic polarity hemispheres. As α varies with solar activity, the changing HCS has a significant effect on CR modulation through the drift motion of CR particles (e.g., Jokipii & Thomas 1981; Potgieter & Moraal 1985; Ferreira et al. 2003; Zhao et al. 2014). The period from 2001 to 2012 May (including the unusual 2006–2009 minimum) had a negative magnetic polarity (A < 0; see http://wso.stanford.edu).

The CRs are subjected to four distinct transport processes in the solar wind plasma flow with its imbedded HMF: (a) convection because of the outward-directed solar wind velocity; (b) gradient, curvature, and current sheet drifts; (c) adiabatic energy changes depending on the sign of the divergence of the solar wind velocity; and (d) spatial diffusion caused by the scattering of random magnetic irregularities (for reviews, see, e.g., Heber 2013; Kota 2013; Potgieter 2013). The drift process has different effects in each solar activity cycle; however, the diffusion and convection processes are dependent on the solar activity cycle rather than the solar magnetic polarity.

According to drift models, the positively charged CR particles (protons, positrons, helium, etc.) drift inward mainly along the heliospheric equatorial regions and outward via the polar regions of the heliosphere during the A < 0 polarity phase. On the other hand, negatively charged CR particles (electrons, antiprotons) will drift downward from the poles and outward through the equatorial regions during this polarity phase. The particle drift direction reverses for both positively and negatively charged CR particles during the opposite polarity configuration (A > 0; Kota & Jokipii 1983). This means that when Galactic positrons drift inward through the equatorial regions, the wavy HCS plays a prominent role in their modulation during the A < 0 polarity phase, as happened during the minimum of 2006–2009.

The detailed study and illustration of the response of Galactic protons during this unusual solar minimum period using a comprehensive three-dimensional (3D) heliospheric modulation model including gradient, curvature, and HCS drift applied to the PAMELA proton observations (Adriani et al. 2013a) was carried out by Potgieter et al. (2014) and Vos & Potgieter (2015). Similar detailed heliospheric modulation of Galactic electrons was also studied by Potgieter et al. (2015) and Potgieter & Vos (2017) by utilizing the PAMELA electron observations (Adriani et al. 2015; Munini et al. 2015) over this period. These authors concluded that all four modulation processes played important roles, but diffusion was relatively dominant during this unusual minimum, in contrast to the paradigm for the modulation of CRs that drift should be dominating during solar minima, but that diffusion is considered to become the dominant process during solar maximum. The current study is a continuation of these studies, now focusing on the Galactic positron observations made by PAMELA from mid-2006 to the end of 2009 (Adriani et al. 2013b; Munini 2015; Munini et al. 2017). Such a comprehensive study of positrons has not been done before.

The PAMELA detector is a spaceborne particle spectrometer, specially designed to measure the spectra of the primary and secondary components of CRs at Earth (Adriani et al. 2014, 2017; Boezio et al. 2017). The main focus was the indirect study of dark matter through the detection of the antiproton and positron spectra up to 200 GeV. The Carrington rotation–averaged proton spectra from PAMELA for the period 2006–2009 over an energy range of ≈80 MeV to 50 GeV were presented by Adriani et al. (2013a). The PAMELA electron spectra for the same period over the energy range ≈80 MeV to 40 GeV were 6 month–averaged, instead of Carrington rotation–averaged, and also reported by Adriani et al. (2015). The first positron spectra above 1 GeV from PAMELA were published by Adriani et al. (2013b). Later, the 6 month–averaged positron spectra for the period 2006 July to the end of 2009 down to 80 MeV were reported by Munini (2015), with the abundance of the positron intensity roughly 10% of the electron intensity and around 1% of the proton intensity at 10 GeV.

Figure 1 shows the positron energy spectra measured by PAMELA, averaged over 6 month periods, from 2006 July to 2009 December (adapted from Munini 2015). The total period is divided into seven semesters of 6 months and labeled "a," January–June, and "b," July–December (i.e., 2006b, the second semester of the year, represents the average over 2006 July–December, and 2007a, the first semester of the year, represents the average over 2007 January–June). The positron spectra presented here have an energy range from ≈80 MeV to 30 GeV. It should be noted that the flux data below 200 MeV are affected by significant systematic uncertainties. Evidently, the highest intensity was observed during the 2009b period.

Zoom In Zoom Out Reset image size

Figure 2 shows the positron intensity ratio relative to that for 2006b as a function of kinetic energy in GeV. Note that the spectra up to 200 MeV are the most responsive to changes in modulation conditions and have undergone an increase of a factor of ≈6. This increase comes down to a factor of ≈2 at 1 GeV and become less pronounced above 10 GeV as the solar modulation of CRs subsides; above ≈30 GeV, the ratio shows very little change (see also Strauss & Potgieter 2014).

Zoom In Zoom Out Reset image size

The prime objective of this work is to reproduce the PAMELA positron spectra observed during the period 2006 July–2009 December using the comprehensive 3D modulation model described below.

As an initial condition, a Galactic positron spectrum—more specifically, a local interstellar spectrum (LIS)—has to be specified in the model to be used as the input spectrum that is then modulated from a given heliospheric boundary up to the Earth at 1 au.

Figure 3 illustrates the differences (variations) between the LISs of protons, electrons, and positrons computed with the web version of the GALPROP code (see, e.g., Moskalenko & Strong 1998) and the corresponding LISs from previous modulation studies. The proton LIS was constructed by Vos & Potgieter (2015) and is similar to the GALPROP proton LIS; the electron LIS constructed by Potgieter et al. (2015) and Potgieter & Vos (2017) and the positron LIS constructed by Potgieter et al. (2017), however, show deviations as a function of kinetic energy from the LIS produced by the GALPROP code. The latter authors used PAMELA and later also AMS-02 observations at the Earth at high energies where modulation is negligible, and particularly Voyager 1 observations (Gurnett et al. 2013; Stone et al. 2013; Webber & McDonald 2013) at very low energies beyond the heliopause to construct the presented LIS of protons and electrons; and PAMELA 2006–2009, also AMS-02 2011–2013 observations to construct the positron LIS.

Zoom In Zoom Out Reset image size

For our modeling, we first revisited electron modulation and reproduced the PAMELA electron spectra for seven semesters from 2006b to 2009b, similar to what was done by Potgieter et al. (2015). For this study, we used the electron LIS computed with the GALPROP code. We use this approach to verify our modulation approach and constrain all the modulation parameters in order to repeat the same process for positrons. Arguing that the only two differences between electron and positron modulation are their respective LISs and the particle drift that they experience, we set out to reproduce the PAMELA positron spectra at the Earth using the LIS for positrons as computed with GALPROP. However, it soon became evident that we could not reproduce the modulated PAMELA positron spectra when using this LIS, despite various tuning of the GALPROP code, so we had to empirically construct a different LIS but still based on the GALPROP LIS. In Figure 4, we illustrate the differences between our LIS and the one from GALPROP (black lines) as a function of kinetic energy and the corresponding differences between the subsequent modulated spectra for the periods 2006b (blue lines) and 2009b (red lines) only. Evidently, the differences are meaningful. This modified LIS for positrons is used in the rest of our study and could reproduce the observed positron spectra, as will be shown and discussed below. The tuning of the GALPROP code as mentioned above is described in detail by Bisschoff (2017); see also Bisschoff & Potgieter (2014, 2016) and Bisschoff et al. (2019).

Zoom In Zoom Out Reset image size

A full 3D model based on the numerical solution of Parker's heliospheric transport equation (TPE; Parker 1965) is used to compute the differential intensity of CR positrons over an energy range from 1 MeV to 70 GeV at the Earth (1 au) and at different radial distances up to 122 au, the position of heliopause in the model where the LIS is specified as an initial condition. No modulation is considered beyond the heliopause; for such an approach, see Luo et al. (2015).

Parker's transport equation is described as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial f}{\partial t} & = & -{{\boldsymbol{V}}}_{\mathrm{sw}}\cdot {\rm{\nabla }}{\text{}}f-\langle {{\boldsymbol{v}}}_{D}\rangle \cdot {\rm{\nabla }}{\text{}}f+{\rm{\nabla }}\cdot ({{\boldsymbol{K}}}_{s}\cdot {\rm{\nabla }}{\text{}}f)\\ & & +\displaystyle \frac{1}{3}({\rm{\nabla }}\cdot {{\boldsymbol{V}}}_{\mathrm{sw}})\displaystyle \frac{\partial f}{\partial \mathrm{ln}P},\end{array}\end{eqnarray} \tag{ 1 }$$

where $f({\boldsymbol{r}},P,t)$ is the CR distribution function, P is rigidity, t is time, and ${\boldsymbol{r}}$ is the vector position in 3D, with the three coordinates r, θ, and ϕ specified in a heliocentric spherical coordinate system where the equatorial plane is at a polar angle of θ = 90°. The four terms shown on the right-hand side of Equation (1) represent the four major physical processes that CR particles undergo when they enter and travel through the heliosphere up to the Earth.

The first term represents the outward convection caused by the expanding solar wind with velocity (${{\boldsymbol{V}}}_{\mathrm{sw}}$). The second term represents the averaged particle drift velocity $\langle {{\boldsymbol{v}}}_{D}\rangle$ (pitch angle–averaged guiding center drift velocity), which is described by

$$\begin{eqnarray}&&\langle {{\boldsymbol{v}}}_{D}\rangle ={\rm{\nabla }}\times {K}_{D}\displaystyle \frac{{\boldsymbol{B}}}{B},\end{eqnarray} \tag{ 2 }$$

where K_D is the generalized drift coefficient, and ${\boldsymbol{B}}$ is the HMF vector with magnitude B. Both the geometry and magnitude of the HMF are important for the modulation process, particularly in the case of gradient, curvature, and HCS drift. As a departure point, a straightforward Parker HMF (Parker 1958) is assumed,

$$\begin{eqnarray}&&{\boldsymbol{B}}={B}_{0}A{\left[\displaystyle \frac{{r}_{0}}{r}\right]}^{2}({\hat{e}}_{r}-\tan \psi {\hat{e}}_{\phi })[1-2H(\theta -{\theta }^{{\prime} })],\end{eqnarray} \tag{ 3 }$$

where ${\hat{e}}_{r}$ and ${\hat{e}}_{\phi }$ are unit vectors in the radial and azimuthal directions, A = ±1, expressing the polarity phase of the Sun (A = +1 is A > 0, positive polarity, and A = −1 is A < 0, negative polarity); B₀ is the magnitude of the HMF at Earth (i.e., r₀ = 1 au); and the spiral angle ψ is the angle between the radial direction and the HMF field line at any given position, defined by

$$\begin{eqnarray}&&\tan \psi ={\rm{\Omega }}\displaystyle \frac{(r-{r}_{\odot })}{{V}_{\mathrm{sw}}}\sin \theta ,\end{eqnarray} \tag{ 4 }$$

where r_⊙ is the solar radius (0.005 au) and Ω is the average angular rotation speed of the Sun (2.66 × 10⁻⁶ rad s⁻¹). The magnitude of this Parker HMF is

$$\begin{eqnarray}&&B={B}_{0}{\left[\displaystyle \frac{{r}_{0}}{r}\right]}^{2}\sqrt{1+{(\tan \psi )}^{2}.}\end{eqnarray} \tag{ 5 }$$

Smith & Bieber (1991) suggested a modification to the Parker field based on their observation 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 that would lead to larger spiral angles at larger radial distances. They proposed a modification so that the expression for the HMF spiral angle (Equation (4)) becomes

$$\begin{eqnarray}&&\tan \psi =\displaystyle \frac{{\rm{\Omega }}(r-{r}_{b}){\rm{\sin }}\theta }{{V}_{\mathrm{sw}}(r,\theta )}-\displaystyle \frac{{{rV}}_{\mathrm{sw}}({r}_{b},\theta )}{{r}_{b}{V}_{\mathrm{sw}}(r,\theta )}\left(\displaystyle \frac{{B}_{T}({r}_{b})}{{B}_{R}({r}_{b})}\right),\end{eqnarray} \tag{ 6 }$$

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. Smith & Bieber (1991) showed that the ratio B_T(r_b)/B_R(r_b) ≈ −0.02 at a position r_b = 20 r_⊙. With r_⊙ = 0.005 au as the solar radius, the value r_b = 20 r_⊙ and the ratio B_T(r_b)/B_R(r_b) = −0.02 are constants that determine the HMF modification. This modification has a significant effect on the HMF structure at high-latitude regions; it keeps the basic Parkerian geometry but modifies its magnitude progressively toward the poles of the heliosphere. The motivation for this modification from a modulation modeling point of view was discussed in detail by Potgieter (2013) and Potgieter et al. (2015), while Raath et al. (2016) gave an elaborate discussion of its relevance and CR modulation effects.

For a modified Parker-type HMF (${{\boldsymbol{B}}}_{m}$), such as the Smith-Bieber modification, with magnitude (B_m), the expression for a generalized drift coefficient is

$$\begin{eqnarray}&&{K}_{D}=\displaystyle \frac{\beta P}{3{B}_{m}}{f}_{D}=\displaystyle \frac{\beta P}{3{B}_{m}}\left[\displaystyle \frac{{(\omega \tau )}^{2}}{1+{(\omega \tau )}^{2}}\right],\end{eqnarray} \tag{ 7 }$$

where β = v/c is the ratio of particle speed to the speed of light, ω is the particle gyrofrequency, and τ is the average time between the scattering of CR particles in the HMF. In most of the numerical modeling studies, it is assumed that ωτ ≫ 1 in the heliosphere so that the drift coefficient takes its simplest form,

$$\begin{eqnarray}&&{K}_{D}=\displaystyle \frac{\beta P}{3{B}_{m}},\end{eqnarray} \tag{ 8 }$$

which is known as weak scattering drift. Establishing τ is complicated and controversial; an elaborate turbulence theory is required to understand how ωτ should change as a function of rigidity and space throughout the entire heliosphere (see also the discussion by Ngobeni & Potgieter 2015). 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 $\langle {v}_{D}\rangle$ become zero, meaning that drift effects will vanish from the modulation model; if f_D = 1, drift is at a maximum, so that K_D has the weak scattering value. Then, the drift effects on CR modulation are very large and dominant, as originally applied in numerical models by Jokipii & Thomas (1981) and Kota & Jokipii (1983).

This function can be used to adjust the rigidity dependence of K_D, which is the most effective direct way of suppressing drift effects at low rigidities, so that Equation (7) becomes

$$\begin{eqnarray}&&{K}_{D}=\displaystyle \frac{\beta P}{3{B}_{m}}{f}_{D}={K}_{A0}\displaystyle \frac{\beta P}{3{B}_{m}}\displaystyle \frac{{(P/{P}_{A0})}^{2}}{1+{(P/{P}_{A0})}^{2}}.\end{eqnarray} \tag{ 9 }$$

Here K_A0 is a dimensionless constant that could be ranging from 0 to 1.0; if K_A0 = 1.0, it is called 100% drift (full weak scattering). In this study, we keep K_A0 = 0.90 and P_A0 = 0.90 GV for all seven semesters (mid-2006 to the end of 2009), which means that particle drift is at a 90% level during this unusual solar minimum period, but below ≈1.0 GV, the drift coefficient is reduced with respect to the weak scattering approach in Equation (8). For more details, see Potgieter (2013, 2014) and Nndanganeni & Potgieter (2016, 2018).

The drift scale (λ_A) with respect to rigidity at the Earth (1 au with polar angle θ = 90°) is shown in Figure 6. The fast decrease in λ_A below rigidity ≈1 GV is the direct result of the scaling of f_D. The weak scattering drift (Equation (8)) and modified drift coefficient (Equation (9)) have the same rigidity dependence above ≈1 GV. The reason this is required in numerical modeling was discussed in detail in the review by Heber & Potgieter (2006); see also Potgieter (2014). However, in the case of low-energy positrons and electrons, particle drift becomes negligible because their transport in the heliosphere is dominated by the diffusion process; see the review by Potgieter (2017). In what follows, these modulation effects will be shown explicitly for positrons. Unlike CR protons and heavy nuclei, the electrons and positrons experience far less adiabatic energy losses at low energies, so that they respond directly to changes of the diffusion coefficients down to very low kinetic energy. Previous numerical studies of Galactic protons and electrons by Potgieter et al. (2015) and Potgieter & Vos (2017) used K_A0 = 1.0 and P_A0 = 0.55 GV for the 2006–2009 minimum period.

Reducing drifts using B and f_D is a way of explicit drift reduction, whereas implicit drift reduction is a way of reducing drift effects without changing K_D. This is made possible by reducing the diffusive process by increasing any of the appropriate diffusion coefficients of the diffusion tensor K_s. This was illustrated in detail for CR electrons by Nndanganeni & Potgieter (2016).

The third term on the right-hand side of Equation (1) indicates the spatial diffusion caused by the scattering of CRs, where K_s is the symmetry diffusion tensor, and the last term represents the adiabatic energy change, which depends on the sign of the divergence of the ${{\boldsymbol{V}}}_{\mathrm{sw}}$. If $({\rm{\nabla }}\cdot {{\boldsymbol{V}}}_{\mathrm{sw}})\gt 0$, adiabatic energy losses occur, as is the case in most of the heliosphere, except inside the heliosheath, where we assume that $({\rm{\nabla }}\cdot {{\boldsymbol{V}}}_{\mathrm{sw}})=0$. Adiabatic energy loss is one of the major mechanisms and is important for CR modulation in most of the heliosphere but far less dominant in the case of electrons and positrons than for protons and other heavy CR species; see also the illustrations by, e.g., Moraal & Potgieter (1982).

When we write the TPE given by Equation (1) in the heliocentric spherical coordinate system, there are nine elements of the diffusion tensor K_s−K_rr, K_rθ, K_rϕ, K_θr, K_θθ, K_θϕ, K_ϕr, K_ϕθ, and K_ϕϕ−based on a Parker-type HMF with a radial solar wind speed V_sw. Here K_rr, K_rϕ, K_θθ, K_ϕr, and K_ϕϕ describe the diffusion process and K_rθ, K_θr, K_θϕ, and K_ϕθ describe the gradient, curvature, and HCS drift. These diffusion coefficients can be expressed in a 3D heliocentric spherical coordinate system as follows:

$$\begin{eqnarray}&&{K}_{{rr}}={K}_{\parallel }{\cos }^{2}\psi +{K}_{\perp r}{\sin }^{2}\psi ,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{K}_{\theta \theta }={K}_{\perp \theta },\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{K}_{\phi \phi }={K}_{\perp r}{\cos }^{2}\psi +{K}_{\parallel }{\sin }^{2}\psi ,\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{K}_{\phi r}=({K}_{\perp r}-{K}_{\parallel })\cos \psi \sin \psi ={K}_{r\phi }.\end{eqnarray} \tag{ 13 }$$

The effective radial diffusion coefficient K_rr is a combination of the parallel diffusion coefficient (K_∥) and the radial perpendicular diffusion coefficient (K_⊥r), with ψ the spiral angle; K_θθ = K_⊥θ is the effective perpendicular diffusion coefficient in the polar direction. Here K_ϕϕ and K_ϕr describe the effective diffusion in the azimuthal direction and diffusion coefficient in the ϕr plane, respectively.

Potgieter et al. (2014, 2015) described the 3D model used for this study in detail, with recent comprehensive reviews of the underlying theory for the global modulation of CRs during a quiet heliosphere given by Potgieter (2013, 2014, 2017).

4.1. Calculation and Selection of Modulation Parameters

The 3D numerical model for this modulation study is for a steady state, so that for the left side of Equation (1), ∂f/∂t = 0; this means that transient events such as Forbush decreases cannot be studied with this model. Determination of the averaged values for time-varying modulation parameters is an important part of the modeling method when working with such a model; e.g., the tilt angle α of the HCS and magnitude of the HMF at Earth B_e are varying continuously according to the solar activity. The time variations in α and B_e are included in the numerical model to set up realistic modulation conditions, which is reasonable for the relatively quiet state of the heliosphere. For a similar approach applied to AMS-02 data, see Corti et al. (2019).

The average solar wind speed is assumed to be ≈430 km s⁻¹ in the slow solar wind region and ≈750 km s⁻¹ in the fast solar wind region. During minimum to moderate solar activity periods, the HCS is mostly confined to the ecliptic region (within ≈30°), remaining in the slow solar wind region. The slow solar wind velocity is thus considered to calculate the time it takes changes in α to travel from the Sun to the heliopause. However, the HMF is not confined to the elliptic region, so the weighted average of both slow and fast solar wind velocities is used to calculate the time taken to reach the heliopause from the Sun. The propagation time calculated for α is ≈16 months (17 solar rotations), and for the HMF, it is ≈10 months (11 solar rotations). The moving averages over this propagation time for both α and B_e are used as intrinsic parameters in this modeling approach, similar to what Potgieter et al. (2015) and Vos & Potgieter (2015) assumed. The observed α and B_e are shown in Figure 5, together with this averaging process, and are summarized in Table 1.

Zoom In Zoom Out Reset image size

Table 1.  Summary of the Calculated Intrinsic and Modulation Parameters at 1 GV as Used to Reproduce the Positron Measurements from PAMELA for the Period 2006 July to 2009 December

Parameters	2006b	2007a	2007b	2008a	2008b	2009a	2009b
α (deg)	16.80	14.43	14.08	15.56	14.38	11.01	9.50
Be (nT)	4.95	4.72	4.36	4.27	4.11	3.97	3.91
rTS (au)	88.0	87.0	86.0	85.0	84.0	82.0	80.0
λ∥ (au)	0.438	0.465	0.506	0.539	0.557	0.574	0.593
KA0	0.90	0.90	0.90	0.90	0.90	0.90	0.90
PA0 (GV)	0.90	0.90	0.90	0.90	0.90	0.90	0.90
${K}_{\perp r}^{0}$	0.02	0.02	0.02	0.02	0.02	0.02	0.02
${K}_{\perp \theta }^{0}$	0.02	0.02	0.02	0.02	0.02	0.02	0.02
c1	0.00	0.00	0.00	0.00	0.00	0.00	0.00
c2∥	2.25	2.25	2.25	2.25	2.25	2.25	2.25
c2⊥	1.688	1.688	1.688	1.688	1.688	1.688	1.688
c3	2.50	2.50	2.50	2.50	2.50	2.70	2.70
Pk (GV)	0.585	0.575	0.565	0.58	0.58	0.57	0.57
d⊥θ	6.00	6.00	6.00	6.00	6.00	6.00	6.00

Download table as:  ASCIITypeset image

Richardson & Wang (2011) reported that because of the dynamic nature of the heliosphere, the position of the termination shock (TS) also varies with solar activity. So, in addition to the solar activity–related variables mentioned above, the changing position of the TS is also incorporated in the model, since it affects the modulated intensities at the Earth (see also, e.g., Manuel et al. 2014). In the model, the position of the TS is changed from 88 au in 2006 to 86 au in 2007, 84 au in 2008, and 80 au in 2009 following Potgieter et al. (2015). The calculated moving averages over the propagation time of α and B_e for each semester, along with the position of TS, are tabulated in Table 1. Note that the value of α was less than 10° during the 2009b period (see Figure 5).

The general expression for the diffusion coefficient parallel to the average background HMF is given by

$$\begin{eqnarray}&&{K}_{\parallel }={({K}_{\parallel })}_{0}\beta \left(\displaystyle \frac{{B}_{0}}{B}\right){\left(\displaystyle \frac{P}{{P}_{0}}\right)}^{{c}_{1}}{\left[\displaystyle \frac{{\left(\tfrac{P}{{P}_{0}}\right)}^{{c}_{3}}+{\left(\tfrac{{P}_{k}}{{P}_{0}}\right)}^{{c}_{3}}}{1+{\left(\tfrac{{P}_{k}}{{P}_{0}}\right)}^{{c}_{3}}}\right]}^{\tfrac{{c}_{2\parallel }-{c}_{1}}{{c}_{3}}},\end{eqnarray} \tag{ 14 }$$

where ${({K}_{\parallel })}_{0}$ is a scaling constant in units of 10²² cm² s⁻¹, with the rest of the equation written to be dimensionless with P₀ = 1 GV and B₀ = 1 nT (in order to keep the same cm² s⁻¹ unit). Here c₁ is a power index that may change with time if required; c_2∥ and c_2⊥ (c_2⊥ = 0.75 c_2∥), together with c₁, determine the slope of the rigidity dependence, respectively, above and below a rigidity with the value P_k, which may change with time if required (see Table 1); and c₃ determines the smoothness of the transition. The rigidity dependence of K_∥ is thus a combination of two power laws; P_K determines the rigidity where the break in the power law occurs, and the value of c₁ determines the slope of the power law at rigidities below P_k. The radial dependence of the diffusion coefficients in the inner heliosphere (less than 5 au) was adjusted according to Vos & Potgieter (2015, 2016) in order to reproduce the radial gradients as observed by Ulysses (Gieseler & Heber 2016), which requires an increase in the diffusion coefficients at these radial distances.

Perpendicular diffusion in the radial direction is assumed to scale spatially similarly to Equation (14) but with a different rigidity dependence at higher rigidities; in the case of positrons and electrons, the value of c₁ = 0.0, so the term ${(P/{P}_{0})}^{{c}_{1}}$ will be unity, so that

$$\begin{eqnarray}&&{K}_{\perp r}=0.02{({K}_{\parallel })}_{0}\beta \left(\displaystyle \frac{{B}_{0}}{B}\right){\left[\displaystyle \frac{{\left(\tfrac{P}{{P}_{0}}\right)}^{{c}_{3}}+{\left(\tfrac{{P}_{k}}{{P}_{0}}\right)}^{{c}_{3}}}{1+{\left(\tfrac{{P}_{k}}{{P}_{0}}\right)}^{{c}_{3}}}\right]}^{\tfrac{{c}_{2\perp }}{{c}_{3}}}.\end{eqnarray} \tag{ 15 }$$

This means that K_⊥r = 0.02 K_∥ is a straightforward, reasonable, and widely used assumption.

On the other hand, the polar perpendicular diffusion (K_⊥θ) is anisotropic in the inner heliosphere, with the consensus that K_⊥θ > K_⊥r away from the equatorial regions, as discussed by Potgieter (2000) and Potgieter et al. (2014).

The perpendicular diffusion coefficient in the polar direction is given by

$$\begin{eqnarray}&&{K}_{\perp \theta }=0.02{K}_{\parallel }{f}_{\perp \theta }={K}_{\perp r}{f}_{\perp \theta },\end{eqnarray} \tag{ 16 }$$

with

$$\begin{eqnarray}&&{f}_{\perp \theta }={A}^{+}\mp {A}^{-}\tanh [8({\theta }_{A}-90^\circ )\pm {\theta }_{F}].\end{eqnarray} \tag{ 17 }$$

Here ${A}^{\pm }=({d}_{\perp \theta }\pm 1)/2$, θ_F = 35°, θ_A = θ for θ ≤ 90° but θ_A = 180°–θ for θ ≥ 90°, and d_⊥θ = 6.0 (in this study; see Table 1). This means that the polar perpendicular diffusion K_⊥θ can be enhanced toward the heliospheric poles by a factor d_⊥θ; we assumed d_⊥θ = 6 with respect to the value K_∥ in the equatorial region of the heliosphere. For the motivation of this approach, see Potgieter (2000), Ferreira et al. (2003), Ngobeni & Potgieter (2011), and Potgieter et al. (2015).

The modulation parameter values required to reproduce the PAMELA positron spectra, discussed next, for the period 2006b to 2009b are summarized in Table 1. The perpendicular diffusion in the radial and polar directions is represented by ${K}_{\perp r}^{0}$ and ${K}_{\perp \theta }^{0}$ in Table 1. The rigidity dependence of the diffusion coefficients and drift coefficient at the Earth are illustrated together in Figure 6. The diffusion coefficients (K) are related to the corresponding mean free paths (MFPs) (λ) by K = λ(v/3), where v is the speed of the positrons. Figure 6 shows the rigidity dependence of λ_∥ and ${\lambda }_{\perp }\equiv {\lambda }_{\perp r}\equiv {\lambda }_{\perp \theta }$, as well as the drift scale (λ_A) at Earth. We have attempted to keep most of these variables the same for all seven semesters (2006b–2009b), except λ_∥, c₃, and P_k, as shown in Table 1.

Zoom In Zoom Out Reset image size

Concerning the diffusion of positrons in the heliosphere, directly related to what is depicted in Figure 6, Teufel & Schlickeiser (2003) showed that for positrons, the parallel MFP at Earth differs from protons in terms of their rigidity; it remains essentially independent of rigidity below almost 100 MV, whereas for protons, it keeps decreasing strongly with decreasing rigidity. We assume that these differences persist throughout the heliosphere. There is, however, no difference between the MFPs for electrons and positrons. For reasons of a fundamental nature on why these MFPs differ for positrons and protons, based on various turbulence theory approaches, see the original work by Bieber et al. (1994), with an update by Engelbrecht & Burger (2013) and a comprehensive review by Shalchi (2009).

4.2. Comparison of Modeling Results with Observations

The Parker TPE is solved for each semester (2006b to 2009b) by using the calculated average values of the HCS tilt and HMF magnitude at the Earth and by adjusting the diffusion and drift coefficients. We successfully reproduced the 6 month–averaged PAMELA positron spectra from 2006 July to 2009 December by carefully adjusting the diffusion and drift coefficients.

Figure 7 shows the observed positron spectra from Figure 1 (colored circles) overlaid by the corresponding computed spectra (colored solid lines) at the Earth with respect to the LIS at 122 au. The modulation of positrons becomes significant below 30 GeV and increases gradually with decreasing energies, similar to other CR particles. The effect of continuously varying heliospheric modulation conditions on galactic positrons is also evident from below 10 GeV. The positron spectrum for 2006 July–December (2006b; blue line) is affected the most, decreasing gradually up to 2006 July–December (2009b; red line), in accord with the gradual decrease in solar activity. These modulated positrons exhibit a characteristic peak in each spectrum just below 1 GeV in all seven semesters. The kinetic energy where the peaks occur gradually shifts to lower values with decreasing modulation as more low-energy positrons reach the Earth. The spectra decrease with lowering energies to become a minimum around 100–200 MeV; for 2006b, this energy value is meaningfully lower than for 2009b. Below 100 MeV, the spectra turn up as the modulation becomes diffusion-dominated. The match between the computed and observed spectra is most reasonable down to about 200 MeV but clearly less so at lower energies.

Zoom In Zoom Out Reset image size

The computed intensity ratios of each semester with respect to 2006b (solid lines), along with the observed positron intensity ratio with respect to 2006b (data points with error bars) adopted from Figures 2 and 7, are shown in Figure 8. Consistent with the results in the previous figure, the intensity ratio with respect to 2006b stays ≈1 down to about 20 GeV, then deviates quickly and significantly below ≈5 GeV. The confidence level of the observed positron spectra below 200 MeV is relatively low, as the total error bars indicate. The intensity ratio increases gradually up a maximum value around 30–40 MeV. Below this energy, it shows a continuous relatively small decrease. As expected, the ratio with respect to 2006b is highest for 2009b (end of the solar activity minimum), and it decreases for each semester (2009a to 2007a) to unity for 2006b. For example, at 30 MeV, the intensity of 2009b was increased by a factor of ≈6.5 with respect to 2006b; at 100 MeV, it was ≈4.5; and at 1 GeV, it was ≈2.0.

Zoom In Zoom Out Reset image size

The ratio between the positron LIS and the modulated spectra at Earth tells how much modulation happens inside the heliosphere; this ratio is called the modulation factor (MF), as shown in Figure 9. The MF is calculated as a function of kinetic energy for each semester by taking the ratio of the computed intensity at the Earth with respect to the LIS at 122 au. Evidently, the spectrum for 2006b has been modulated the most, with the largest MF about 0.0003 around 60 MeV. At lower energies, the MF becomes steady as the modulated spectra assume the spectral shape of the LIS.

Zoom In Zoom Out Reset image size

As a further illustration of positron modulation, the MF calculated for the computed modulated spectrum of 2006b is shown in Figure 10 at different radial locations, from the Earth (1 au) to 100 au from the Sun, with the heliopause at 122 au.

Zoom In Zoom Out Reset image size

4.3. The Effect of Drifts and HMF Polarity on Positrons

The solar minimum between solar cycles 23 and 24 is a so-called negative (A < 0) polarity period. The Sun reversed this polarity during the solar maximum period so that it completely turned to an A > 0 polarity in 2014 April (Sun's polar magnetic field data source: http://wso.stanford.edu). According to drift theory, when this happens, the drift patterns of both positively and negatively charged CR particles change direction; during an A > 0 cycle, the positrons drift inward toward the Earth mainly through the polar regions of the heliosphere and then outward mostly along the HCS.

Figure 11 illustrates what happens to the modulated positron spectra at the Earth when this HMF direction reversal is implemented in the model. It shows the computed positron spectra for 2006b (blue solid line) and 2009b (red solid line), both being in an A < 0 epoch, compared to the corresponding spectra for an A > 0 epoch (blue and red dashed lines), assuming similar modulation conditions, except that only the HMF direction was changed and, consequently, the drift directions. Evidently, a higher spectral intensity for positrons is predicted during this A > 0 period compared to the A < 0 period, down to about 20 MeV, where the modulation process become diffusion-dominated for positrons (see Figure 6) and the drift effects subside. Note that for the modulated spectra above 10 GeV, drift effects are still there but become increasingly insignificant to gradually subside. In the case of positrons, drift effects are a maximum between about 100 MeV and 2 GeV, depending on the level of modulation.

Zoom In Zoom Out Reset image size

We are also interested in looking at how the positron spectra vary radially from the Earth to the heliopause during both A < 0 and A > 0 polarity. This is illustrated in Figure 12, showing the computed positron intensity at different radial distances from 1 au to the heliopause at 122 au during the 2006b period: at 1 au (Earth), 10, 50, and 100 au (colored solid lines). The colored dashed lines are the corresponding computed positron intensity for an A > 0 polarity epoch with similar modulation conditions as in 2006b.

Zoom In Zoom Out Reset image size

The drift effect on positrons decreases gradually from Earth to the heliopause. The difference between the A > 0 and A < 0 spectra is a maximum at Earth, decreasing with increasing radial distance, becoming much less at 100 au. In terms of energy, the maximum difference shifts to higher energies with increasing radial distances. Note that the computed A > 0 and A < 0 spectra cross between 1 and 2 GeV at the Earth. This crossover energy shifts to lower energies with increasing radial distance. In principle, this crossing of spectra is required to explain why the intensity at NM energies during A > 0 epochs is lower than during A < 0 epochs, in contrast to what occurs for protons at lower energies, where A > 0 intensities are always higher than during A < 0 cycles. Qualitatively, this means that more low-energy positrons reach Earth during A > 0 cycles, so that the corresponding spectrum is softer than during A < 0 cycles, causing the peak in the spectral intensity to shift to lower energies. This drift effect was also shown for protons and electrons by Potgieter & Vos (2017; see their Figure 6). For a detailed discussion of this intriguing effect for Galactic protons, see Potgieter & Moraal (1985) and Reinecke & Potgieter (1994). As seen before, the drift effects vanish at lower energies, whereas at high energies, the effects gradually subside.

The primary objective of this study is the reproduction of 6 month–averaged positron spectra as observed by PAMELA from 2006 July to 2009 December. This is accomplished by using a well-established and comprehensive 3D numerical modulation model that had been used previously to study the heliospheric modulation of monthly averaged Galactic proton spectra observed by PAMELA for the same period over the energy range 80 MeV to 30 GeV (Potgieter et al. 2014; Vos & Potgieter 2015; Potgieter & Vos 2017) and for the modulation of 6 month–averaged Galactic electron spectra observed by PAMELA for the same period over the same energy range, which is important for solar modulation (Potgieter et al. 2015; Potgieter & Vos 2017).

When simulating the modulation of CRs, the tilt angles (α) of the HCS and HMF strength at Earth (B_e) are considered very good proxies for solar activity. The moving average values of α and B_e are used as the essential time-varying modulation parameters, used as input to the numerical model along with the position of the TS to reproduce the positron PAMELA spectra; these parameters are shown in the first three rows of Table 1. Changing B_e directly affects the value of the diffusion and drift coefficients, as in Equations (9), (14), and (15). However, in order to reproduce the observed PAMELA positron spectra, the value of the diffusion coefficients had to be increased for each semester from 2006b to 2009b, that is, in addition to the changes in the three mentioned parameters. A similar result was reported by Potgieter et al. (2014) for the proton spectra.

The ratio of the perpendicular diffusion coefficients to the parallel diffusion coefficient is kept at 0.02, which gives K_⊥r = K_⊥θ = 0.02 K_∥, but an additional enhancement of a factor of d_⊥θ = 6 for K_⊥θ in the polar direction is also required. We increased the value of λ_∥, and consequently also λ_⊥r and λ_⊥θ, by a factor of ≈1.35 from 2006b to 2009b, as shown in the fourth row of Table 1.

The value of K_A0 = 0.90 and P_A0 = 0.90 GV for all seven periods means that we reduced drift in an implicit way through increasing K_∥, K_⊥r, and K_⊥θ without changing K_D (in Equation (9)). Furthermore, c₁, together with c₂ (c_2∥ and c_2⊥), determines the slope of rigidity dependence above and below P_k.

In the case of positrons, c₁ = 0, c_2∥ = 2.25, and c_2⊥ =1.688 (0.75c_2∥) are kept constant for all seven periods; only the value of P_k is changed when required, to reproduce the seven positron spectra. The smoothness of the transition in rigidity dependence is determined by c₃ = 2.50 for 2006b to 2008b and c₃ = 2.70 for 2009a and 2009b, as shown in Figure 6 and listed in Table 1.

Although particle drift played a significant role (Di Felice et al. 2017; Munini et al. 2018) in the modulation of electrons and positrons during the unusual 2006–2009 minimum period, the process was diffusion-dominated. As such, there was an unexpected result for an A < 0 solar minimum epoch (phase). Shown in Figure 11 is that drift effects below 20 MeV and above 10 GeV are found to be negligible. Also shown is that the drift effect on spectra at Earth is found to be the lowest during the 2009 July–December period (when solar activity was also at its lowest).

Figures 10 and 12 illustrate the total modulation of positrons from the heliopause to the Earth. Our 3D model predicts that the intensity of positrons over the energy range 1 MeV to at least 1 GeV will be higher during the upcoming A > 0 solar minimum period than during 2006–2009, assuming similar quiet modulation conditions.

The main modulation features of Galactic positrons in the heliosphere are shown in Figures 11 and 12. First, we notice that the very LIS flattens off below about 30 MeV. The modulated spectra at the Earth depict the same trend but at a lower energy, because these electrons also experience a modest amount of adiabatic energy loss (Moraal & Potgieter 1982). The total amount of modulation as a function of energy between the LIS, specified at the heliopause (122 au), and the modulated spectra (1 au) is thus almost the same at these low energies, also shown quantitatively in Figure 9. This is true for positron modulation throughout the heliosphere, as illustrated in Figure 10. The main reason for this is because the MFPs for positrons are independent of rigidity at these low energies, so that the diffusion process effectively dominates the total modulation. Above this energy, this is no longer the case; particle drift also begins to play a role, so that the shape of the modulated spectra becomes significantly different, causing relatively more modulation to occur, and is significantly influenced by drift, the contribution of which varies from the A > 0 to A < 0 cycles. Above a few GeV, particle drift, as well as the total amount of modulation, declines rapidly with increasing kinetic energy.

The next phase of this study is to apply this model to Galactic electron and positron spectra observed by AMS-02 (Aguilar et al. 2018) for the period 2011–2017, as was done for PAMELA protons (Martucci et al. 2018). These AMS-02 observations, together with the PAMELA observations for 2006–2009, can provide electron and positron spectra over a full solar cycle. It will help to understand how the modulation of positrons is different from that for electrons over a complete solar cycle. This may also clarify the role of drift over a complete solar cycle, especially how large the drift effects are during an HMF reversal period (see also Heber et al. 2002; Adriani et al. 2016).

The 6 month–averaged CR positron spectra of energy range 80 MeV to 30 GeV, recorded by PAMELA from 2006 July to 2009 December (Munini 2015), were used together with a comprehensive 3D numerical modulation model, which includes particle drift and a solar cycle–dependent dynamic heliosheath, to understand how the solar modulation of galactic positrons is affected by the unusual solar minimum conditions of solar cycles 23/24. For this purpose, a new very LIS for Galactic positrons was constructed.

It follows from the observations that the positron intensity has increased gradually from 2006 July–December to reach its highest level during 2009 July–December, in accord with decreasing solar activity. This recovery toward solar minimum modulation is effectively reproduced by the model using the same approach as for the modulation of protons and electrons observed by PAMELA. An effort was made to keep the modeling approach as simple as possible within the context of a 3D drift model but to adhere to the main features of CR modulation as also observed by previous space missions, such as Ulysses. Essentially, all three diffusion coefficients scale as 1/B, with a small modification over the first 5 au for the radial dependence of K_∥ but a much larger modification to K_⊥θ in terms of its polar dependence. Apart from changing the tilt angle of the HCS, the HMF magnitude, and the position of the TS in a time-varying manner, we found that it was necessary to additionally increase, e.g., λ_∥ at the Earth from 0.438 au to 0.593 au at 1 GV in order to reproduce the positron spectra from 2006b to 2009b.

Comparison of simulated spectra with precise observations makes it possible to extend the modeling with confidence down to 1 MeV for positrons where observations are unavailable. It follows that the LIS becomes flat, essentially independent of energy below 50 MeV. Because positron modulation is dominated by diffusion at these low energies (see Figure 6), similar to electron modulation, the spectral shape of the LIS at these energies is maintained to deep inside the heliosphere (see Figures 4, 11, and 12). This implies that if the positron spectrum could be observed at these energies at the Earth, it would tell what the spectral shape of the very LIS is. This is in sharp contrast to protons, antiprotons, and all CR nuclei, the spectra of which are determined by adiabatic energy losses at low energies at the Earth.

It is interesting to note that the highest positron spectrum was observed by PAMELA during the second half of 2009 (2009 July–December), similar to Galactic electrons and protons (highest intensity in 2009 December), for the latter the highest since the beginning of the space era. This was unexpected for positively charged CR particles (protons, positrons, etc.) during an A < 0 polarity epoch because drift theory generally predicts lower spectra for positively charged CRs during A < 0 cycles. As such, the unusual 2006–2009 minimum solar activity period seems to have broken all predictions and expectations caused by the Sun, which unexpectedly became relatively very quiet and seems to continue to be so. If similar modulation conditions prevail during the upcoming solar minimum period, an A > 0 cycle with different drift patterns than during 2006–2009, the spectra for protons, positrons, and all heavy nuclei should be even higher than in 2009 December. If the Sun is even quieter than in 2006–2009, new record CR intensities will then be set.

The authors from SA express their gratitude for the partial funding by the South African National Research Foundation (NRF) under grant 98947. The authors wish to thank the GALPROP developers and their funding bodies for access to and use of the GALPROPWebRun service. We also acknowledge the use of HCS tilt data from the Wilcox Solar Observatory's website, http://wso.stanford.edu, and HMF data from NASA's OMNIWEB data interface, http://omniweb.gsfc.nasa.gov. O.P.M.A. and D.B. acknowledge the partial financial support from the postdoctoral program of the North-West University, South Africa. O.P.M.A. also thanks Jan-Louis Raath for many useful discussions.