New Very Local Interstellar Spectra for Electrons, Positrons, Protons, and Light Cosmic Ray Nuclei

Because of solar modulation, the local interstellar spectra (LISs) for cosmic rays (CRs) are still not fully determined over all energy ranges. A breakthrough came when Voyager 1 (V1) made in situ CR measurements beyond the heliopause in 2012 August for the first time (Gurnett et al. 2013; Stone et al. 2013; Webber & McDonald 2013; Cummings et al. 2016). With its unique positioning, the V1 observations have become central to recent lower energy LIS studies, allowing for galactic electrons, protons, and light nuclei (specifically helium, carbon, boron, and oxygen) LISs to be determined with vastly more confidence down to a few MeV/nuc (about 3 MeV/nuc for protons and nuclei, and about 2.7 MeV for electrons Cummings et al. 2016).

For much higher energies very precise CR spectra are observed at the Earth by experiments such as the satellite-borne PAMELA instrument (Adriani et al. 2011a, 2014b; Boezio et al. 2015, 2017). PAMELA made highly valuable observations during the 2006–2009 solar minimum, with specifically the end of 2009 showing unusually high CR intensities (see, e.g., Heber et al. 2009; Mewaldt et al. 2010; Aslam & Badruddin 2012; Lave et al. 2013; Potgieter et al. 2014). These PAMELA observations can be considered ideal for estimating the CR LISs, as they cover a large range of energies, have high statistics and include not only measurements of CR electrons and protons, but also antimatter and light nuclei. The observations made by PAMELA are over an energy range of about 80 MeV/nuc–50 GeV/nuc for protons and 80 MeV to 40 GeV for electrons. For energies above about 30 GeV/nuc, the CR spectra can be considered directly observed LISs because modulation is largely negligible here, while below this energy the effects of modulation become increasingly pronounced and cannot be neglected. The PAMELA observations of interest to this study are those made at the end of 2009 for electrons (Adriani et al. 2015), positrons (Munini 2016), and protons (Adriani et al. 2013b), as well as those presented as averaged over the 2006–2009 period for positrons (Adriani et al. 2013a), and over the 2006–2008 period for helium (Adriani et al. 2011b), carbon, and boron (Adriani et al. 2014a).

A comprehensive galactic propagation model, such as the GALPROP code (Strong & Moskalenko 1998; Strong et al. 2007; Porter et al. 2017), can be used to compute galactic spectra for a wide set of CR species. GALPROP has been used for various aspects of CR spectrum studies, including the study of CR electrons (Strong et al. 2011b), CR antiparticles (Moskalenko et al. 2002), galactic turbulence, and CR diffusion (Ptuskin et al. 2006), and the galactic broadband luminosity spectrum (Strong et al. 2010). Other propagation codes have also been used to similarly compute galactic spectra, assumed to be LISs. A simple leaky box model by Webber & Higbie (2009) has been used to compute electron (Webber 2015a; Webber & Villa 2017) and light nuclei (Webber 2015b) LISs separately. Extensive studies on CR isotopes were done by Lave et al. (2013) using this type of model. Other more advanced models focused on specific features such as the effects of discrete CR source distributions by Büsching et al. (2005) and the effect of the galactic spiral arm structure investigated by Büsching & Potgieter (2008) and Kopp et al. (2014). Improving the implementation of physical processes has also been investigated, such as incorporating fully anisotropic diffusion by Effenberger et al. (2012). While codes, such as those presented by Kissmann (2014), are being developed to improve on the numerical capabilities of galactic propagation models, in our opinion GALPROP remains the most sophisticated and comprehensive galactic propagation model. The history of GALPROP, together with the accessibility offered by GALPROP's dedicated WebRun service (Vladimirov et al. 2011), lead the authors to believe that it is the best suited model for an LIS study with a wide scope of CR particles such as presented here.

After implementing the propagation model, the numerically computed LISs can then be matched against available observations in order to verify the spectra. For the V1 and higher energy PAMELA observations the comparison can be done directly, but for Earth-based observations below about 30 GeV/nuc the effects of solar modulation need to be taken into account. In order to do this and compute the CR spectra inside the heliosphere at the Earth's position, the sophisticated 3D modulation code presented by Potgieter & Vos (2017) is used. This model has previously been used to compute electron (Potgieter et al. 2015) and proton (Vos & Potgieter 2015) spectra with older LISs. Modulation models require a spectrum to be specified at the chosen modulation boundary, the heliopause, as the initial input condition (see the reviews by Potgieter 2013, 2017). The propagation models compute a galactic spectrum for a given CR species, which is strictly speaking not an exact representation of the CR intensities in the local vicinity of the heliosphere. When not including local CR sources and local galactic structures, a computed galactic spectra may not be the same as a very local interstellar spectrum (the required heliopause spectrum). However, for this study the computed galactic spectra will be taken as a heliopause spectra for input into the modulation model, and the GALPROP parameters presented here take this into account. For further discussion on galactic spectra, LISs and heliopause spectra, from a solar modulation point of view, see Potgieter (2014b).

This study is considered to be the first to concurrently implement these comprehensive models for CR propagation in the Galaxy and heliosphere simultaneously for electrons, positrons, protons, helium, carbon, boron, and oxygen. Using the relevant observations made by V1 and PAMELA, LISs can be computed with more confidence than before in the energy range of 3 MeV/nuc–100 GeV/nuc, even for the energies that are not directly observed.

The GALPROP code computes the propagation of relativistic charged particles through the Galaxy, by describing this propagation as a diffusive process. The GALPROP model implements the following equation for the propagation of a particular particle species:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial \psi }{\partial t} & = & S({\boldsymbol{r}},p)+{\rm{\nabla }}\cdot (K{\rm{\nabla }}\psi -{\boldsymbol{V}}\psi )\\ & & +\,\displaystyle \frac{\partial }{\partial p}\left[{p}^{2}{K}_{p}\displaystyle \frac{\partial }{\partial p}\displaystyle \frac{1}{{p}^{2}}\psi +\displaystyle \frac{p}{3}({\rm{\nabla }}\cdot {\boldsymbol{V}})\psi -\dot{p}\psi \right]-\displaystyle \frac{\psi }{\tau },\end{array}\end{eqnarray} \tag{ 1 }$$

where ψ = ψ(r, p, t) is the CR density per unit of total particle momentum p at position r (Strong et al. 2007). The source term S(r, p) includes contributions from primary sources, as well as spallation and decay. Spatial diffusion is represented by the coefficient K. The spatially dependent convection velocity is represented by V and determined by the gradient in the galactic wind dV/dz. The ∇ · V term represents the adiabatic momentum gain or loss in the nonuniform flow of gas with a frozen-in magnetic field whose inhomogeneities scatter the CRs. Reacceleration is described as diffusion in momentum space and determined by the coefficient K_p, while $\dot{p}$ is the momentum loss rate and comprises all forms of energy losses, mostly by synchrotron radiation for electrons or ionization loss for protons and heavier nuclei. The catastrophic particle losses are represented by the timescale τ, this includes the timescale for fragmentation (τ_f), which depends on the total spallation cross-section, and the timescale for radioactive decay (τ_r).

In GALPROP the spatial diffusion coefficient K is assumed to be independent of radius r and height z, and is taken as being proportional to a power law in rigidity P:

$$\begin{eqnarray}&&K=\beta {K}_{0}{\left(P/{P}_{0}\right)}^{\delta },\end{eqnarray} \tag{ 2 }$$

where δ = δ₁ for rigidity P < P₀ (the reference rigidity), while δ = δ₂ for P > P₀. Here β = v/c is the speed of particles v at a given rigidity relative to the speed of light c. The magnitude of the diffusion coefficient K₀ is in effect a scaling factor for diffusion, generally with units of 10²⁸ cm² s⁻¹. When considering reacceleration, the momentum-space diffusion coefficient K_p is estimated as related to K so that ${K}_{p}K\propto {p}^{2}{V}_{A}^{2}$, with V_A the Alfveń wave speed with units of km s⁻¹.

The CR sources are assumed to be concentrated near the galactic disk and have a radial distribution similar to that of supernova remnants, with the distribution assumed to be the same for all CR primaries. The primary contribution to the sources requires an injection spectrum and relative isotopic compositions to be specified. The injection spectrum for nuclei, as input to the source term, is assumed to be a power law in rigidity so that:

$$\begin{eqnarray}&&S(P)\propto {\left(P/{P}_{\alpha 0}\right)}^{\alpha },\end{eqnarray} \tag{ 3 }$$

for the injected particle density and usually contains a break in the power law with index α = α₁ below the source reference rigidity P_α0 and α = α₂ above. Values for α₁ and α₂ are negative and nonzero, thus giving a rigidity dependent injection spectrum. For isotopes considered wholly secondary, such as for positrons, no input spectrum is given and they are set to zero at the sources. For the source abundance values as used in GALPROP, this study keeps the values unchanged from what is used by Ptuskin et al. (2006). The source spectrum is the same for all CR nuclei, but differs for electrons. Similar to that of the CR nuclei, the injection spectrum for electrons is given by:

$$\begin{eqnarray}&&S(P)\propto {\left(P/{P}_{\alpha e0}\right)}^{{\alpha }_{e}},\end{eqnarray} \tag{ 4 }$$

where α_e = α_e1 below the source reference rigidity P_αe0 and α_e = α_e2 above (Strong et al. 2007).

The Galaxy is described as a cylindrical volume for CR propagation studies. This includes a galactic halo, in which CRs have a finite chance to return to the galactic disk. Assuming symmetry in azimuth leads to a two-dimensional (2D) model of space that depends only on galactocentric radius and height, additionally, neglecting time dependence leads to a steady-state model. When implemented in the GALPROP code this gives a 2D model with radius r, the halo height z above the galactic plane, and symmetry in the angular dimension in galactocentric-cylindrical coordinates. The propagation region is bounded by r = R and z = ±H, beyond which free escape is assumed. For this study the halo height is fixed at 4 kpc, because varying its size can simply be counteracted by directly varying the diffusion coefficient. While the GALPROP code has been designed for the propagation of CRs on either a 2D or 3D spatial grid (Strong & Moskalenko 2001), only the 2D model is considered for this study, that is, two spatial dimensions and momentum, giving the basic coordinates (r, z, p) for the rotationally symmetric cylindrical grid. Symmetry is assumed above and below the galactic plane in order to save on the computational requirements of the code. The GALPROP code solves the propagation equation, for each of the CR species that are taken into account, using a Crank–Nicholson implicit second-order scheme. The processes are described by differential operators in the propagation equation and these operators are implemented as finite differences for each dimension (r, z, p) in the numerical scheme. For extensive details on solving the propagation equation, the numerical scheme and differential operators see Strong & Moskalenko (1998), Strong et al. (2007), and Strong et al. (2011a). The computational runs for this study are done via the GALPROP WebRun service, which offers the benefits of running the most recent version of GALPROP, with error detection, powerful computing power, and user support. The service can be accessed at http://galprop.stanford.edu/webrun. The details on the implementation of the WebRun service, the updated features of the code and the computer cluster specifications are presented by Vladimirov et al. (2011).

In order to compute the spectra of CRs as they arrive at the Earth, after their transport through the heliosphere, a numerical modulation model is used. This full 3D model, as described by Potgieter & Vos (2017), can compute modulated differential intensities throughout the heliosphere by implementing the CR transport equation first derived by Parker (1965). This model takes into account the major modulation mechanisms of convection and adiabatic energy losses due to the expanding solar wind, particle diffusion, and drifts due to the heliospheric magnetic field (HMF). This steady-state model also includes a wavy current sheet and a heliosheath, but does not consider shock acceleration at the termination shock. For the purposes of this work, we restrict the modulation, with a few exceptions, to the PAMELA observations made during the solar minimum of 2009, which was an A < 0 solar magnetic field epoch.

The CR transport equation by Parker (1965) can be written in terms of rigidity (P) as:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial f}{\partial t} & = & -{{\boldsymbol{V}}}_{\mathrm{sw}}\cdot {\rm{\nabla }}f-\langle {{\boldsymbol{v}}}_{D}\rangle \cdot {\rm{\nabla }}f+{\rm{\nabla }}\cdot ({{\boldsymbol{K}}}_{s}\cdot {\rm{\nabla }}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{ 5 }$$

where $f({\boldsymbol{r}},P,t)$ is the CR distribution function at time t and at vector position r. For the left side of the equation $\tfrac{\partial f}{\partial t}=0$ as we consider only a steady-state solution for the solar minimum conditions where the modulation parameters only gradually change over time. The first term on the right gives the outward convection due to the solar wind. The second term represents the averaged particle drift velocity $\langle {{\boldsymbol{v}}}_{D}\rangle$

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

where K_A is the generalized drift coefficient, and B is the HMF vector with magnitude B. The third term describes the spatial diffusion caused by the scattering of CRs, where K_s is the symmetric diffusion tensor, and the last term represents the adiabatic energy change, which depends on the sign of the divergence of V_sw. If (∇ · V_sw) > 0 adiabatic energy losses occur, as is the case in most of the heliosphere, except inside the heliosheath where we assume that (∇ · V_sw) = 0. For a modified Parker-type HMF, B_m with magnitude B_m, such as the Smith–Bieber modification (Smith & Bieber 1991), the drift coefficient K_A can be written as:

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

The dimensionless constant K_A0 ranges from 0.0 to 1.0, where K_A0 = 1.0 is called 100% drift or full weak scattering. In this study K_A0 is kept at 0.90, effectively setting particle drift to a 90% level, as this has been found in previous studies to best fit the observations (see, e.g., Aslam et al. 2019). For detailed discussions of this process, see also, for example, Ngobeni & Potgieter (2015), Nndanganeni & Potgieter (2016), and Raath et al. (2016).

The symmetric diffusion tensor K_s is comprised of three diffusion coefficients, ${K}_{\parallel }$, ${K}_{\perp r}$, and K_⊥θ. The expression for the diffusion coefficient parallel to the average background HMF is given by:

$$\begin{eqnarray}&&{K}_{\parallel }={\left({K}_{\parallel }\right)}_{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{ 8 }$$

with ${({K}_{\parallel })}_{0}$ a scaling constant in units of 10²² cm²s⁻¹, P₀ = 1 GV, and B₀ = 1 nT. The power indices c₁ and ${c}_{2\parallel }$ respectively determine the slope of the rigidity dependence of ${K}_{\parallel }$ above and below the rigidity P_k, while c₃ determines the smoothness of the transition. Perpendicular diffusion in the radial direction (K_⊥r) is assumed to scale spatially similar to Equation (8), but with a different rigidity dependence at higher rigidities so that:

$$\begin{eqnarray}&&{K}_{\perp r}=0.02{\left({K}_{\parallel }\right)}_{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\perp }-{c}_{1}}{{c}_{3}}},\end{eqnarray} \tag{ 9 }$$

which is close to the widely used assumption of ${K}_{\perp r}=0.02\,{K}_{\parallel }$. The polar perpendicular diffusion coefficient (K_⊥θ) is given by:

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

The latitudinal dependence (f_⊥θ) in the above equation is given by:

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

where ${A}^{\pm }=({d}_{\perp \theta }\pm 1)/2$, θ_F = 35°, θ_A = θ for θ ≤ 90°, and θ_A = 180°–θ for θ ≥ 90°. With this expression K_⊥θ can be enhanced toward the heliospheric poles by a factor of d_⊥θ. For examples of modulated spectra computed using this model, see the reviews by Potgieter (2013, 2017), for a comprehensive discussion of charge-sign dependent modulation in the heliosphere, see Potgieter (2014a), and for details on the 3D model described in short above, see also Potgieter et al. (2014, 2015), Vos & Potgieter (2015), and Aslam et al. (2019).

With this model the proton LIS can be modulated to compute the corresponding spectrum at the Earth. The relevant modulation parameter values required to reproduce the observed PAMELA spectra for the last half of 2009, indicated as 2009b, are summarized in Table 1 (for 1.0 GV at the Earth) and are adjusted from those presented by Vos & Potgieter (2015). The diffusion coefficients are related to the corresponding mean free paths (MFPs) by, K = λ(v/3), where v is the particle speed, giving rigidity dependent MFPs: ${\lambda }_{\parallel }$ (blue curve), λ_⊥r (red curve), λ_⊥θ (green curve), and the drift scale (yellow curve) as shown in Figure 1. Similarly the light nuclei LISs can be modulated using the same parameters and accounting for the specific particle masses and charges. For electrons and positrons the parameters are similarly implemented, as presented by Potgieter et al. (2015), and have been adjusted to the values listed in Table 1 with the MFPs shown for protons in Figure 1 and for electrons in Figure 2. The modulated spectra will be shown together with the corresponding LISs in the rest of the figures in the next section.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Summary of Parameters Used in the Modulation Model (at 1.0 GV at the Earth) for the 2009b Time Period

Parameters	Electrons and	Protons and
	Positrons	Light Nuclei
${\lambda }_{\parallel 0}$ (au)	0.593	1.185
KA0	0.90	0.90
PA0 (GV)	0.90	0.90
c1	0.00	0.70
${c}_{2\parallel }$	2.25	1.52
c2⊥	1.688	1.14
c3	2.70	2.50
Pk (GV)	0.57	4.00
d⊥θ	6.00	6.00

Download table as:  ASCIITypeset image

The LISs for electrons, positrons, and protons were computed first using the values listed by Ptuskin et al. (2006) as a reference. The parameters were then systematically adjusted in order to find LISs that directly match the mentioned observations (V1 and PAMELA 2009b period) for electrons, positrons, and protons. The observed B/C ratio was also used as a constraint to the parameter sets, while verifying that the helium, carbon, boron, and oxygen LISs also give reasonable matches to the corresponding observations. Changes to the electron and nuclei source indices for lower energies were found necessary in order to update the GALPROP models to take into account the V1 observations. With diffusion having the largest effect on the LISs, changes to the diffusion parameters were also required to reproduce the observations well, instead of just a rough reproduction. Initial studies showed that a GALPROP plain diffusion model was sufficient when studying electrons (as presented by Bisschoff & Potgieter 2014), protons and helium LISs separately (as presented by Bisschoff 2018). Modeling protons, helium, carbon, and the B/C ratio simultaneously, as shown by Bisschoff & Potgieter (2016), required the inclusion of reacceleration in galactic space in the plain diffusion model. To further include electrons and positrons, and thus having a model that can simultaneously compute the required set of CR LISs and match the observations, convection also needed to be taken into account.

The LISs computed with GALPROP for electrons, positrons, protons, and the light nuclei are shown in Figures 3 to 10. The parameters used in GALPROP are listed in Table 2, these values were iteratively found and based on the parameter values presented by Ptuskin et al. (2006) and Strong et al. (2010). Additional changes in GALPROP are made to the source abundances of helium and carbon. The PAMELA observations necessitate decreasing the ¹²C₆ abundance by 10% and increasing the ⁴He₂ abundance by 5%.

Zoom In Zoom Out Reset image size

The computed electron LIS, shown in Figure 3 together with the corresponding V1 observations and the modulated spectrum at the Earth in comparison with PAMELA observations, is approximated over the energy range of 4 MeV–100 GeV by the following expression:

$$\begin{eqnarray}\begin{array}{rcl}{J}_{\mathrm{elec}}(E) & = & 255.0\,\displaystyle \frac{1}{{\beta }^{2}}{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{-1}{\left(\displaystyle \frac{E/{E}_{0}+0.63}{1.63}\right)}^{-2.43}\\ & & +\,6.4{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{2}{\left(\displaystyle \frac{E/{E}_{0}+15.0}{16.0}\right)}^{-26.0},\end{array}\end{eqnarray} \tag{ 12 }$$

where the CR intensity J_elec(E) (given in part.m⁻² s⁻¹ sr⁻¹ GeV⁻¹) is a function of kinetic energy E (given in GeV), E₀ = 1 GeV and again with β = v/c. The modulation parameters are listed in Table 1.

The computed positron LIS and the corresponding modulated spectrum at the Earth, in comparison with PAMELA observations, are shown in Figure 4. The LIS is approximated by the following expression:

$$\begin{eqnarray}\begin{array}{rcl}{J}_{\mathrm{posi}}(E) & = & 25.0\,\displaystyle \frac{1}{{\beta }^{2}}{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{0.1}{\left(\displaystyle \frac{{\left(E/{E}_{0}\right)}^{1.1}+{0.2}^{1.1}}{1+{0.2}^{1.1}}\right)}^{-3.31}\\ & & +\,23.0{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{0.5}{\left(\displaystyle \frac{E/{E}_{0}+2.2}{3.2}\right)}^{-9.5}.\end{array}\end{eqnarray} \tag{ 13 }$$

Zoom In Zoom Out Reset image size

The computed proton LIS and the corresponding modulated spectrum at the Earth, shown in Figure 5 in comparison with the corresponding V1 and PAMELA observations, is approximated over the energy range 4 MeV/nuc–100 GeV/nuc by the following expression:

$$\begin{eqnarray}\begin{array}{rcl}{J}_{{\rm{p}}}(E) & = & 2620.0\,\displaystyle \frac{1}{{\beta }^{2}}{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{1.1}{\left(\displaystyle \frac{{\left(E/{E}_{0}\right)}^{0.98}+{0.7}^{0.98}}{1+{0.7}^{0.98}}\right)}^{-4.0}\\ & & +\,30.0{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{2}{\left(\displaystyle \frac{E/{E}_{0}+8.0}{9.0}\right)}^{-12.0},\end{array}\end{eqnarray} \tag{ 14 }$$

where the CR intensity J_p(E) (given in part.m⁻² s⁻¹ sr⁻¹(GeV/nuc)⁻¹) is a function of kinetic energy E (given in GeV/nuc), E₀ = 1 GeV/nuc, and with β = v/c.

Zoom In Zoom Out Reset image size

Similarly, the computed LISs for helium, carbon, boron, and oxygen are approximated by the following expressions. The helium LIS, shown in Figure 6 together with the corresponding modulated spectrum at the Earth and relevant observations, is approximated by:

$$\begin{eqnarray}&&{J}_{\mathrm{He}}(E)=163.4\,\displaystyle \frac{1}{{\beta }^{2}}{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{1.1}{\left(\displaystyle \frac{{\left(E/{E}_{0}\right)}^{0.97}+{0.58}^{0.97}}{1+{0.58}^{0.97}}\right)}^{-4.0}.\end{eqnarray} \tag{ 15 }$$

The carbon LIS, shown in Figure 7 together with the corresponding modulated spectrum at the Earth and relevant observations, is approximated by:

$$\begin{eqnarray}&&{J}_{{\rm{C}}}(E)=3.3\,\displaystyle \frac{1}{{\beta }^{2}}{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{1.22}{\left(\displaystyle \frac{{\left(E/{E}_{0}\right)}^{0.9}+{0.63}^{0.9}}{1+{0.63}^{0.9}}\right)}^{-4.43}.\end{eqnarray} \tag{ 16 }$$

The boron LIS, shown in Figure 8 together with the corresponding modulated spectrum at the Earth and relevant observations, is approximated by:

$$\begin{eqnarray}\begin{array}{rcl}{J}_{{\rm{B}}}(E) & = & \displaystyle \frac{1}{{\beta }^{2}}{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{1.7}{\left(\displaystyle \frac{E/{E}_{0}+0.685}{1.685}\right)}^{-4.8}\\ & & +\,(3.0\times {10}^{-4}){\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{3}{\left(\displaystyle \frac{E/{E}_{0}+0.204}{1.204}\right)}^{-11.0}.\end{array}\end{eqnarray} \tag{ 17 }$$

These two LISs are used to calculate the B/C ratio for the LISs as shown in Figure 9, in comparison with the ratio at the Earth calculated from the modulated spectra. Lastly, the oxygen LIS, shown in Figure 10 together with the corresponding modulated spectrum at the Earth and relevant observations, is approximated by:

$$\begin{eqnarray}&&{J}_{{\rm{O}}}(E)=3.3\,\displaystyle \frac{1}{{\beta }^{2}}{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{1.23}{\left(\displaystyle \frac{{\left(E/{E}_{0}\right)}^{0.86}+{0.62}^{0.86}}{1+{0.62}^{0.86}}\right)}^{-4.62}.\end{eqnarray} \tag{ 18 }$$

A more detailed description of these LISs can be found in Bisschoff (2018).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We present new LISs for CR electrons, positrons, and protons in Figures 3–5 with the corresponding expressions for these LISs in Equations (12)–(14). Similarly the new LISs for the light CR nuclei helium, carbon, boron, and oxygen are shown in Figures 6–10; together with their corresponding expressions Equations (15)–(18). The corresponding B/C ratio is shown in Figure 9. These LISs were all computed with the GALPROP propagation code using a single model and optimized parameter set as described above. The parameter values are tuned to closely match the respective observations made by V1 and PAMELA, simultaneously for all the above listed CRs. In the energy range at which solar modulation has a significant effect (below about 30 GeV/nuc), the LISs cannot be directly compared to observations at the Earth. A 3D solar modulation model that takes into account all relevant heliospherical modulation effects, such as particle diffusion, convection, and drifts, is therefore used to compute the corresponding CR spectra at the Earth (1 au), which are shown with respect to appropriate LISs and relevant observations as indicated in Figures 3–10. In summary, all the new LISs are shown together in Figure 11. From this comparison the relative differences between the computed LISs can be seen. Electrons have the largest intensity below about 100 MeV, while above this energy protons have the largest intensity. This is a consequence of the V1 observations suggesting a power law for the electron LIS, while the observed proton spectrum levels out instead. The positron LIS has an intensity lower than that of carbon and oxygen above about 5 GeV, but with our current results may exceed that of helium below 200 MeV. As expected from the V1 observations, the carbon LIS is seen to be nearly identical to that of the oxygen LIS.

Zoom In Zoom Out Reset image size

In the past the GALPROP plain diffusion model was sufficient when studying electrons (Bisschoff & Potgieter 2014), protons, and helium LISs separately. Including carbon, the B/C ratio, and positrons in the study, the constraints placed on the LISs by observations necessitated also considering reacceleration and convection in galactic space. As with the study by Bisschoff & Potgieter (2016), the B/C ratio of Figure 9 more closely matches the ratio observed by PAMELA above 1 GeV/nuc after taking into account reacceleration. The inclusion of positrons proved the greatest challenge even after also considering the effects of convection, indicating that GALPROP in general might not yet be optimally suited to compute positron LISs. For future studies antiproton observations by PAMELA (Adriani et al. 2010) and/or AMS-02 (Aguilar et al. 2016) and the ability of GALPROP to fully compute secondary antiparticles, need to be investigated in more detail. If follows from Figure 8 that the LIS for boron as calculated with GALPROP needs some refinement below 30 MeV/n when compared to V1 observations.

Isotopic CR observations from V1 have been studied extensively by Cummings et al. (2016) and such observations at the Earth from ACE by Lave et al. (2013). Isotopic observations made by PAMELA, such as deuterium, and the extensive range by the AMS-02 experiment could expand these previous studies and what is presented in this work. Heavier CR nuclei can also be considered for follow-up studies, as well as ratios such as p/He, ³He₂/⁴He₂, and ¹⁰Be₄/⁹Be₄.

The results presented here are valuable in addition to other LIS studies, such as those by Webber (2015a, 2015b), who did not directly consider solar modulation, thus not quantitatively taking into account the observations at the Earth at lower energies. Studies that account for modulation, such as that by Webber & Higbie (2009), use the force-field approximation, which does not include important modulation effects such as charge-sign dependence (particle drifts) and perpendicular diffusion. These force-field models typically overestimate proton spectra below about 1 GeV during A < 0 solar magnetic cycles, while underestimating the intensity during A > 0 cycles (see, e.g., Potgieter & Moraal 1985; Caballero-Lopez & Moraal 2004; Potgieter & Vos 2017). The GALPROP code has also been used by Cummings et al. (2016) to produce LISs compared to the same observations, but they too only apply a force-field approximation. Studies have also been presented that attempt to improve on the basic force-field approximation, such as shown by Corti et al. (2016). A more complicated modulation model has also been presented by Boschini et al. (2017), who, similarly to the work presented here, use the GALPROP model together with their HELMOD solar modulation model. While these various modeling studies have shown the ability to produce CR LISs, we believe our 3D solar modulation model to be more refined, and certainly to exceed the force-field approximation. Together with our larger scope of CR LISs, which also includes positrons, we believe the new LISs presented here will hopefully form an additional integral part of future CR and solar modulation studies.

It was announced that Voyager 2 crossed the heliopause in a different direction than Voyager 1 on 2018 November 5, so that we look forward to see new observational very LISs from this mission (press release: https://voyager.jpl.nasa.gov/news/details.php?article_id=112).

The authors 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 GALPROP WebRun service. D.B. and O.P.M.A. acknowledge the partial financial support from the postdoctoral program of the North-West University, South Africa. This work is dedicated to the memory of William R. Webber who passed away at the end of 2018 November.