Long- and Short-term Variability of Galactic Cosmic-Ray Radial Intensity Gradients between 1 and 9.5 au: Observations by Cassini, BESS, BESS-Polar, PAMELA, and AMS-02

LEMMS was a double-ended energetic charged particle telescope that belonged to Cassini's Magnetospheric Imaging Instrument (MIMI) suite (Krimigis et al. 2004). The signal of GCRs in LEMMS, which appears as low-intensity background noise in many of its 57 channels, can be easily isolated. We refer the reader to Roussos et al. (2018a, 2019) for the details on how this signal is extracted, and we focus on the data processing aspects that are most relevant or novel for the present investigation. Additional information is provided in the Appendix.

A new element compared to the past, the Cassini-based GCR studies that we cite above, is that we extended our data set back to the first switch-on of the LEMMS instrument on 1999 January 3. The main purpose was to analyze the LEMMS GCR data around the time of the Earth flyby, compare them with well-defined measurements of GCR proton spectra at 1 au (Section 3.2), and obtain an absolute GCR flux calibration for LEMMS.

Here we use data from LEMMS channel E7. Away from sources of planetary electrons and protons in Saturn's radiation belts, E7 is responding to >300 MeV protons. The signal-to-noise ratio easily exceeds 50 when few-days averaging is applied to its GCR time series. In addition, the E7 calibration settings were identical at 1 au and for the Saturn tour, while the channel shows no sensitivity to neutrons and gamma rays from Cassini's Radioisotope Thermal Generators (RTGs). At the energies of E7, protons are dominated by GCRs for periods with quiet solar wind. Solar proton fluxes become comparable or stronger only at energies below ∼10–20 MeV.

The channel's most important property is that contributions to its signal from sideways-penetrating GCR protons is the lowest compared to all LEMMS proton channels. Sideways penetrators may propagate through the body of the Cassini spacecraft before reaching LEMMS, an effect not taken into account in our channel-response function simulations. In that respect, we consider that the modeled responses of E7 (Figure 1) are the most reliable compared to other LEMMS channels. The channel's directionality and capability to reject sideways penetrators have been confirmed through measurements of the harsh environments of Earth's and Saturn's proton radiation belts (Roussos et al. 2018b). The processing steps for E7 data were the following:

• 1.  
  We excluded time periods through Earth's and Saturn's radiation belts, as well as data for an extended period around the Cassini flyby of Jupiter (day 360/2000–day 060/2001) where Jovian relativistic electrons (>7 MeV) were detected with E7 in and out of the planet's magnetosphere (Krupp et al. 2004).
• 2.  
  We excluded time periods of interplanetary coronal mass ejections (ICMEs). Even though the high-energy proton threshold at 300 MeV makes E7 insensitive to the direct detection of ICME protons, its GCR signal captures the associated Forbush decreases (Lockwood 1971). Lists of ICMEs were obtained from Roussos et al. (2018a, 2018c) and Lario et al. (2004, 2005). As these three studies did not fully cover 1999, 2003, and early 2004, we made an additional survey to add missing events.
• 3.  
  We excluded time periods where the calibration settings for E7 were different compared to those at 1 au and for the Saturn tour. This amounts to several months of data in 2001 (Figure 10 in the Appendix).
• 4.  
  After all these filtering steps, we applied a 26 day average to the E7 count rates in order to suppress solar periodicities induced by corotating interaction regions (CIRs; Roussos et al. 2018c). Even though shorter averaging intervals were possible, it is expected that CIR-driven oscillations at 1 au and at Cassini would be out of phase, leading to unphysical solar periodicities of G_r through application of Equation (1).

Zoom In Zoom Out Reset image size

The resulting data, along with relevant ephemeris information, are shown in Figure 2.

The measurements by Cassini are compared with GCR observations at Earth (to which we may refer by Earth's heliocentric distance at 1 au). For that purpose, we use data obtained from several advanced GCR observatories (Table 1), namely BESS, BESS-Polar, PAMELA, and AMS-02. All these observatories comprise large magnetic spectrometers capable of retrieving precision differential energy spectra of GCRs with high energy and species resolution throughout the energy range of Cassini's >300 MeV proton measurements. These energy-resolved spectra offer several advantages for our analysis compared to the use of integral flux measurements (e.g., van Allen 1988) or compared to spectral reconstructions based on neutron monitors and theory, as will be presented in Section 3.3.

Among the different observatories, BESS and BESS-Polar were high-altitude balloons operated for 12 short-duration missions between 1993 and 2008. Only the five missions after 1999 (1999, 2000, 2002, 2004, 2008), when Cassini's MIMI/LEMMS was functional, are relevant here. PAMELA and AMS-02 operated or still operate on long-duration orbital missions, offering time series of GCR spectra nearly continuously since 2006. Certain measurements by BESS-Polar, PAMELA, and AMS-02 overlap in time such that small differences between the respective spectra can be evaluated and cross-calibrated. The BESS flight in 1999 overlaps both in space and time with Cassini's Earth flyby, as mentioned in Section 2.

GCR spectra from the aforementioned experiments were subject only to minor processing described below:

• 1.  
  All GCR spectra were interpolated at common time and energy tags, such that simultaneous measurements can be compared and cross-calibrated, when available. Linear interpolation was used for GCR time series at a fixed energy, while logarithmic interpolation was applied for GCR fluxes as a function of energy.
• 2.  
  AMS-02 spectra where extrapolated from their minimum energy at 430 MeV (Table 1) down to the lowest energy where LEMMS E7 responds (∼300 MeV) as this is necessary for reconstructing the LEMMS count rate in one of our analysis steps introduced in Section 4.2. Extrapolation was performed through an analytical spectrum function (Section 3.3), requiring that its flux at 430 MeV matches the one measured by AMS-02. This extrapolation has negligible impact on our results because of the relatively narrow energy range involved and the low GCR fluxes between 300 and 430 MeV and is primarily done for the application of several numerical integrations starting from the ∼300 MeV threshold.
• 3.  
  PAMELA fluxes were divided by an energy-dependent scaling factor, P(T) (T: energy). This factor was estimated by time-averaging the energy-dependent ratio of simultaneous PAMELA and AMS-02 spectra after 2011. For the energy range of interest, that ratio ranged between 0.95 and 1.1, indicating good agreement between the two data sets even without correction. For most energies, the PAMELA fluxes were slightly higher than that of AMS-02 (P(T) > 1). The reason we adjusted PAMELA fluxes to match AMS-02 and not the other way around is that PAMELA fluxes were similarly larger than those estimated from the second BESS-Polar flight in 2008. This choice on how to match PAMELA and AMS-02 fluxes has some minimal impact on our final results, which will be discussed in Sections 5.1 and 6.
• 4.  
  Since proton spectra from BESS are provided until 20 GeV (Shikaze et al. 2007), we extended them manually to 100 GeV. To do that, we estimated the spectral index from measurements between 15.0 and 20.0 GeV for each of the BESS spectra, and we used that to extrapolate them as a power law above that energy. This extrapolation has a negligible impact on our results and was only performed such that several numerical integrations described in Section 4 can be performed. This procedure does not concern spectra from BESS-Polar that are available up to 156 GeV (Abe et al. 2016).

Example time series of 0.5 and 2.0 GeV proton fluxes are shown in Figure 3. By adjusting slightly only the fluxes of PAMELA, simultaneous measurements by the three experiments are in excellent agreement, with differences in fluxes rarely exceeding 1%.

Zoom In Zoom Out Reset image size

While our analysis relies on measured GCR fluxes, an analytical approximation for the GCR proton spectra can provide useful, practical context. This context information is important for assessing aspects of LEMMS's calibration or results from combined Cassini and 1 au measurements before 2006, which are only available for four isolated intervals (Figure 3). They also serve to complete the missing spectral information of AMS-02 below 430 MeV (Section 3.2).

The reconstruction that we employ here is based on the force-field approximation (Gleeson & Axford 1968), which offers a sufficiently good prediction of the GCR spectral shape and fluxes above 300 MeV that it is applicable for the heliospheric distance range of 1–10 au and serves well the purposes for which we intend to use reconstructed spectra in the present study (Usoskin et al. 2017; Roussos et al. 2019). The GCR differential proton flux spectrum, j(t, T), at a given time (t) and kinetic energy (T) is controlled through a single free parameter, the so-called modulation potential (ϕ[GV]). In the case of protons, it can be represented in units of energy as Φ[GeV] = e ϕ, with e the electron charge. The regulation of this parameter by the solar activity is the way in which force-field spectra are modulated. The equation of the GCR spectrum is as follows:

$$\begin{eqnarray}&&j(t,T)={j}_{\mathrm{LISM}}(T+{\rm{\Phi }}(t))\displaystyle \frac{{E}^{2}-{T}_{o}^{{}^{2}}}{{\left(E+{\rm{\Phi }}(t)\right)}^{2}-{T}_{o}^{2}},\end{eqnarray} \tag{ 2 }$$

with E = T + T_o[GeV] being the total proton energy and T_o = 0.938 GeV the proton rest mass. The spectrum of the local interstellar medium (j_LISM) in Equation (2) is given by

$$\begin{eqnarray}&&{j}_{\mathrm{LISM}}=2.7\cdot {10}^{3}\displaystyle \frac{{T}^{1.12}}{{\beta }^{2}}{\left(\displaystyle \frac{T+0.67}{1.67}\right)}^{-3.93}\ {\left({{\rm{m}}}^{2}{\rm{s}}\mathrm{sr}\mathrm{GeV}\right)}^{-1}\end{eqnarray} \tag{ 3 }$$

with β = v/c the ratio of the proton speed to the speed of light. Time series of ϕ have been calculated at 1 au by Usoskin et al. (2017) and are regularly updated at http://cosmicrays.oulu.fi/phi/phi.html. The calculation of ϕ is obtained through neutron monitor observations that have been cross-calibrated with PAMELA (Usoskin et al. 2017).

While the force-field approximation achieves a good, first-order representation of the GCR fluxes in the inner heliosphere, it also shows nonnegligible deviations from the measured spectra with amplitudes that depend on the solar-cycle phase. These deviations, which partly stem from the approximation's incomplete physics (Caballero-Lopez & Moraal 2004; Moraal 2013), peak around solar maximum, with fluxes that can be overestimated by as much as a factor of two (Gieseler et al. 2017; Corti 2019). This time-dependent deviation may translate into an artificial solar-cycle dependence of the radial gradient, which is why this analytical reconstruction is mostly used here for context, as explained in the beginning of this subsection.

In this method, we try to find if LEMMS's channel E7 can be described by a characteristic effective energy (T_eff), for which we can convert its count rate to differential flux. In that case, the gradient can be directly estimated through Equation (1), since fluxes at 1 au are energy-resolved. We describe this approach below.

The simplest method for converting channel count rates (R) to differential fluxes (j) is through Equation (4):

$$\begin{eqnarray}&&j=\displaystyle \frac{R}{G\,{\rm{\Delta }}T},\end{eqnarray} \tag{ 4 }$$

where G is the average, omnidirectional geometry factor within the channel's energy interval ${\rm{\Delta }}T={T}_{\max }-{T}_{\min }$. The characteristic energy for j is at the geometric mean of the channel boundaries, that is, at $\sqrt{{T}_{\max }\,{T}_{\min }}$ (Kronberg & Daly 2013). Since the response of LEMMS's E7 channel is integral in energy, this approximation cannot be used. Instead, because for a differential energy spectrum, j(T), the count rate can be approximated as

$$\begin{eqnarray}&&R={\int }_{{T}_{\min }}^{{T}_{\max }}j(T)G(T){dT},\end{eqnarray} \tag{ 5 }$$

then from Equation (4) we obtain the count rate to flux conversion factor, g:

$$\begin{eqnarray}&&g=G\,{\rm{\Delta }}T=\displaystyle \frac{{\int }_{{T}_{\min }}^{{T}_{\max }}j(T)G(T){dT}}{j(T)}.\end{eqnarray} \tag{ 6 }$$

In order to evaluate g, we need to use a predefined spectral form for j(T). This form is given through Equation (2), which is a valid approximation of the GCR spectrum between 1 and 9.5 au. For the E7 geometry factor G(T), we use the profile of Figure 1, integrated over the polar angle (see Figure 9 in the Appendix). Integration limits (${T}_{\min },{T}_{\max }$) are between 100 MeV and 100 GeV. A limitation is that for all of Cassini's positions besides its crossing from 1 au, Φ is unknown. We therefore evaluate the energy profile of g(T) for different, realistic values of Φ. The results are shown in Figure 4 (left). The same calculation can be performed using the spectra measured by BESS, PAMELA, and AMS-02. These spectra are detailed enough that the shape of j(T) can be resolved through measurements only. Here we obtain different spectra by sampling them from different time periods. Evaluation of Equation (6) for this case is shown in Figure 4 (right). We stress that for both cases, estimation of g(T) depends only on the spectral shape, not on the absolute GCR fluxes. Assuming that the spectral shapes within 9.5 au are similar, we may rely on 1 au measurements and approximations to estimate g(T) for any time period of the Cassini mission.

Zoom In Zoom Out Reset image size

For both estimations, we observe that around 2.0–2.2 GeV, all g(T) curves converge toward a single energy and g values that are independent of the chosen Φ (left) or date (right). An effective energy and flux conversion factor for channel E7 (T_eff, g_eff) can be evaluated by finding the centroid of all these curve intersections. Based on the analysis of the measured spectra (Figure 4, right), we obtain that for GCR-type spectra, channel E7 can be described with T_eff = 2.0 ± 0.2 GeV and g_eff = (1.0 ± 0.1)10⁻⁴ m² sr GeV. The method of determining these two parameters is a generalization of the "bow-tie analysis" (Van Allen et al. 1974), typically applied with the use of power-law spectra where the variable input is that of the spectral index.

We can now estimate the differential flux at 2.0 GeV measured by Cassini, requiring also that its value should agree with the 2.0 GeV flux measured by the BESS99 experiment, which occurred close in time to the Cassini flyby. In that case, we obtain

$$\begin{eqnarray}&&{j}_{\mathrm{Cassini}}(2.0\,\mathrm{GeV})=\displaystyle \frac{R}{4.086\,{g}_{\mathrm{eff}}},\end{eqnarray} \tag{ 7 }$$

with the factor 4.086 accounting for the aforementioned normalization at 1 au. For this normalization, we compared Cassini data between 0.75 and 1.25 au with those of BESS99. The magnitude of this correction is discussed in the Appendix. We note here that BESS or BESS-Polar measurements require an adjustment that accounts for a small reduction of GCR fluxes at the operational altitude of these spectrometers within Earth's upper atmosphere (Sanuki et al. 2000; Shikaze et al. 2007). This kind of adjustment may introduce some systematic errors in the obtained spectrum. It is thus important that the proton spectrum from BESS99 is in excellent agreement with fluxes from neutron-monitor-based reconstructions or other GCR flux proxies (Gieseler et al. 2017; Usoskin et al. 2017). We therefore do not find a reason to apply any additional correction factor to the BESS99 fluxes, before scaling Cassini's measurements to them.

Applying Equation (7) for Cassini (j_Cassini) and retrieving from the measurements the 2.0 GeV fluxes at 1 au (j_Earth), we obtain the time series shown in the top panel of Figure 5. Following 2001, j_Cassini is steadily higher than j_Earth, as expected for a positive GCR radial intensity gradient. Unfortunately, because GCR spectra at 1 au are not continuously available before 2006, we cannot utilize all of the Cassini observations for our analysis, which are available also for 1999–2006 (Figure 2).

Zoom In Zoom Out Reset image size

A careful inspection of the j_Cassini time series reveals small differences from the count rate profile of Figure 2, even though the two quantities, based on Equation (7), should be proportional. That is because the profile of j_Cassini has been time-shifted to account for solar wind propagation effects between 1 au and Cassini and then interpolated back to the time tags of j_Earth, such that Equation (1) can be applied. The time shift was performed assuming an average solar wind velocity of 431 km s⁻¹ (Echer 2019; Roussos et al. 2019),

With the data shown in Figure 5 (top), we can apply Equation (1), where j₁ = j_Earth and j₂ = j_Cassini (Section 4.3). A possible caveat of this method is that the resulting profiles depend significantly on the value of the effective energy, T_eff. This is a calibration-sensitive value that affects the time series of j_Earth that we choose to compare with, since T_eff controls the amplitude of their sinusoidal solar-cycle variation of GCR fluxes that decreases with energy. For this reason, in the following subsection we also apply a method that does not depend on the choice of an effective energy.

The second method is driven by the following question: what would have been the count rate of LEMMS's channel E7 had Cassini remained continuously at 1 au? As proton spectra at 1 au are known, we can reconstruct time series of this count rate (R_Earth) by folding them into Equation (5). At Cassini's position, the count rate (R_Cassini) is the actual one measured by channel E7 (Figure 2). Similarly, the simulated values of R_Earth were normalized so that they match the ones measured by Cassini around its Earth flyby. Overall, we obtain

$$\begin{eqnarray}&&{R}_{\mathrm{Earth}}=4.116{\int }_{{T}_{\min }}^{{T}_{\max }}j(T)G(T){dT},\end{eqnarray} \tag{ 8 }$$

with 4.116 being the normalization factor. Notice that the normalization factor is similar to that of Method 1 (Equation (7)), suggesting that the two methods are complementary. The resulting time series of E7 count rates are plotted in Figure 5 (bottom).

In principle, Equation (1) can be applied even if we set j₁ ≡ R_Earth and j₂ ≡ R_Cassini, assuming that the energy dependence of G_r in the range that E7 receives most of its signal from is secondary. This energy range is rather broad, as we illustrate in Figure 6. The GCR contribution to the E7 signal peaks at ∼0.6 GeV (black curve), but other energies are also important. If we apply Equation (5) but with a variable minimum energy integration limit, we obtain the cumulative count rate (Figure 6, bottom, red). This shows that 50% of the total count rate is contained below T = T_eff (Section 4.1), while 90% of the E7 signal comes approximately between 0.5 and 10 GeV. Gieseler & Heber (2016) and de Simone et al. (2011) indicate a slightly decreasing G_r between 0.43 and 1.2 GeV. At higher energies, the models predict that the decrease in G_r is stronger. A decrease in G_r at high GCR energies can be deduced directly from observations, since above ∼10–20 GeV, heliospheric GCR fluxes are similar as in the LISM. Alternatively, we may consider G_r in this case to be an integral gradient, as in many past studies (van Allen 1988; Honig et al. 2019).

Zoom In Zoom Out Reset image size

The calculated radial intensity gradient time series using Methods 1 and 2 are shown on the top panel of Figure 7. Time series of averaged G_r values over 4 month intervals and context information from past studies are also shown (de Simone et al. 2011; Gieseler & Heber 2016; Vos & Potgieter 2016; Honig et al. 2019). The monthly sunspot number (http://www.sidc.be/silso/monthlyssnplot), the heliospheric current sheet (HCS) tilt (http://wso.stanford.edu/Tilts.html), and the interplanetary magnetic field magnitude at 1 au from the Advanced Composition Explorer are added for context in the lower panels. We now focus on the the top panel of Figure 7. We may identify two extended time periods, before and after 2006.

Zoom In Zoom Out Reset image size

4.3.1. Radial Gradient: 1999–2006

4.3.2. Radial Gradient: 2006–2017

The radial gradient after 2006, when PAMELA and AMS-02 spectra are available almost continuously, displays several distinct features. Some of these features may be more easily visible in Figure 8, which zooms into the 2006–2017 period. Shaded areas in Figure 8 are overplotted to guide the eye for identifying short-duration structuring of the G_r profile and several of the following discussion points. We recommend that readers use Figures 7 and 8 in a complementary way.

Zoom In Zoom Out Reset image size

Since both calculation methods give nearly identical results for G_r, we show the long-term average profile from the second method only (Figure 7, top, red line, or Figure 8). Between 2006 and 2014, the gradient is quasi-constant with an average of G_r = 3.5 ± 0.3%/au, with evidence of a broad, weak minimum present between 2008 and 2011. After 2014, G_r experiences a clear, steady drop, reaching 2.0%/au by 2017 September. The dropout begins to develop around the time of the maximum of solar cycle 24 and the full transition to A > 0. Superimposed on this long-term 11 yr profile are hints of faster variations, each lasting between a few hundred days and ∼2 yr. In the following section, we expand on the validation of these results, based on which we proceed with initial physical interpretations of our findings.

LEMMS was not designed to measure GCRs, which is why we have invested a significant fraction of this study into the absolute GCR flux calibration of LEMMS (see also the Appendix), from the processing of the data and to the validation of our results. Here, we extend the validation efforts by comparing our findings with several independent G_r estimations. Some of these estimates and their uncertainties (de Simone et al. 2011; Gieseler & Heber 2016; Honig et al. 2019) are overplotted as shaded or hatched boxes in the top panel of Figure 7.

Most of these earlier studies provide G_r measurements between 2006 and 2009 July, with Honig et al. (2019) covering several more intervals in 2005, 2009–2011, and 2014. A comparison with the 2006–2011 cluster shows that our G_r values are similar, although they are generally toward the high end of the uncertainty range of those earlier studies. Given the multiple data processing steps and the diversity of the data sets compared, we consider such differences insignificant. A small source of discrepancy, for instance, may be the energy dependence of G_r. The values shown for de Simone et al. (2011) and Gieseler & Heber (2016) are for ∼1.2 GeV, the highest proton energies resolved in these two investigations, while our GCR measurements are more representative for T_eff = 2.0 GeV protons (Section 4.1).

In terms of data analysis pipeline, Table 2 summarizes how sensitive our results can be if we treat the cross-calibration between PAMELA and AMS-02 in two different ways (last two rows of Table 2). We find that in both cases, gross features of the G_r variability profile shown in Figures 7 and 8 are retained (e.g., fast drop of G_r after 2014), excluding the transient G_r minimum around the 2009 solar minimum, which cannot be resolved if correction factors to PAMELA and AMS-02 fluxes are obtained by averaging the two data sets at their overlapping time periods (2011–2014). Furthermore, all changes discussed in Table 2 lead to an average drop of the G_r amplitudes by 0.2%–0.5%/au, improving agreement with certain independent estimates, particularly for 2006–2008. We nevertheless consider the differences between the three cases small and the profile plotted in Figures 7 and 8 as representative. Any subsequent discussion points to the profile of those two figures.

The biggest discrepancies in G_r are with the study of Honig et al. (2019), who derive much lower gradients for 2005–2006 July and for 2014 January to March. While for the 2005–2006 period we do not have simultaneous estimates, it would seem unlikely that the gradient would change so drastically, from ∼1.6 to 3.4%/au within a few months after 2006 July. We note that the aforementioned calculation of Honig et al. (2019) took place 0.43–0.75 au away from Earth. It is possible that local solar wind dynamics may partly mask global-scale gradients of GCRs over such scales, leading to the small G_r values observed. Other estimates by Honig et al. (2019), farther away from Earth (2008 March–2011 July), are in much better agreement with our derivation.

The same argument cannot explain the even larger discrepancy in 2014. Honig et al. (2019) found that for the period 2014 January to March the gradient between 1 and ∼4.36 au was G_r = 2.13%/au, as opposed to 4.05%/au from our analysis. This measurement by Honig et al. (2019) is based on Rosetta observations from the vicinity of its target comet 67P. After 2014 March, GCR rates measured by Rosetta progressively dropped until they became smaller than those at 1 au, leading to a negative radial gradient and hinting that the comet locally modifies GCR fluxes by some process that is yet to be identified. Even though Honig et al. (2019) considers that before 2014 March the GCR fluxes may be unaffected by the comet, in view of our findings, we advocate that a partial GCR flux reduction could have been present at Rosetta also in that earlier period. What also favors such an interpretation is that Honig et al. (2019) presents GCR proxy data from Mars Odyssey at 1.5 au, indicating a signal comparable to that of Rosetta at 4.36 au before 2014 March. That would imply a considerable radial gradient between 1 au and the orbit of Mars and a vanishing one between Mars and 4.36 au, which would be surprising. We therefore consider that our estimate is more representative of the global radial gradient for 2014. We cannot exclude that additional dependencies, such as a heliocentric distance dependence of G_r, contribute to the differences, but we consider this unlikely, because if that was true, Cassini and Rosetta estimates would have diverged also in earlier time periods (2007–2011).

Based on all these arguments, we are confident that our G_r estimations after 2006 July (Figure 7, top panel) are valid within 0.5%/au. The comparison with independent G_r calculations and possible variations of our data analysis steps indicate that a slight G_r overestimation from our side is more likely than an underestimation. Also, GCR transport models generally predict slightly lower G_r values for few-gigaelectronvolt protons compared to what we find here (Vos & Potgieter 2016). But since our analysis relies on the combination of measurements at a quasi-constant heliocentric distance, with a single data set at 9.5 au and three accurately cross-calibrated data sets at 1 au (Figure 3), we consider that there is little ambiguity regarding the relative temporal G_r variations that we derive compared to past studies. The discussion that follows focuses on the resulting G_r variability pattern.

We now turn our attention to the G_r variability trends (Section 4.3). Most of our G_r observations are in solar cycle 24, such that we cannot make extensive comparisons with solar cycle 23. For this study it is more sensible to separate the G_r time series into their A < 0 and A > 0 phases as well as compare solar minimum against solar maximum observations.

A prediction of drift-modulated GCR transport models for the inner heliosphere is that G_r values during the A < 0 phase are larger than those of A > 0 (Burger & Hattingh 1998; Potgieter et al. 1989; Potgieter 2017). This theoretical prediction is in agreement with our observations, where we find that G_r progressively decreases into the A > 0 phase, reaching the lowest estimated values in 2017. This is even more obvious when comparing periods "A" and "F": at similar solar-cycle phases and GCR flux levels (Figure 8), their G_r differs by 1.5%/au, with the highest values for the A < 0 period in 2006–2007. In terms of solar-cycle phase control, our observations also indicate a stronger control for A > 0, when the G_r dropout develops over several years and in line with changes in the various solar indices shown in the lower three panels of Figure 7. For the A < 0 phase, the apparent alignment of the broad G_r minimum with the solar minimum also suggests a solar-cycle control, albeit a much weaker one than for A > 0. Correlation coefficients between G_r and the solar indices are generally positive and highest with the sunspot number and the HCS tilt (values ∼0.5). Due to the long, 4 month averaging used and the limited G_r data set, introducing time lags for exploring higher correlation coefficients does not yield clear results.

Our observation of a peaked and a transient G_r minimum toward solar maximum and solar minimum, respectively, is consistent with several G_r calculations for previous solar cycles (van Allen 1988; Webber & Lockwood 1991). Such a solar-cycle control may have a straightforward explanation: assuming a static GCR LISM spectrum, the difference between LISM fluxes and those within the heliosphere are the greatest during solar maximum, when heliospheric GCR fluxes reach their minimum levels (or the contrary). A difference with those earlier studies, however, is that they report similar G_r values on each side of the solar maximum, indicating a stronger control by solar-cycle phase (e.g., sunspot number), rather than from the polarity of the solar magnetic field. We cannot assess whether these disagreements signify issues in certain observations, differences between (odd/even) solar cycles, or a combination of the above.

As briefly mentioned in Section 4.3.2, the G_r time series display variations also on yearly or biennial timescales, most clearly visible in Figure 8. Several G_r changes before 2015 are apparent as enhancements bounded by clear minima (marked as A–D). After 2015, enhancements E and F are more subtle but still visible. The enhancements have a visible correspondence in GCR fluxes for events A, B, E, and F. Events C and D are more obvious in lower energy GCR measurements by LEMMS and at 1 au (Figure 3, Roussos et al. 2019).

Pizzella (2018) observed that many such features in GCR intensities at 1 au tend to recur every several hundred days. Since a similar recurrence was seen at 9.5 au, Roussos et al. (2019) argued that these are heliospheric modulations corresponding to the so-called quasi-biennial oscillations (QBOs). These are quasi-periodic variations of the solar magnetic field with a broad range of periods in the 0.6–2.0 yr range (Bazilevskaya et al. 2014) that propagate in the heliosphere and may modulate accordingly the GCR spectra (Krivova & Solanki 2002; Obridko & Shelting 2007; Mandal et al. 2017). Uncertainties and the long-averaging of our G_r time series prevent us from exploiting the full QBO periodicity range with accuracy. A Lomb–Scargle periodicity analysis (not shown here) leads to a noisy spectrum with several peaks in the expected range for QBOs. The durations of individual G_r enhancements and their alignment with corresponding features in GCR intensities are what gives some credence to the link between QBOs and the short-term variability of G_r.

There are several ways that QBOs may leave signatures in radial GCR gradients. One possibility is that their appearance is an artifact that is due to the fact that for the alignment of GCR observations at 1.0 and 9.5 au we are using an average solar wind velocity of 431 km s⁻¹ (Echer 2019; Roussos et al. 2019) for the whole time period analyzed. Large deviations from this velocity may lead to a misalignment of QBO events in GCR intensity and thus lead to an artificial gradient roughly at the time of each event's occurrence. We do not believe this to be the case here, since uncertainties in solar wind propagation times between 1 and 9.5 au may be as large as a week but only for rare extreme events. Propagation uncertainties are also small compared to our 26 day averaging of the GCR intensities, the 3–4 month averaging of the G_r series, and the duration of the QBO events.

For instance, the peak of event "A" in GCR fluxes at 9.5 au (Figure 8, top) is misaligned by several months to the corresponding 1 au feature (Figure 5). This misalignment is too large to be attributed to an error in solar wind propagation. The much better aligned event "B" features much sharper changes in GCR fluxes at 9.5 au, explaining the similarly sharp decrease in G_r at its outset in early 2010. These two examples suggest that QBO properties (spatial extent, intensity) evolve from the inner to the outer heliosphere, leading to transient variations of G_r across large spatial scales, in the same way that the merging of interaction regions in the heliosphere may lead to similar effects (Potgieter et al. 1993; Potgieter & Le Roux 1994; Wang & Richardson 2002). It is, however, questionable whether global merged interaction regions had a big impact during the latter two quiet solar cycles (Luo et al. 2019).