The Modulation of Anomalous and Galactic Cosmic-Ray Oxygen over Successive Solar Cycle Minima

The 2009 solar minimum was characterized by record-setting galactic cosmic-ray (GCR) intensities. Mewaldt et al. (2010) reported GCR intensities, in the energy range ∼70–450 MeV nuc⁻¹, to be 20%–26% greater in late 2009 than in the previous 1997–1998 minimum as well as all the previous solar minima of the space age (1957–1997). Interestingly, the 2009 solar minimum was in an A < 0 magnetic polarity cycle, which is usually associated with lower GCR levels than a similar A > 0 cycle. This led Strauss & Potgieter (2014) and Moloto et al. (2018), among others, to postulate that even higher GCR intensities would be observed in the 2020 A > 0 solar minimum if similarly quiet solar conditions prevailed just due to the drift effect. This was confirmed by Fu et al. (2021) showing that GCRs intensities in 2020 were ∼6% higher than in 2009 and ∼25% higher than the previous A > 0 drift cycle in 1997 at energies of ∼50–500 MeV nuc⁻¹. Record-setting intensities were also reported in 2009 from ground-based neutron monitors (NMs), which are sensitive mostly to protons with energies >500 MeV (e.g., Moraal & Stoker 2010). This trend continued into 2020, as illustrated in Figure 1. Here, we show long-term NM observations, normalized to 100% in 2009, by the South African NM network, extending back to 1957. In this data set, the 2009 intensities showed an enhancement of ∼1% over that of the 1987 A < 0 cycle, while being only slightly higher than the 1965 A < 0 cycle. However, the recent 2020 intensities, in an A > 0 drift cycle, are ∼2% higher than previous A > 0 maxima in 1997 and 1975. It should be noted that, because of drift effects, high-energy GCR observations at NM energies are lower in an A > 0 cycle than an A < 0 cycle, in contrast to low-energy space-based measurements, leading to well-known crossover of the energy spectra (e.g., Reinecke & Potgieter 1994).

Zoom In Zoom Out Reset image size

Possible reasons for the very high GCR levels in the recent solar cycles have been debated over the past decade. Ultimately, the high GCR levels are due to extraordinarily low levels of modulation caused by the very quiet Sun in both 2009 and 2020. The top panel of Figure 2 shows, as a solar proxy, the sunspot number, which exhibits a downward trend over recent cycles, with a long and prolonged minimum in both 2009 and 2020. The main process influencing GCR modulation is still debated, but it is probably driven by a combination of factors, including a weaker and less turbulent magnetic field (Smith & Balogh 2008), as well as a slowly declining heliospheric current sheet (HCS) tilt angle (Wang et al. 2009).

Zoom In Zoom Out Reset image size

In general, the intensities of anomalous cosmic rays (ACRs) are observed to track their GCR counterparts rather well. At least, this was the situation for most of the Space Age. A striking feature of the 2009 solar minimum is that, although GCRs reached record-setting intensities, ACR intensities only reached 90% of their 1997 levels (Leske et al. 2009, 2013; McDonald et al. 2010). This trend continued in 2020 (Fu et al. 2021), as is also illustrated in Figure 1. A natural interpretation for both the record-setting GCR intensities and the fact that ACRs were lower in recent solar cycles was presented by Moraal & Stoker (2010): these authors argue that the prevailing quiet solar conditions in 2009 (and therefore also in 2020) led to less GCR modulation due to lower levels of turbulence that increase the diffusion coefficients (and therefore the mean free paths), because of less efficient particle scattering. For ACRs, however, a larger diffusion coefficient leads to less efficient particle acceleration at the solar wind termination shock. While ACRs are therefore also less modulated in 2009 and 2020, their source intensities may also be less, creating an intricate interplay between acceleration efficiency and modulation. This idea has not been robustly tested before, and we aim to do so in this current work by simulating the modulation of ACR and GCR oxygen, simultaneously in a numerical modulation model, for different levels of solar modulation.

For the simulations presented here, we use the 2D acceleration and transport model of Strauss et al. (2010a, 2010b, 2011), which solves the isotropic Parker (1965) transport equation in the meridional plane of the heliosphere. This model is based on earlier work by le Roux et al. (1996) and Langner et al. (2004) where the Parker transport equation is integrated numerically by a finite difference numerical scheme discussed in detail by Langner (2004). Diffusive shock acceleration (DSA) of a pick-up ion source at the solar wind termination shock (TS; assumed to be located at 90 au) is included, with a latitude-independent compression ratio of s = 2.2. Incompressible plasma flow is assumed in the heliosheath, leading to no adiabatic energy losses/gains in this region, while second-order Fermi acceleration (i.e., momentum diffusion) is also neglected. A solar wind speed of 400 km s⁻¹ is assumed at Earth in the equatorial plane, increasing to 800 km s⁻¹ at the polar regions. A Parker (1958) heliospheric magnetic field (HMF) is used in the equatorial regions, but it is modified in the polar regions in accordance with Jokipii & Kóta (1989). This field is assumed to also be frozen into the slowing heliosheath plasma. Gradient, curvature, and HCS drifts are included, and we use a constant HCS tilt angle of 10° for all simulations, as only solar minimum periods are considered.

The transport (diffusion and drift) coefficients are discussed in Section 2.1. It is important to note that the only differences between the GCR and ACR simulations presented here are the charge of the particle population under consideration (fully ionized GCRs and single-ionized ACRs) and the particle source (for ACRs, a pick-up ion source function is specified at the TS, while for GCRs, a local interstellar spectrum, LIS (see Section 2.2), is specified at the heliopause, HP, assumed to be located at 120 au).

2.1. Transport Parameters

The transport coefficients used in this study are shown, as a function of rigidity, in the left panel of Figure 3, and briefly discussed below. The parallel and perpendicular coefficients are also compared to the so-called Palmer consensus values (as taken from Bieber et al. 1994) in order to validate the assumptions used in this study. Note that the diffusion coefficient, κ, is related to the corresponding mean free path, λ, as κ = v/3λ, where v is the particle speed. The right panel of the figure relates the kinetic energy to the rigidity for the different oxygen species considered here. Note that, due to the different charge-to-mass ratio, this relationship is different for the different species.

Zoom In Zoom Out Reset image size

For the parallel mean free path, we use the results of Burger et al. (2008) and Engelbrecht & Burger (2010), based on the calculations of Teufel & Schlickeiser (2003):

$\begin{eqnarray}&&{\lambda }_{| | }=\zeta \displaystyle \frac{3}{\pi }\displaystyle \frac{{B}^{2}}{\delta {B}_{\mathrm{slab}}^{2}}\displaystyle \frac{{k}_{\min }{r}_{{\rm{L}}}^{2}}{\nu -1}\left[\displaystyle \frac{1}{4}+\displaystyle \frac{2{\left({k}_{\min }{r}_{{\rm{L}}}\right)}^{-\nu }}{(2-\nu )(4-\nu )}\right].\end{eqnarray} \tag{ 1 }$

In this expression, B is the HMF magnitude, δ B² is the variance associated with the slab component of the turbulent magnetic field, r_L is the maximal Larmor radius, k_min is the wavenumber associated with the transition from the energy to the inertial range of the slab power spectrum, and ν = 5/3 is the spectral index of the inertial (Kolmogorov) turbulence range. Following Burger et al. (2008), we use ${k}_{\min }=32\ {\mathrm{au}}^{-1}\sqrt{r}$, where r is radial distance. The value of k_min decreases until the TS is reached and kept constant in the heliosheath. Similarly, we assume that $\delta {B}_{\mathrm{slab}}^{2}=13.2\ \mathrm{nT}{r}^{-2.5}$ inside the TS but $\delta {B}_{\mathrm{slab}}^{2}\propto {B}^{2}$ farther into the heliosheath. The free parameter ζ is introduced to account in an ad hoc fashion for time-dependent changes in λ_∣∣ (see, e.g., the discussion by Strauss & Potgieter 2010) and is fixed by comparing model results to observations.

For the perpendicular diffusion coefficients, κ_{⊥,r θ} , we use the estimate of Burger et al. (2000):

$\begin{eqnarray}&&{\kappa }_{\perp /r}=\eta {\kappa }_{| | }{\left[\displaystyle \frac{P}{1\ \mathrm{GV}}\right]}^{-1/3}\ \mathrm{and}\ {\kappa }_{\perp ,\theta }=\eta g(\theta ){\kappa }_{| | }{\left[\displaystyle \frac{P}{1\ \mathrm{GV}}\right]}^{-1/3}.\end{eqnarray} \tag{ 2 }$

Note that the perpendicular diffusion coefficient in the polar region, κ_⊥,θ , is enhanced with respect to the radial direction, κ_⊥,r , via the function g(θ), leading to isotropic perpendicular diffusion in the equatorial plane but becoming increasingly anisotropic, so that g = 4 at the polar regions. This enhancement has been implemented in order to reproduce observed latitudinal gradients (e.g., Potgieter et al. 1997). The parameter η is adjusted in the following sections in order to reproduce 1 au observations.

The efficiency of drift effects can be characterized by the so-called drift scale, λ_d , which reduces to the Larmor radius if scattering is negligible (e.g., Minnie et al. 2007; Burger & Visser 2010; Engelbrecht et al. 2017; van den Berg et al. 2021). If, however, solar wind turbulence is present, λ_d = f_s r_L, with f_s ∈ [0, 1] referred to as the drift reduction factor. Here, we use a parameterized form of this reduction factor from Burger et al. (2000) and Langner et al. (2004),

$\begin{eqnarray}&&{f}_{s}=\chi \displaystyle \frac{{\left(P/{P}_{0}\right)}^{2}}{{\left(P/{P}_{0}\right)}^{2}+1},\end{eqnarray} \tag{ 3 }$

with P₀ = 1/10 GV. Again, χ will be fixed by comparing simulation results with observations.

In contrast to incompressible turbulence observed in the supersonic solar wind, turbulence in the heliosheath has been observed to be mostly compressible (e.g., Burlaga & Ness 2009). At present, a complete theory of particle scattering (both parallel and perpendicular to the mean field) in compressive turbulence does not exist, and as such, Equations (1) and (2) are assumed to be valid throughout the heliosphere and are also applied in the heliosheath.

2.2. The Oxygen LIS

For the GCR oxygen simulations, we assume an LIS that is comparable to the available observations while having a rather simplistic mathematical form,

$\begin{eqnarray}&&{j}_{\mathrm{LIS}}={j}_{0}{\beta }^{\gamma }{\left(\displaystyle \frac{E}{1\ \mathrm{GeV}\,{\mathrm{nuc}}^{-1}}\right)}^{-\delta },\end{eqnarray} \tag{ 4 }$

where β is the ratio of particle speed to the speed of light, E is the kinetic energy per nucleon, γ = 6.25, δ = 2.625, and j₀ = 6250 cm⁻² s⁻¹ sr⁻¹ MeV⁻¹ nuc. In Figure 4, this LIS is compared to other estimates by Boschini et al. (2018), Bisschoff et al. (2019), and Shen et al. (2019), where the latter should rather be interpreted as a source function, specified in their simulations at 80 au. Also shown in the figure are high-energy AMS-O2 data (Aguilar et al. 2017), taken between 2011 and 2016, and low-energy Voyager measurements beyond the HP (Cummings et al. 2016) used to constrain the LIS.

Zoom In Zoom Out Reset image size

We start by using GCR oxygen energy spectra at Earth to constrain the remaining free parameters in the transport model. To do so, we use the highest-energy channels of the ACE CRIS instrument (Stone et al. 1998a) from 2008 and 2020, along with recent PSP ISOIS (McComas et al. 2016) data from the EPI-Hi sensor in 2020 (Rankin et al. 2022). The model is then run for different values of the transport parameters until a reasonably good fit is obtained for both polarity cycles. We find optimal agreement with ζ = 1.08, η = 0.05, and χ = 0.6. Here, ζ > 1 suggests that both the 2009 and 2020 solar minimum were less active than previous cycles, and it is generally consistent with lower turbulence levels. A value of η = 0.05 for the ratio of perpendicular to parallel diffusion is about a factor of 2 higher than previous estimates (e.g., Ferreira & Potgieter 2004), but this is still reasonable considering how much this ratio depends on the level of turbulence (see, e.g., Giacalone & Jokipii 1999), which could be very different for the recent quiet solar cycles. A drift reduction factor of χ = 0.6 is in line with previous estimates of this quantity from, e.g., Langner et al. (2004), who found a value of 0.55 for this coefficient during solar minimum conditions.

With all free parameters now determined, we run the model for ACR oxygen in both polarity cycles. As the exact level of the ACR oxygen seed population at the TS is not known, the simulations must, unfortunately, be normalized to 1 au levels. Note that the same normalization factor is used in simulations in both polarity cycles. For the normalization, we use A > 0 measurements in 2020 from PSP ISOIS (EPI-Hi; McComas et al. 2016) and ACE SIS (Stone et al. 1998b). However, as shown by Rankin et al. (2022), ACR oxygen spectra for different solar cycles compare extremely well, allowing us to also use data from the 1990s, including SOHO measurements from Marsden et al. (1999) and SAMPEX measurements from Mazur et al. (2000). For A < 0 measurements, we use ACE and STEREO measurements (Cummings et al. 2009).

Our best-fit scenario, for both ACR and GCR oxygen, is shown in Figure 5. Here, A > 0 simulations and measurements are shown in red, while A < 0 simulations and measurements are shown in blue. The individual ACR and GCR simulations are added in gray. The simulations compare very well with both the ACR and GCR measurement for both polarity cycles. For the GCR Oxygen spectra presented here, we do not observe the crossover for different drift cycles alluded to earlier. While this can potentially occur at higher energies than considered here, its absence may point to an inaccurate form of our assumed transport coefficients at the highest energies, or it could be that such a crossover does not exist for oxygen. However, such high-energy discrepancies will not influence the rest of the results presented here.

Zoom In Zoom Out Reset image size

To show the effect of changing diffusion conditions on the resulting ACR and GCR spectra at 1 au, we vary the value of ζ by ±50% and show the resulting spectra in Figure 6. Here, the top panels correspond to GCR simulations and the bottom panels to ACR simulations. The left panels are for the A > 0 polarity cycle, while the right panels are for the A < 0 polarity cycle. The result for GCR oxygen is rather straightforward; a larger diffusion coefficient, which corresponds to less turbulent conditions, leads to higher GCR levels with slight differences between different drift cycles. Most notably, the GCR intensity is much more sensitive to changes in the assumed turbulence conditions during A < 0 cycles than during A > 0 cycles.

Zoom In Zoom Out Reset image size

For ACR oxygen, however, an interesting interplay is obtained between the ACR source spectrum at the TS and the resulting modulated spectrum at Earth. At the TS, a smaller diffusion coefficient, which corresponds to more turbulent conditions, leads to an accelerated spectrum extending to higher energies. This is a well-known result for DSA (see, e.g., Steenberg & Moraal 1999; Moraal & Stoker 2010; Prinsloo et al. 2019); a smaller diffusion coefficient leads to better particle confinement at/near the shock, allowing the DSA process to continue to higher energies before an exponential cutoff is observed. In general, therefore, more turbulence (i.e., more particle scattering) leads to more efficient ACR acceleration at the TS. However, this resulting harder energy spectrum will now be modulated more efficiently by the additional turbulence as the ACR particles propagate toward Earth. In the bottom panels of Figure 6, this interplay between more efficient acceleration versus more effective modulation is illustrated for both drift cycles. More turbulent conditions (i.e., a smaller diffusion coefficient) lead to higher ACR intensities at the highest ACR energies (≳5 MeV nuc⁻¹) at both the TS and at Earth. For low-energy ACRs at Earth, however, a smaller diffusion coefficient leads to lower ACR intensities.

We can now use the preceding simulations to address the observations by Leske et al. (2013), who found higher GCR but lower ACR intensities during the recent, quiet solar minima. The simulations presented in Figure 5 are representative of the quiet solar conditions that were observed in both 2009 and 2020, and we refer to these as results for recent solar cycles. Next, we proposed that early solar cycles, i.e., the solar cycles before 2009, were generally more active (i.e., more turbulent) than these recent cycles, and we assume that the transport coefficients (both the diffusion and drift coefficients) were 25% smaller in early cycles as compared to recent solar cycles. The resulting simulations are shown in the top panel of Figure 7 for the A > 0 (left panels) and A < 0 (left panels) polarity cycles. Results for recent solar cycles are shown in red, while blue corresponds to results for early solar cycles with a 25% reduction in all transport coefficients. As expected from the previous simulations, this results in GCR intensities that are higher in recent solar cycles but ACR intensities for E ≳ 5 MeV being higher in earlier cycles. At the lowest energies, we again find higher intensities in recent cycles.

Zoom In Zoom Out Reset image size

This behavior is best illustrated in the bottom panels of Figure 7, where we show the ratio of modeled intensities for recent/early cycles. We do so for ACRs (black dashed line) and GCRs (dashed–dotted line) separately but also for the combined ratio (solid red line). For the sake of completeness, we also repeat the calculation assuming only a 10% reduction in the transport coefficients, and we show this as the red dotted line. Also shown, as the green curve, is the measured ratio from ACE (ratio of 2018–2020 to 1997–1998), showing a very good qualitative comparison with the model results.

We show modeling results where the same set of transport coefficients is used to reproduce combined ACR and GCR measurements, at Earth, for two recent solar minima. As the GCR oxygen LIS has now been directly observed by the Voyager spacecraft (Cummings et al. 2016; Stone et al. 2019), a comparison between GCR simulations and observations allows us to constrain the remaining free parameters in the modulation model. The model is then used to also simulate ACR spectra, which we find to be consistent with recent observations.

We study the effect of changing diffusion conditions on energy spectra and find that GCRs behave as expected, with larger diffusion coefficients, associated with more quiet solar conditions, leading to higher GCR fluxes during recent solar cycle minima. For ACRs, this is more complex, as the acceleration of ACRs at the TS is also affected by the magnitude of the diffusion coefficient, with more turbulent conditions leading to more efficient ACR acceleration, while the ACR will be more effectively modulated while propagating toward Earth. Our simulations suggest that high-energy ACR intensities will therefore be lower during very quiet solar minimum conditions. At the lowest ACR energies, however, ACR intensities are again higher, illustrating the interplay between acceleration and modulation and their opposite depending on turbulence (scattering) conditions.

By assuming that recent solar cycles are associated with less scattering, we calculate the combined ACR and GCR spectral ratio between recent and early solar cycles. Our results are in agreement with observations presented by Leske et al. (2013) showing higher GCR and lower ACR fluxes during the recent quiet conditions. This is also in line with the proposed explanation by Moraal & Stoker (2010), which is now robustly tested with simulations.

This work is based on research supported in part by the National Research Foundation of South Africa (NRF grant Nos. 119424, 120345, and 120847). Opinions expressed and conclusions arrived at are those of the authors and are not necessarily to be attributed to the NRF. The responsibility of the contents of this work is with the authors. R.A.L. acknowledges support from NASA grant 80NSSC18K0223. This work was also supported as a part of the Integrated Science Investigation of the Sun on NASA's Parker Solar Probe mission, under contract NNN06AA01C. Figures were prepared with Matplotlib (Hunter 2007) and certain calculations were done with NumPy (Harris et al. 2020). Sunspot number measurements were provided courtesy of the Royal Observatory of Belgium, Brussels. HCS tilt angle and solar polar magnetic field measurements were provided courtesy of the Wilcox Solar Observatory (WSO). WSO data used in this study were obtained via the website http://wso.stanford.edu at 2022:11:04_18:40:53 PDT, courtesy of J. T. Hoeksema. The WSO is currently supported by NASA.