1 Introduction

Galactic cosmic rays (GCRs) are charged particles and atomic nuclei with energies spanning the range from a few MeV up to approximately \(10^{21}~\mbox{eV}\), which encroach upon the Earth from all directions (Giacalone, 2010). They mainly originate outside the solar system, within the Milky Way; however, they are also expected to originate from other galaxies (Aab et al., 2017). GCRs at the top of the atmosphere are mostly composed of protons (\(\approx 87\)%) and \(\alpha \)-particles (\(\approx 12\)%), with a smaller contribution (\(\approx 1\)%) from heavier nuclei (Dunai, 2010).

When cosmic rays (CRs) enter the atmosphere, they interact with atmospheric atoms and produce cascades of secondary particles, which at ground level are primarily neutrons and muons. Neutron monitors (NMs) and muon detectors (MDs) located at different locations on Earth have been used since the 1950s to observe GCRs. Information on GCRs prior to the modern epoch of NMs and MDs, and the space age, rely on the studies of cosmogenic isotope records from ice cores and tree rings (Owens and Forsyth, 2013).

It has long been established that there exists an anti-correlation between GCR intensity and the level of solar activity over a cyclic 11-year period, with perhaps some time lag (Forbush, 1958; Parker, 1965; Usoskin et al., 1998; Van Allen, 2000). Figure 1 clearly shows the anti-correlation between GCRs and sunspot number (SSN).

Figure 1
figure 1

SSN (top), with vertical lines showing the beginning of each solar cycle. CR intensity recorded by NMs (bottom), with vertical lines showing the approximate epochs of solar magnetic-field polarity reversals. (MCMD = McMurdo, NEWK = Newark, SOPO = South Pole, THUL = Thule).

Full size image

It is well known that the 11-year solar-activity cycle is in fact a 22-year cycle – the Hale cycle – that describes the alternating polarity of the large-scale solar magnetic field (Thomas, Owens, and Lockwood, 2014). The alternating peaked and flat-topped shape of GCR intensity in Figure 1 is a manifestation of this effect in addition to other CR transport processes (Aslam and Badruddin, 2012).

The polarity of the solar magnetic field \([A]\) is taken to be negative when the field axis is aligned with the axis of rotation, and positive when the opposite is true (Thomas, Owens, and Lockwood, 2014). The solar field polarity is conventionally described in combination with the particle charge \([q]\) that is due to the effect of curvature and gradient drift on charged particles; thus it is customary to define the solar polarity as \(qA\). Vertical lines showing the approximate epochs at which the polarity reverses are plotted in Figure 1 (Thomas, Owens, and Lockwood, 2014; Janardhan et al., 2018).

Particle drifts differ during different \(qA\) cycles, with positive CRs (i.e. protons) predominantly arriving into the heliosphere from the heliospheric poles and outwards to Earth during periods when \(qA > 0\), whereas when \(qA< 0\) positive CRs predominantly arrive at Earth inwards along the heliospheric current sheet (HCS) (Belov, 2000; Thomas et al., 2014). As the solar magnetic dipole axis is tilted to the solar rotation axis, so is the HCS; the HCS tilt varies with solar cycle and is typically smaller during solar minimum and larger during solar maximum (Owens and Forsyth, 2013). The tilt angle of the HCS has also been shown to be strongly correlated to the GCR intensity and the GCR lag behind the solar activity (Belov, 2000; Mavromichalaki, Paouris, and Karalidi, 2007).

Aslam and Badruddin (2012, 2015) found that the different processes of CR transport have varying levels of importance throughout the solar-activity cycle, but around solar maximum, it is likely that drifts play less of a role, and disturbances in the solar wind (and hence HCS) are the predominant factor of CR modulation. Even-numbered cycles encounter \(qA<0\) polarity during their onset phase and \(qA>0\) during their declining phase, thus experiencing a faster GCR recovery after solar maximum as the GCRs predominantly enter the heliosphere from the heliospheric poles and experience an outwards drift towards Earth. Odd cycles encounter \(qA>0\) polarity during their onset phase and \(qA<0\) during their declining phase and so experience a slower recovery after solar maximum, as the GCRs predominantly enter the heliosphere along the HCS. When the HCS is tilted and disturbed during the declining activity phase, the path length that GCRs must travel to Earth increases, resulting in an increased time lag.

Several studies have demonstrated that the lag between GCR and solar-activity proxies is approximately zero (i.e. no lag) during even solar cycles, and that there exists a lag of around a year or more during odd solar cycles (Usoskin et al., 1998; Mavromichalaki, Paouris, and Karalidi, 2007; Singh, Singh, and Badruddin, 2008).

Van Allen (2000) showed, through plotting the annual mean intensity of GCRs against sunspot number between 1953 and 1999 (covering Solar Cycles 19 – 22), that there is a distinct difference in the plot shapes between the different solar cycles, with Cycles 19 and 21 producing broad ovals, and Cycles 20 and 22 approximately flat lines. The striking difference between odd- and even-numbered cycles is shown for Cycles 19 – 24 in Figure 2. This hysteresis effect is probably caused by the combination of the heliospheric magnetic field (HMF), solar magnetic-field polarity and thus the particle drift, and the tilt of the HCS leading to a slow recovery of GCR intensity after maxima in odd cycles and a fast recovery after maxima in even cycles (Van Allen, 2000; Belov, 2000; Thomas, Owens, and Lockwood, 2014).

Figure 2
figure 2

Hysteresis plots between yearly averaged SSN and yearly averaged GCR intensity for each of the four main NM stations over Cycles 19 – 24.

Full size image

An extension of this work has since been carried out by Inceoglu et al. (2014), showing that the even-numbered solar-activity cycles can be best modelled using a linear fit due to the narrow shape of the hysteresis loops, whereas odd-numbered solar-activity cycles are better represented by ellipses due to their shape.

There has been speculation in the literature on the behaviour of Cycle 24 compared to recent odd and even cycles. It has been suggested that there exists a lag between SSN maxima and GCR intensity minima in excess of ten months (Kane, 2014; Mishra and Mishra, 2016), which does not follow the pattern of previous even cycles having a near-zero lag and suggestions that Cycle 24 behaved similarly to previous odd cycles; however, these studies do not make use of a complete cycle of data and thus may draw inaccurate conclusions about the behaviour of the whole cycle because of the unusually extended nature of the declining phase of Cycle 23 and the low amplitude of the Cycle 24 maximum (Broomhall, 2017). Mishra and Mishra (2016) made use of a more complete data set for Cycle 24 that is still incomplete, however, and concluded that it is also likely that a four-month lag could exist between GCRs and SSN.

The presented work aims to provide a timely update on the statistical relationship between GCR intensity and solar activity during Solar Cycle 24, since the cycle has now almost declined to a minimum. These aims have been achieved through a time-lag analysis and hysteresis-effect analysis between SSN and GCR intensity.

In Section 2 we provide a brief description of the data that were used throughout this study for both CRs and SSN.

We show in Section 3 through a correlative time-lag analysis that there exists a small time-lag between the SSN and GCR intensity over Solar Cycle 24, which is slightly longer than preceding even-numbered cycles but not as high as observed in previous odd-numbered cycles. We also discuss whether the time lag between SSN and GCR shows a dependence on the rigidity cut-off of the observing station.

In Section 4 we model the shapes of hysteresis plots between GCR intensity and SSN. We show that the behaviour of the hysteresis loops for Cycle 24 follows the preceding even-numbered solar-activity cycles and is better represented by a straight-line fit rather than an elliptical model.

2 Data

For the majority of the work in this study, we have considered the pressure-corrected count rates measured by four NM monitor stations as acquired from the NM data base (NMDB) event search tool (NEST: nmdb.eu/nest/ ). The four stations are McMurdo (MCMD), Newark (NEWK), South Pole (SOPO), and Thule (THUL), i.e. the same NM stations as were used by Inceoglu et al. (2014) to provide a comparison to existing literature. Table 1 details the basic characteristics of the NM stations we used here.

Table 1 Neutron-monitor stations used in this study and their vertical geomagnetic cut-off rigidity [\(R_{c}\)], longitude, latitude, and altitude acquired from NEST. The first 4 stations have been used for all of the analyses, while the lower 12 stations have been used exclusively for the investigation into the dependence of \(R_{c}\) on the time lag.
Full size table

We have investigated the long-term GCR modulation in the heliosphere from 1964 – 2018, spanning Solar Cycles 20 – 24, for the cycle epochs of Cycle 20 (October 1964 – March 1976), Cycle 21 (March 1976 – September 1986); Cycle 22 (September 1986 – May 1996), Cycle 23 (May 1996 – December 2008), and Cycle 24 (December 2008 – March 2018). Early predictions on Solar Cycle 25 suggest that Solar Cycle 24 is unlikely to reach a minimum earlier than the middle of 2019 up to as far as early 2021 (Howe et al., 2018; Upton and Hathaway, 2018; Pesnell and Schatten, 2018). The data used in this study are therefore of an incomplete Cycle 24; however, we believe this to have a minimal effect on the results as Cycle 24 draws to a minimum. Cycle 19 was omitted from this study due to the incomplete dataset for this period (see Figure 1 and Figure 2).

During the time-lag correlation analysis, as our results suggested there may be a rigidity dependence on the time lag, we introduced a further 12 NM stations with data acquired from NEST spanning Cycles 20 – 24 to increase the rigidity spectrum that we used. These stations and their basic characteristics are also detailed in Table 1. These stations are not included in the rest of the results, however, as the results from these stations do not change the conclusions of this study.

Furthermore, we have also used monthly and yearly averaged SSN, as collected by WDC-SILSO ( sidc.be/silso/ ), for the time-lag analysis and hysteresis analysis, respectively, as our chosen proxy of solar activity.

3 Time-Lag Analysis

3.1 Method

To investigate the time delay between the modulation of GCRs compared to the solar activity, a time-lag cross-correlation analysis was performed between monthly mean GCR intensity and monthly mean SSN for each station, following the approach of Usoskin et al. (1998). We used a temporal window of width \(T\) centred on a time \(t\), i.e. shifting within the interval \(t-T/2\) to \(t+T/2\). Here we used \(T = 50\) months.

The window was shifted in steps \(\Delta t = 1\) month within this interval, and for each step the Spearman’s rank correlation coefficient [\(\rho \)] between GCR intensity and SSN was calculated. The lag between GCR and SSN was then estimated by finding the peak correlation coefficient within the time interval \(T\).

The results from the 4 main NM stations used suggested that there may be a relationship between the rigidity cut-off [\(R_{c}\)] of an NM station and the time lag for GCRs; hence 12 additional NM stations were introduced to determine whether this was so, as detailed above.

3.2 Results

The correlation [\(\rho \)] between monthly averaged GCR counts and SSN for different time lags was calculated for Cycles 20 – 23. The variation in \(\rho \) is presented in Figure 3, showing that for each cycle there is a time-lag corresponding to peak anti-correlation between GCR intensity and SSN. Table 2 summarises the time lag with the highest correlation and the corresponding correlation coefficient for all stations in each individual solar cycle.

Figure 3
figure 3

Variation in the correlation coefficient with time-lag NM station GCR intensity and SSN during Solar Cycles 20 – 23.

Full size image
Table 2 Time lags and the corresponding cross-correlation coefficient between NM CR count and SSN for Solar Cycles 20 – 23.
Full size table

As previously reported in the literature, we see here that all of the NM stations clearly exhibit almost no lag during even solar cycles, and a longer lag varying between 11 – 21 months during odd solar cycles. There is a strong agreement between the results presented in Table 2 and those of Mavromichalaki, Paouris, and Karalidi (2007), Kane (2014), and Paouris et al. (2015), thus providing further evidence on the distinction between odd and even solar cycles due to particle transport in the heliosphere. The agreement with existing literature provides evidence of a suitable methodology in this study.

The same cross-correlation technique was then applied to Cycle 24 between December 2008 to March 2018, and the results are presented in Figure 4 and Table 3.

Figure 4
figure 4

Variation in the correlation coefficient with time lag between NM GCR intensity and SSN during Solar Cycle 24.

Full size image
Table 3 Time lags and the corresponding cross-correlation coefficient between NM GCR intensity and SSN for Solar Cycle 24.
Full size table

Cycle 24 is seen here to follow the pattern of almost no lag for even cycles; however, Cycle 24 does display a lag that is larger than that seen in the previous two even-numbered cycles, despite not being as large as that seen in the two previous odd-numbered cycles. The cause for the increased time lag in Cycle 24, as compared to the previous two even-numbered cycles, is likely the combined effects of the unusually deep and extended minimum between Solar Cycles 23 and 24, which delayed the decline in GCR intensity and caused record-breaking high GCR intensities (Pacini and Usoskin, 2015), and the small amplitude of the Cycle 24 maximum.

The results presented in this study, using data for a near-complete Cycle 24, show that the results of Kane (2014) and Mishra and Mishra (2016) were likely unduly influenced by the unusually deep and extended declining phase of Cycle 23 given that they had a limited dataset. Mishra and Mishra (2016) used data for just over half of Cycle 24, and this resulted in a time lag of four months, which agrees with the results of this study.

As a further note on time lag, Tomassetti et al. (2017) showed that through the introduction of time lag as a parameter in the CR transport calculations of CR spectra that there exists a time lag of \(8.1 \pm 1.2\) months during the period 2000 – 2012 spanning across Cycles 23 and 24. We performed the time-lag analysis for the first four NM stations detailed in Table 1 for the period between 2000 – 2012 to investigate whether these results can be reproduced. The results of this analysis are presented in Figure 5 and Table 4.

Figure 5
figure 5

Variation in the correlation coefficient with time lag between NM GCR intensity and SSN between 2000 – 2012.

Full size image
Table 4 Time lags and the corresponding cross-correlation coefficient between NM GCR intensity and SSN during 2000 – 2012.
Full size table

From the time-lag analysis of the four stations that we presented, there is a mean lag of \(10.00 \pm 1.47\) months, which is in good agreement with the results of Tomassetti et al. (2017).

Finally, allowing for the odd- and even-cycle dependence, we see in Figures 3, 4, and 5 that the time lag appears to depend on the rigidity of the NM station used for observation. Such a dependence may impact the conclusions depending on the choice of NM station. We expect that if a dependence exists, a station with a higher rigidity cut-off [\(R_{c}\)] would have a shorter lag as this station observes higher-energy CRs, which are affected less by solar modulation and thus are able to recover faster from solar maximum. Conversely, a station with a lower cut-off rigidity observing lower-energy CRs, which are more influenced by solar modulation, would recover more slowly from solar modulation and therefore experience a longer time-lag. This is supported by Figures 3, 4, and 5, but in order to provide more conclusive evidence of such a relationship, we introduced the additional NM stations detailed in Table 1 to provide a rigidity range spanning 0 – 13 GV. We present in Figure 6 a plot of the time-lag versus station \(R_{c}\) for all 16 stations over Cycles 20 – 24. To acquire uncertainties on the time lag, we ran 1000 Monte Carlo simulations of the time-lag analysis, sampling from the uncertainty distributions for each of the monthly averaged SSN and GCR counts; however, the uncertainties in the data propagated in the Monte Carlo simulations produced no appreciable scatter in the overall results.

Figure 6
figure 6

Variation in time lag plotted against NM station rigidity cut-off for the 16 NM stations detailed in Table 1.

Full size image

The results of this analysis do not show a clear dependence of rigidity on the time lag between SSN and GCR intensity; the sampling of higher \(R_{c}\) is too low because of the availability of high-\(R_{c}\) stations to reasonably conclude on such a dependence at high rigidities, although Cycle 21 suggests the existence of a dependence for low-\(R _{c}\) stations as we expected. For low-\(R_{c}\) stations there appears to be a more pronounced distinction between the time lag observed between odd- and even-numbered cycles than for higher-\(R_{c}\) stations, but again, this is not definitive because of the low sampling at higher-\(R _{c}\). We therefore conclude that there will be no significant dependence of the time-lag analysis on the \(R_{c}\) of the observing station.

4 Hysteresis Effect Analysis

4.1 Method

To investigate the hysteresis effect, we have adopted the approaches of Van Allen (2000), Singh, Singh, and Badruddin (2008), and Inceoglu et al. (2014). Plots of the annual mean SSN versus the annual mean GCR intensity were generated for Cycles 20 – 24 and analysed by fitting different models to the data.

As highlighted by Inceoglu et al. (2014), even-numbered solar cycles can be suitably modelled by a linear fit because of their narrow hysteresis shape, and odd-numbered solar cycles were shown to be better modelled by ellipses because of their broadened hysteresis shape. Here, we first repeated for Solar Cycles 20 – 23 the linear and ellipse fitting to confirm that the method reproduces the results reported by Inceoglu et al. (2014), before we applied the method to Cycle 24.

For even cycles, which display narrow hysteresis loops, an unweighted least-squares linear regression was used to reconstruct estimates of the GCR intensity from SSN. As odd-numbered solar cycles display a broader hysteresis loop, they were separately modelled using unweighted linear regression and ellipse fitting to determine the model that provided the better fit.

The equation of the ellipse fitting took the form

$$ \left [ \textstyle\begin{array}{c} x \\ y \end{array}\displaystyle \right ] = \left [ \textstyle\begin{array}{c} x_{0} \\ y_{0} \end{array}\displaystyle \right ] + R(\phi ) \left [ \textstyle\begin{array}{c} a \, \cos {\theta } \\ b \, \sin {\theta } \end{array}\displaystyle \right ] , $$
(1)

where \(x\) is the GCR intensity, \(y\) is the SSN, \((x_{0}, y_{0})\) are the centroid coordinates of the fitted ellipse, \(R(\phi )\) is the rotation matrix, \(\phi \) is the ellipse tilt angle, \(a\) and \(b\) are the semi-major and semi-minor axes, respectively, and \(0 \leq \theta \leq 2\pi \) is the polar angle measured anti-clockwise from the semi-major axis.

In order to reconstruct the GCR intensity from the model, where linear regression was used to model the data, GCR intensity was derived directly from the SSN for each year. For the ellipse model, the GCR intensity was derived from the model as a function of time using \(\theta \), where the time lag calculated from the analysis in Section 3 was used to correctly phase the ellipse allowing \(\theta \) to be calculated using standard equations of ellipses.

The GCR intensity values predicted by the linear regression and ellipse models were compared to the measured GCR intensity using Spearman’s rank correlation following Inceoglu et al. (2014).

4.2 Results

The hysteresis loops between yearly averaged SSN and GCR intensity for each station were first modelled with a linear regression for both odd and even solar-activity cycles, and then the odd cycles were separately re-modelled by ellipse fitting to show that this provides a more representative fit, as suggested by Inceoglu et al. (2014). The correlation between measured CR intensities and modelled CR intensities for Cycles 20 – 23 is presented in Table 5.

Table 5 Correlation coefficients of the linear regression and ellipse modelling of the hysteresis plots for Solar Cycles 20 – 23.
Full size table

There is a consistent and good agreement between the measured and modelled CR intensities for even solar cycles modelled through linear regression because the hysteresis loops are quite narrow, as shown in Figure 7. These results support the findings of Inceoglu et al. (2014). We note that discrepancies in the correlation coefficients between this study and Inceoglu et al. (2014) are likely due to a number of reasons: Inceoglu et al. (2014) used data-smoothing processes, while we used raw data; Inceoglu et al. (2014) made use of monthly mean data, while we used annual mean data; and Inceoglu et al. (2014) interpolated missing data, whereas we left gaps untreated.

Figure 7
figure 7

Hysteresis plots for even Solar Cycles 20 and 22 and the linear regression fit to the data.

Full size image

The linear relations for odd solar cycles are less consistent in their agreement with observed CR intensities, with the correlation during Cycle 21 as low as 0.34 for South Pole, and as high as 0.91 for Thule. Across both of the odd cycles considered in this study, linear regression is not as good a representation of the data as for even cycles. Figures 8 and 9 show the wider hysteresis loops, which is a characteristic of odd solar cycles and visually shows that a linear fit does not provide a good representation of the data.

Figure 8
figure 8

Hysteresis plots for odd Solar Cycles 21 and 23 and the linear regression fit to the data.

Full size image
Figure 9
figure 9

Hysteresis plots for odd Solar Cycles 21 and 23 and the ellipse fit to the data.

Full size image

In agreement with the results of Inceoglu et al. (2014), it can be seen from the results in Table 5 and Figure 9 that the ellipse models provide estimates of the CR intensity that are in good agreement with the measured intensities as a result of the increased correlation coefficient for each station during Cycles 21 and 23; for South Pole during Cycle 21, the increase in \(\rho \) is seen to be 0.58, which proves the benefit of the ellipse model.

If Cycle 24 follows the pattern of different behaviour between odd and even cycles, it is expected that the best fit will be provided by the linear model; however, Figure 2 shows that Cycle 24 appears to display a wider hysteresis loop than the two preceding even-numbered cycles, as shown in Figure 10. The linear model and the ellipse model were both applied to the hysteresis plots for Solar Cycle 24 to determine which model would provide the better fit; the correlation between measured CR intensities and modelled CR intensities are presented in Table 6.

Figure 10
figure 10

Hysteresis plot for Solar Cycle 24, and the linear regression fit to the data (left) and ellipse fit to the data (right).

Full size image
Table 6 Correlation coefficients of the linear regression and ellipse modelling of the hysteresis plots for Solar Cycle 24.
Full size table

The linear model for Cycle 24 shows a good correlation between the observed and modelled CR intensities, providing evidence to suggest that Cycle 24 follows the two preceding even-numbered cycles. The ellipse model does improve the relation between the observed and modelled CR intensities for two out of the four stations, however: McMurdo and Thule. For South Pole and Newark, the ellipse model was not able to provide a fit at all, which is believed to be due to the Newark data points crossing where the semi-major axis would be defined for the ellipse model, causing the calculation of the semi-major and semi-minor axes to not be defined, and the South Pole has data missing at the beginning of Cycle 24.

The results for Solar Cycle 24 do not provide a conclusive answer as to whether Cycle 24 behaves like past even- or odd-numbered cycles from this dataset alone; however, the ellipse model does not provide as significant an improvement for the two modelled NM stations as it does for odd cycles. The small improvement in the ellipse fit is again likely due to the effects of the extended declining phase of Cycle 23 and the unusually low activity of Cycle 24.

We repeated the analysis for the additional four NM stations that featured in this study. Again, the linear model provided a good fit to the data, but the ellipse model was not able to provide a fit; this favours the conclusion that Cycle 24 is best represented by a simple linear model, as was true for the preceding even-numbered cycles.

Despite Cycle 24 having not yet declined to a minimum, it is clear from the observations shown in the hysteresis plots that additional data in Cycle 24 are unlikely to broaden the loop any further. The hysteresis loop begins to tighten up after 2016 following the broadening between 2014 – 2016; hence it appears unlikely that by the end of Cycle 24, further observations will support the ellipse model.

5 Conclusions

As CRs are modulated by the heliosphere during the 11-year solar activity cycle, and this effect has been studied for previous solar cycles, the principal aim of this study was to investigate the nature of GCRs during the current activity of Cycle 24 as it draws to a minimum.

We presented a time-lag analysis between GCR intensity and SSN, which showed that Cycle 24 has a longer lag (two to four months) than the preceding even-numbered solar activity cycles (typically zero to one month); however, its lag is not as large as those of preceding odd-numbered cycles, and Cycle 24 follows the trend of a short or near-zero lag for even-numbered cycles. We suggest here that the extended lag in Cycle 24 compared to previous even-numbered cycles is likely due to the deep, extended minimum between Cycles 23 and 24, and the low maximum activity of Cycle 24 (Broomhall, 2017).

It has previously been shown in the literature that there is a striking difference in the shape of the relation between SSN and GCR intensity between odd- and even-numbered solar cycles. Because of the difference in the shape of the hysteresis plots for odd- and even-numbered cycles, we have modelled the hysteresis plots using both a simple linear model and an ellipse model. The results of this study tend to support that Cycle 24 follows the same trend as preceding even-numbered cycles and is best represented by a straight line rather than an ellipse; such is the case for odd-numbered activity cycles.

We emphasise that although Cycle 24 has not ended yet, the shape of the hysteresis plots suggests that we are now past the main broadening region, and the inclusion of further data for Cycle 24 will very likely only support the linear model. This study will continue to follow the evolution of Cycle 24 until the onset of Cycle 25, in around 2019 – 2021 (Howe et al., 2018; Upton and Hathaway, 2018; Pesnell and Schatten, 2018), when an update on the final results of Cycle 24 should be provided.