Tracking Cosmic-Ray Spectral Variation during 2007–2018 Using Neutron Monitor Time-delay Measurements

Our techniques to examine GCR spectral variations using time-delay observations were mainly developed using data from PSNM at the summit of Doi Inthanon, Thailand's highest mountain, with geographic coordinates 1859 N, 9849 E and altitude 2560 m. The PSNM has taken data since 2007 August, and since 2007 December 9, it has operated with a complete set of 18 BP-28 (¹⁰BF₃) neutron-sensitive proportional counter tubes with the standard NM64 configuration (18NM64). The 18 tubes are arranged in a single row to reduce edge effects. Further details are available from Ruffolo et al. (2016).

We analyze time-delay observations from 2007 December 9 to 2018 April 18. During theat 10 yr period, we made various changes to PSNM firmware, software, and hardware. Based on these, the data were divided into nine normalization periods, DI-1 to DI-9. Details of changes in the data acquisition system and normalization periods can be found in Appendix A.

Electronic firmware to record time-delay histograms has also been deployed at the South Pole NM. This is located at an altitude of 2828 m. It effectively comprises three 1NM64 detectors, each in a separate temperature-controlled outdoor box, using an ³He neutron-sensitive proportional counter tube. Details of changes in firmware and software in the data acquisition system used to collect time-delay histograms from the South Pole NM, as well as the four normalization periods SP-1 to SP-4, are provided in Appendix B.

Although we do not analyze other data in the present work, for completeness, we point out that these specialized electronics to collect time-delay histograms have been deployed on other NMs. They were first deployed with a shipborne 3NM64 for latitude surveys during 2000–2007, as analyzed by Mangeard et al. (2016a). More recently, they have been deployed on some of the counter tubes at the following NM stations: McMurdo, Antarctica (from 2015 until its closure in 2017); Newark, USA (since 2015); and Jang Bogo, Antarctica (since 2015). They were also deployed on two shipborne NM latitude surveys in 2019 and 2019–2020, and we plan to deploy them at the NM at Mawson, Antarctica, in 2020. In future work, these data should be analyzed to provide a more detailed picture of GCR spectral variations.

The leader fraction L is defined statistically as the fraction of neutron counts that did not follow another count in the same tube due to the same primary cosmic ray, i.e., did not follow another count in the same tube with a short time delay that is not consistent with chance coincidence. More energetic primary cosmic rays lead to more energetic secondary particles that produce more neutrons when interacting inside the NM (typically in the lead producer), which can then be detected as sequences of more counts in the same tube with shorter time delays. Thus, a harder cosmic-ray spectrum with a lower spectral index leads to a lower leader fraction.

In this work, we used the exponential tail fit parameters from the time-delay histogram, which Ruffolo et al. (2016) described as their first method to statistically measure chance coincidences and determine L. These values differ mainly due to each tube's position and electronic dead time (for details, see Mangeard et al. 2016b; Ruffolo et al. 2016). For the PSNM, tubes T1 and T18 at the edges of the row of 18 counter tubes have a lower value of L than the middle tubes, because the middle tubes have a larger fraction of neutron counts from interactions of atmospheric secondary particles in lead rings around neighboring tubes, leading to a lower neutron multiplicity at the tube of interest and a higher leader fraction. The variation of L from all counter tubes clearly shows a yearly wave, which can be attributed to the yearly variation in water vapor and atmospheric pressure (see Appendix A).

The 600 series firmware recorded time delays for only 16 pulses s^–1, while the count rate in the 18NM64 at Doi Inthanon is ≈34 Hz per tube and that at the South Pole is ≈100 Hz per tube, resulting in poorer statistics and a "first-pulse bias" in the time-delay distribution. All L values from 600 series firmware were corrected for the first-pulse bias in the manner described by Ruffolo et al. (2016). The 700 and 800 series firmware record almost all time delays, so the first-pulse bias was negligible for such data.

Another issue before we determine L from the time-delay histograms is the overflow in the electronic recording system. For the 600 and 700 series firmware, the time delay was recorded modulo t_o = 142 ms, and we derived L from the following formula (Ruffolo et al. 2016):

$$\begin{eqnarray}&&L=\displaystyle \frac{1-{e}^{-\alpha {t}_{o}}}{\alpha {e}^{\alpha {t}_{d}}}A,\end{eqnarray} \tag{ 1 }$$

where α and A are the parameters from fitting ${{Ae}}^{-\alpha t}$ to the exponential tail of the normalized measured histogram of time delay t, and t_d is the effective electronic dead time. When the time delays were calculated by the data acquisition software for the 800 series, a count was excluded from the histogram when the time delay was larger than t_o. Hence, for the 800 series firmware, we can use Equation (1) without the overflow correction:

$$\begin{eqnarray}&&L=\displaystyle \frac{A}{\alpha {e}^{\alpha {t}_{d}}},\end{eqnarray} \tag{ 2 }$$

normalizing the histogram to account for missing values at t > t_o.

It is well known that any NM count rate is affected by variations in the atmospheric pressure, P, which serves as a proxy for the overlying airmass. This is due to energy-dependent absorption of secondary neutrons by the atmosphere, so it is unsurprising that the leader fraction L has been found to vary with P as well. Furthermore, previous work inferred that atmospheric water vapor pressure strongly affects the count rates of bare neutron counters at multiple locations and L at Doi Inthanon (Rosolem et al. 2013; Ruffolo et al. 2016). Water vapor even has a small but noticeable effect on the PSNM count rate C at Doi Inthanon, with a peak-to-peak seasonal variation of 0.3% (Mangeard et al. 2018). In this section, we discuss our techniques to correct both C and L as measured at Doi Inthanon and the South Pole for variations in P and (for Doi Inthanon) water vapor pressure.

For the PSNM at Doi Inthanon, it is possible to distinguish the effects of P and water vapor pressure. While both exhibit a yearly (seasonal) variation, there are also strong 12 and 24 hr variations in P due to atmospheric tides. In contrast, water vapor pressure is usually relatively constant over small timescales, and when it does vary rapidly, C and L do not exhibit corresponding variations, as will be discussed shortly. Because the South Pole is an extremely dry location (precipitable water vapor (PWV) typically below 1 mm), the effect of atmospheric water vapor is negligible, and we only correct for variations in P.

First, we focus on the effect of pressure variations, which result in anticorrelated variations in C and correlated variations in L. For C of either the PSNM or the South Pole NM, we use the standard barometric correction developed previously, by which C is multiplied by $\exp \left[\beta (P-{P}_{0})\right]$, where β is a pressure coefficient (0.641% hPa⁻¹ for PSNM and 0.740% hPa⁻¹ for the South Pole) and P₀ is a reference pressure (750.5 hPa for PSNM and 680.0 hPa for the South Pole).

The leader fraction L is quite sensitive to the energy of secondary nucleons (see Figure 10 of Mangeard et al. 2016a), with a much lower L value (higher multiplicity) for GeV-range nucleons that typically originated in the upper atmosphere after at most a few GCR-induced interactions. With increasing atmospheric pressure, GeV-range nucleons are less likely to survive, and L increases. To observationally determine the effect on L, we fit data of ${\rm{\Delta }}\mathrm{ln}{L}_{i}$ versus ΔP to straight lines of slope b_i for each tube i. (Note: Given that the pressure has a positive correlation with L, in contrast with the negative correlation with C, to avoid confusion, we use the quantity b = −β for the pressure correction of L.) A quantity labeled by Δ refers to the residual with respect to a running average, i.e., a short-term fluctuation. For the PSNM, we use an hourly value minus a 24 hr average, so as to examine the effect of short-term P fluctuations, effectively filtering out the effects of water vapor and more gradual variations in the GCR flux. At the South Pole NM, the pressure mainly varies over ∼1 month (and the water vapor effect is negligible), so we use an hourly value minus a 30 day average.

For a given type of firmware, the b values for different counter tubes had minor variations with no clear dependence on the tube position or electronic dead time, so we use the mean value b for all counter tubes. For the PSNM, b for the 800 series firmware (0.0193% hPa⁻¹) was higher than that for the 700 series (0.0154% hPa⁻¹). Because of the similarities between the 600 and 700 series firmware, while L is determined with better statistical accuracy for the latter, we used b from the 700 series for the 600 series as well. For the South Pole NM, b for the 800 series was 0.0191% hPa⁻¹, and that for the 600 series was 0.00922% hPa⁻¹. Once we determined a value of b for the relevant firmware, then L was corrected for pressure variation by multiplying by $\exp \left[-b(P-{P}_{0})\right]$.

Next, let us examine the relationship between atmospheric water vapor and pressure-corrected C and L from Doi Inthanon. Because the effect of water vapor on C is much smaller and, over subannual timescales, is difficult to distinguish from the effect of a changing GCR flux, we analyze the effect on L and expect that the effect on C should be qualitatively similar. Note that previous work found a nonlinear relationship between the surface water vapor pressure E_w and pressure-corrected L, where E_w was derived using data of relative humidity and temperature from the Global Atmospheric Data Assimilation (GDAS) database,⁷ interpolating values from the four surrounding degree-scale spatial grid points to the geographic location of the PSNM, for a pressure level close to the reference pressure at the summit of Doi Inthanon (Ruffolo et al. 2016). The nonlinear relationship is consistent with a saturation of surface water vapor during the rainy season (roughly from April through October), while the effect of water vapor on L continues even as E_w saturates. Note also that Mangeard et al. (2018; see their Figure 2) modeled the effect of E_w on the count rate C using a linear relationship; however, the effect on C is quite small, as noted above, so a possible nonlinearity would be difficult to measure.

Physically, it seems more reasonable that L is affected by water vapor throughout the atmosphere (hence the continued effect on L as the surface water vapor pressure E_w saturates), and that the effect of the total atmospheric water vapor should be linear. Therefore, we used GDAS data to determine the PWV at Doi Inthanon, i.e., the height of liquid water (in mm) that would result if all of the water vapor in the atmosphere were to condense at the surface, to determine whether this would provide a more reasonable model for the water vapor effect on L. The PWV was determined by numerically integrating the water vapor pressure for the atmosphere above Doi Inthanon as inferred from GDAS data, which have a 6 hr cadence. The PWV for the entire period of PSNM data can be found in Appendix A.

Figure 2 shows scatter plots of the pressure-corrected L averaged over all 18 counter tubes with 24 hr smoothing (to reduce the statistical uncertainty) versus (a) E_w and (b) PWV, with the 6 hr cadence of the GDAS data, for the normalization period DI-2 (with 700 series firmware) from 2011 January 15 to 2014 February 8. This time period was chosen because of the good statistical quality provided by 700 series firmware, the constant configuration for 3 yr (see Appendix A), and the nearly constant level of solar modulation for the PSNM at this time (Mangeard et al. 2018), so that time variations in L can mostly be attributed to the water vapor effect, except for occasional Forbush decrease effects that are unrelated to the water vapor pressure. It can be seen that the relationship with E_w is nonlinear, in agreement with Ruffolo et al. (2016), while the relationship with PWV is apparently linear. Presumably, this is due to moderation and absorption of low-energy atmospheric secondary neutrons by water vapor throughout the atmosphere. Thus, a higher PWV can lead to a lower fraction of low-energy secondary neutrons, a higher average multiplicity, and a lower leader fraction. The surface water vapor E_w can serve as a proxy for PWV but often saturates during the rainy season (with nearly 100% relative humidity) so that L has a nonlinear relationship with E_w. In that case, PWV should provide a more direct measure of water vapor content for correcting L and C from the PSNM.

Zoom In Zoom Out Reset image size

Figure 3 shows PSNM data of C, L, and PWV as a function of time for the same time period (DI-2). For C and L, the gray line indicates daily tube-averaged values that are corrected for pressure but not yet for water vapor, while in Figure 3(c), the blue line indicates PWV determined from GDAS data with a 6 hr cadence. Note that there are some relatively sharp decreases in L in association with changes in C, reflecting the effect of GCR spectral and flux variations, respectively. For example, there are sharp Forbush decreases in C (Forbush 1937) associated with the passage of a shock or possibly the ejecta from a solar storm, such as those on 2012 March 8 and 12, and some Forbush decreases are also associated with decreases in L, i.e., GCR spectral variations above the ∼17 GV cutoff (Ruffolo et al. 2016).

Zoom In Zoom Out Reset image size

The remaining nonstatistical variations in L are closely but not perfectly associated with variations in PWV. In addition to the large seasonal wave, L also varies according to shorter-term variations in PWV. Some clear changes in L of this type, while C was nearly constant, occurred from late 2011 December to early 2012 January and during 2013 October. These were associated with changes in PWV that lasted ≳10 days and provide convincing evidence of an association of L with atmospheric water vapor as opposed to another seasonal variation. However, L does not exactly follow more rapid changes in PWV as determined from GDAS. For example, during the rainy season in 2013 June, there was a sharp decline of PWV to a dry-season level, followed by a rapid recovery back to a nearly maximal level. At that time, there was only a small change in L, nowhere near the dry-season level. If we were to correct L based directly on this PWV time series, that would induce spurious spikes in L from overcorrection of shorter-term PWV variations that do not have a proportional effect on L.

We have solved this problem by smoothing the GDAS-derived PWV time series for use in correcting the leader fraction. After testing smoothing of various types and durations, we found a good correspondence between variations of L_Pcorr (as corrected for pressure but not water vapor variations) and those of PWV smoothed with a triangular function that peaks at the actual time and goes to zero over ±5 days. This 10 day triangle-smoothed GDAS-derived PWV time series is shown as the black curve in Figure 3(c) and can be visually seen to be better anticorrelated with the vapor-uncorrected leader fraction (gray) in Figure 3(b) during times of nearly constant C.

To confirm that GDAS provides a useful description of PWV, we compared PWV as determined from GDAS data with PWV determined in other ways. For example, we determined PWV from another database, the Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2),⁸ during 2013. This comparison was quite interesting, as the two values of PWV agreed very well at most times, with a correlation coefficient r² = 0.933, but occasionally sharp features were different or occurred with a time difference of up to 4 days in the two databases. This implies that there can be a temporal uncertainty in the use of such databases to infer the time profile of PWV. Another comparison was with direct measurements from a radiosonde on balloons launched daily at the airport in Chiang Mai, Thailand, which is located 55 km from Doi Inthanon at an elevation of about 300 m above sea level. We compared PWV inferred from the radiosonde and GDAS data for that location during 2013 January–October and found a very good correspondence between the two (r² = 0.926). Finally, we operated a two-frequency GPS receiver at the PSNM site from 2017 November to 2018 June, determining the zenith time delay and subtracting the delay associated with the dry atmosphere to infer the PWV (Takiguchi et al. 2000). We used the techniques of Bevis et al. (1992), with the refinement of using the mean atmospheric temperature T_m as determined from the GDAS database instead of estimating this from the surface temperature. The PWV values obtained from the GPS data had a different absolute normalization, nearly 40% lower than the PWV from GDAS, but the values from GPS and GDAS data were very well correlated (r² = 0.953).

These comparisons give us confidence that GDAS provides a useful measure of PWV. However, the comparison with PWV from MERRA-2 indicates that PWV features can appear at different times when using different databases. This could be partly attributed to our bilinear interpolation of data from the four nearest 1° grid points, which could make features appear at times when, in reality, the feature was not occurring at Doi Inthanon. This could be a reason why temporal smoothing of the inferred PWV is required to better match the actual effects on L. Another possibility is that PWV is in part a proxy for rainwater on or in the concrete and rock surfaces surrounding the PSNM building. (Note that the building itself has a roof slanted at 45° that does not accumulate rainwater; see Figure 1(a) of Ruffolo et al. 2016.) However, we tested pouring water onto these surfaces on two occasions for up to several hours and were unable to detect a change in PSNM data. Furthermore, the effect of rainwater should lead to a time delay, whereas when we tried different smoothing profiles, a delayed profile provided no better results than the symmetric profile that we adopted. For these reasons, we view the need for temporal smoothing as best attributed to temporal uncertainties in our determination of PWV.

To correct L for the effect of water vapor, we performed linear fits to L_Pcorr after 24 hr smoothing versus the 10 day triangle-smoothed PWV calculated from GDAS data with a 6 hr cadence separately for each tube and series of firmware. Then the leader fraction corrected for pressure and water vapor is determined as

$$\begin{eqnarray}&&{L}_{\mathrm{PWcorr}}(t)={L}_{\mathrm{Pcorr}}(t)\displaystyle \frac{1-a\,{\mathrm{PWV}}_{0}}{[1-a\,\mathrm{PWV}(t)]},\end{eqnarray} \tag{ 3 }$$

where a is the negative slope divided by the intercept of the linear relationship. The reference value PWV₀ = 7.2 mm is the median value of PWV in January, during the dry season. The PWV at the PSNM varies from as low as 1 mm in the dry season to a maximum of ≈55 mm in the rainy season.

For a given firmware series, the value of a was quite different between the two end counter tubes and the 16 middle counter tubes, so we used average values for either the end or middle tubes. The values of a for the 700 series firmware were 2.635 × 10⁻⁵ mm⁻¹ for the middle tubes and 3.685 × 10⁻⁵ mm⁻¹ for the end tubes. These are higher than the values for 800 series firmware, which are 2.412 × 10⁻⁵ mm⁻¹ for the middle tubes and 3.221 × 10⁻⁵ mm⁻¹ for the end tubes. As with the pressure correction, for the 600 series firmware, we used the same values of a as for the 700 series firmware.

The leader fraction after vapor correction (black) in Figure 3(b) does not appear to contain spurious short-term spikes due to over- or undercorrection. There are still some systematic variations over scales of months that are not associated with cosmic-ray flux variations as indicated by the count rate C, such as during 2011 December and 2013 January, when C was nearly constant at ≈2.2 × 10⁶ hr⁻¹. Over such times, the pressure- and vapor-corrected L varied from about 0.7816 to 0.7822. We conclude that our vapor correction involves a residual systematic error in L of up to 3 × 10⁻⁴. By averaging L over 3 month intervals, we can reduce that systematic error, as well as statistical errors.

Mangeard et al. (2018) reported a peak-to-peak annual variation of 0.3% in the pressure-corrected count rate, which is attributed to a water vapor effect, compared with the total peak-to-peak solar modulation of ≈3%, so the effect of water vapor is nonnegligible. Our procedure to correct this effect on C is slightly different than for L, because C mostly varies due to true variations in the cosmic-ray flux, e.g., due to Forbush decreases and solar modulation. Therefore, we follow a procedure similar to that of Mangeard et al. (2018), using PWV in place of the surface water vapor pressure E_w. For various ranges of the 10 day triangle-smoothed PWV (bins of 5 mm each), we calculated the median value over 2007 December to 2018 April of the pressure-corrected count rate summed over all 18 counter tubes, C_Pcorr. The median is used in order to reduce the effects of outlier values due to short-term cosmic-ray flux variations. Then we performed a linear fit to the median values of C_Pcorr versus the smoothed PWV. Finally, the count rate corrected for pressure and water vapor is evaluated as

$$\begin{eqnarray}&&{C}_{\mathrm{PWcorr}}(t)={C}_{\mathrm{Pcorr}}(t)\displaystyle \frac{1-a\,{\mathrm{PWV}}_{0}}{[1-a\,\mathrm{PWV}(t)]},\end{eqnarray} \tag{ 4 }$$

where a is found from the negative slope divided by the intercept of the linear fit to be (10.18 ± 2.41) × 10⁻⁵ mm⁻¹. For time period DI-2, Figure 3(a) shows C_Pcorr (gray) and C_PWcorr (black). The latter, corrected for pressure and water vapor variations, is now seen to have the seasonal variation effectively removed, revealing a nearly constant GCR flux punctuated by sporadic Forbush decreases due to solar storms.

After correction for atmospheric effects, it was necessary to normalize L to account for differences among counter tubes, firmware versions, some software versions, and hardware. Details can be found in Appendix C for the PSNM and Appendix D for the South Pole NM. The resulting time series of C and L at Doi Inthanon are shown in Figure 4. Those at the South Pole are shown in Figure 5, where the leader fraction data from 2013 December to 2014 November are plotted in blue to highlight their uncertain normalization relative to later data (in black; see Appendix D for details).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, let us consider the variation of the count rate C and leader fraction L over the recent sunspot cycle from NMs at Doi Inthanon and (for a half-cycle) the South Pole. In each case, C measures the sum of the GCR differential response function (spectrum times NM yield function) for each species, integrated over rigidity above the cutoff (Simpson et al. 1953), i.e., for rigidity ≳17 GV at Doi Inthanon and ≳1 GV at the South Pole. Thus, variations in an NM count rate indicate variations in the GCR flux above the local cutoff. Similarly, for each station, variations in L indicate variations in the GCR spectral index over the rigidity range of the NM response, as shown by Mangeard et al. (2016a) for observations and Monte Carlo simulations of L during a series of NM latitude surveys.

Figure 6 shows 3 month averaged data for the leader fraction L versus count rate C over the whole time period of analysis from Doi Inthanon. The data are averaged over 3 months to reduce the systematic uncertainty from the water vapor correction and the statistical uncertainty. The number at each point indicates the order of the 3 month period, starting with 1 for 2007 December–2008 February and ending with 42 for 2018 March–April (truncated to 2018 April 19 when our analysis ends). One way to estimate the uncertainty of L after the 3 month averaging is to consider the coherence time of the fluctuations. The statistical fluctuations seen in Figure 4 are independent from day to day, so after averaging over 3 months, they are reduced by a factor of about $\sqrt{90}$, making the statistical uncertainty much smaller than the systematic uncertainty. Section 3.2 reported a systematic effect due to imperfect correction for the effect of atmospheric water vapor, with a maximum deviation of 0.0003. This maximum is rarely reached, and we estimate the standard deviation to be about 0.0002. Based on Figure 3, focusing on times when C was nearly constant, so changes in L due to GCR spectral variations should be small, the coherence time of these fluctuations is estimated as about 0.5 months. Therefore, a 3 month average contains about six independent samples of such fluctuations, so the standard error is a factor of $\sim \sqrt{6}$ lower than the standard deviation of the fluctuations (0.0002). We conservatively estimate the uncertainty in the 3 month averages of L from the PSNM as 0.0001. Another way to estimate this is to simply examine the scatter among L from neighboring intervals in Figure 6, which confirms an uncertainty of about 0.0001.

Zoom In Zoom Out Reset image size

Figure 7 shows monthly averaged data for L versus C from the South Pole. Although the total count rate is ≈1.1 × 10⁶ hr⁻¹ compared with ≈2.2 × 10⁶ hr⁻¹ for Doi Inthanon, the South Pole data are not affected by water vapor with an associated systematic uncertainty, and the total variation of L over the solar cycle is much greater, i.e., about 6 × 10⁻³ for the South Pole compared with 8 × 10⁻⁴ for Doi Inthanon. For these reasons, the monthly L data from the South Pole provide useful data for studying solar modulation, currently for half a solar cycle. The numbers in Figure 7 start at 1 for 2013 December and end at 53 for 2018 April.

Zoom In Zoom Out Reset image size

A basic concept of solar modulation is that with increasing solar activity, the cosmic-ray spectrum "rolls" downward from the local interstellar spectrum, and with decreasing solar activity, the spectrum "unrolls" upward again. The change in flux is much greater at low rigidity, while at sufficiently high rigidity, there is almost no solar modulation (Rao 1972). Thus, we expect that as C decreases, L will also decrease. If the GCR spectral shape were the same during the rolling and unrolling, e.g., as in the force-field model of solar modulation, then L versus C would follow the same path during downward and upward trajectories.

For the South Pole NM, the normalization of the initial data for L during SP-1 (2013 December to 2014 November) relative to later data is uncertain (as explained in D), as indicated in blue in Figure 7. After a 3 month data gap (corresponding to the missing numbers 13–15), the data from 2015 March (point 16) to 2018 April (point 53) follow a consistent, nearly linear trend. This is consistent with the "unrolling" of the GCR spectrum in the 1–10 GV range that dominated the variation in South Pole data with solar modulation as solar activity weakened.

For the PSNM at Doi Inthanon, recording L and C over a full solar cycle, we observe a hysteresis loop, implying a change in the GCR spectral shape above the cutoff rigidity of ∼17 GV. In other words, the GCR spectrum had a different shape in this high rigidity range while rolling down and unrolling back up again. Note that the time period of declining cosmic-ray flux was associated with solar modulation from the A < 0 phase of the solar cycle (i.e., with a mostly inward magnetic field at the northern hemisphere of the Sun) as its effects propagated through the heliosphere, and the rising cosmic-ray rate was associated with the A > 0 phase. Thus, our results imply a change in spectral shape with a change in solar magnetic polarity at a rigidity range above the ∼17 GV cutoff.

Such changes have been reported from comparisons between NM measurements at different locations with different cutoff rigidity (Moraal et al. 1989; Nuntiyakul et al. 2014; Mangeard et al. 2018), a technique that is limited to the Earth's range of atmospheric and geomagnetic cutoff rigidities, ∼1–17 GV. In this work, we find distinct changes in the evolution of the cosmic-ray spectral shape at P ≳ 17 GV during periods of opposite solar magnetic polarity, as detailed below.

During 2007 December to 2015 February (points 1–29 in Figure 6), C and L at Doi Inthanon first increased and then declined along a roughly linear trend, which is again consistent with the rolling of the spectrum with greater variation at lower rigidity, so that the spectral index decreased (hardened) along with the count rate decrease. This time period was mostly associated with A < 0 solar magnetic polarity (Mangeard et al. 2018). However, shortly after the cosmic-ray minimum in early 2015, now with A > 0 solar magnetic polarity, L shot upward to a nearly maximal value while C increased only slightly, a qualitatively different pattern from that observed at the South Pole. This implies little change in the count rate while the spectral index above the cutoff rapidly increased, indicating a rapid recovery of the GCR flux at about 17 GV, while the flux at higher energy had little or no increase or perhaps even declined to result in the nearly constant count rate. This was a significant reorganization of the GCR spectrum above ∼17 GV within a short time period in early 2015.

From 2015 March to 2018 April, the count rate at Doi Inthanon mostly had a gradual increase (punctuated by sharp 27 day variations), while the L and spectral index remained nearly constant, in contrast with the behavior during 2007–2014 in the A < 0 phase of modulation. The nearly constant spectral index after early 2015 implies a nearly rigidity-independent increase in the flux over a significant rigidity range above the cutoff.

We have noted that the leader fraction L serves as an indicator of the GCR spectral index, so it is interesting to estimate a relationship between the two. This leads to the question: what is the rigidity (or energy) for which L indicates the spectral index? For the PSNM, the median primary GCR rigidity is 35 GV, and the 10th–90th percentile range of its response is from 18 to 155 GV (Ruffolo et al. 2016). However, solar modulation has a very low effect on the higher rigidities in that range. Based on measurements by the AMS-02 detector on the International Space Station over the time period of 2011 May to 2013 November (Aguilar et al. 2015a, 2015b), we see that the primary cosmic-ray spectral index was roughly independent of rigidity above 40 GV (within the uncertainties), which could indicate that solar modulation has little effect at P > 40 GV and hardens the spectrum at lower rigidity. This view is supported by AMS-02 measurements indicating little or no time variation in the proton flux at 41.9–45.1 GV from 2011 May to 2017 April (Aguilar et al. 2018). Therefore, we estimate that the variation of L at Doi Inthanon due to solar modulation mainly reflects the variation of the GCR spectral index over ∼20–40 GV. For a polar NM, the response includes a very wide range of rigidity over which solar modulation has strong effects, and the spectral index strongly varies with rigidity, so we do not attempt to relate L at the South Pole NM to a spectral index over a specific rigidity interval.

Even for the PSNM, estimating the relationship between L and a spectral index over ∼20–40 GV brings up further questions. What spectral form should we assume? Can the variation of the GCR spectrum above ∼20 GV due to solar modulation be characterized by varying a single parameter? The second question is clearly answered by our discussion in the preceding subsection: solar modulation of the GCR spectrum above the ∼17 GV cutoff at Doi Inthanon is not well characterized by a single parameter. This is clear from the hysteresis in Figure 6; different times with the same count rate can have different leader fractions, indicating a different spectral shape. Nevertheless, in order to provide a rough calculation of how changes in the spectral index over ∼20–40 GV correspond to changes in L, we have performed a calculation based on broken power-law spectra for protons and for ions with Z ≥ 2. The break in the power law is set at 40 GV. The power-law indices for >40 GV are fixed at values based on results from AMS-02 (Aguilar et al. 2015a, 2015b), while the power-law index γ_<40GV for protons at <40 GV is allowed to vary (i.e., due to solar modulation) and the power-law index for Z ≥ 2 varies in tandem. We stress that we are not proposing this broken power law as a proper description of the GCR spectrum.

Appendix E explains how we obtained yield functions from Monte Carlo simulations and used the model spectrum described above to calculate ${dL}/d{\gamma }_{\lt 40\mathrm{GV}}=0.0042$. Figure 6 indicates that the full variation of L over the past solar cycle was about 0.0008. This implies a variation in ${\gamma }_{\lt 40\mathrm{GV}}$ of ≈0.2 over the solar cycle, as measured over the rigidity range from the ∼17 GV cutoff at Doi Inthanon to ∼40 GV.

Previous work reported short-term spectral hardening as indicated by decreases in L from Doi Inthanon (Ruffolo et al. 2016) during Forbush decreases in C. These short-term decreases are associated with solar storms (Forbush 1937), in particular, the passage of a coronal mass ejection and/or an associated shock (Cane 2000). The decrease is stronger at lower rigidity, hence the spectral hardening. Selected Forbush decreases are indicated in Figures 4 and 5. These are times when C and L decrease very rapidly (typically for ≲1 day) and recover (increase) to near their prior levels for up to several days. Figure 4 confirms decreases in L (spectral hardening) at Doi Inthanon to various degrees for the stronger Forbush decreases in the past solar cycle.

In Figure 5, we see decreases in L much more clearly in data from the South Pole NM at times of Forbush decreases in C. As noted earlier, variations in L are clearer at the South Pole because GCR variations are stronger at low cutoff rigidity, and also because the South Pole is a very dry location where there is no water vapor effect or associated systematic uncertainty.

Figure 5 also exhibits a sharp increase in L on 2017 September 10. This is associated with the ground-level enhancement on that day (Kurt et al. 2018), as indicated in the figure, due to a solar storm that produced a sufficiently high flux of relativistic solar particles that there was a measurable flux enhancement over the GCR flux, especially at low cutoff rigidity, due to the much softer spectrum of solar particles. Interestingly, the ground-level enhancement is more clearly visible in the leader fraction L than in the count rate C, indicating that spectral signatures can assist in the detection of ground-level enhancements that involve only a weak count rate enhancement.

A noteworthy feature of C versus time at both Doi Inthanon and the South Pole (in Figures 4 and 5) is the 27 day variation (Fonger 1953), which was particularly prominent from late 2014 to 2018. This varies with the 27 day period of solar rotation, because the solar wind speed depends on the source location at the Sun. In particular, there can be high-speed solar wind streams in interplanetary space that collide with the preceding slower solar wind to form corotating interaction regions (e.g., Hundhausen et al. 1980), and these corotating features can be associated with recurrent decreases in the GCR flux. In the data for L at the South Pole, for the first time, we can confirm 27 day variations in L together with those in C. The implied spectral hardening during recurrent decreases in C is again expected because the 27 day features are fractionally stronger at lower rigidity (Gil & Alania 2010).

Note that in Figures 4 and 5, C and L are plotted to show the variations with solar modulation on a similar scale. With this scale, the shorter-term variations in L are much smaller than the variations in C. In other words, variations in L relative to the 27 day variations and Forbush decreases in C are much weaker than those for solar modulation. This indicates that the physical processes causing such variations do not roll and unroll the cosmic-ray spectrum as strongly as solar modulation. In other words, 27 day variations and Forbush decreases are not as strongly enhanced at low rigidity as solar modulation.