Extracting Hale Cycle Related Components from Cosmic-Ray Data Using Principal Component Analysis

Although cosmic-ray intensity lags somewhat behind the solar sunspot numbers, we have used the dates mentioned in Table 1 for the cosmic-ray data cycles in the first PCA analysis. We used this choice first because the lags from solar cycles vary between different studies (Mavromichalaki et al., 1997; Gupta, Mishra, and Mishra, 2005; Rybansky, Kudela, and Minarovjech, 2009; Kane, 2014; Porta, 2018; Iskra et al., 2019; Ross and Chaplin, 2019; Koldobskiy et al., 2022). We also varied the starting dates in order to see if this affects to the results of the PC analysis. The other dates, which we use in this paper are the local maxima between different cosmic-ray cycles. These are shown in Table 2. We used in our analyses only cosmic-ray data, which have at least four total cycles, excepting with a few values missing due to gaps in the measurement. For this reason we chose Oulu CR flux as a primary data, because it consists five whole Cycles 20 – 24. We also study Hermanus and Climax (Colorado) CR as another primary data, because the former station is located in the southern hemisphere and the latter is the only reliable station with continuous values for the whole Cycles 19 – 22. Note, however, that the early neutron monitors were more simple, called ’IGY’ type monitors, and only since 1964 the improved NM64 monitors have been used (Väisänen, Usoskin, and Mursula, 2021). In comparison for these data we use data from Kerguelen (20 – 24), McMurdo (20 – 23), Moscow (20 – 23), Newark (20 – 24) and Thule (20 – 24). Figure \iref\thearticlefig:Oulu_Climax_ja_muuta and b show the cosmic-ray intensities at Oulu and Climax neutron monitor (NM) stations for the Cycles 20 – 24 and 19 – 22, respectively. We have added the other station to these figures also (see captions of the Fig. \iref\thearticlefig:Oulu_Climax_ja_muut). Because the cutoff rigidities and amount of NMs differ in different stations, i.e., the levels of the counts vary a lot (Usoskin et al., 2017; Väisänen, Usoskin, and Mursula, 2021), we normalize here the CR intensities such that the maximum of each station is one. The normalizing method does not affect the PCA, because we do analyses separately for each station data. Note also from the Fig. \iref\thearticlefig:Oulu_Climax_ja_muut that the maxima of the CR intensities are lagging from the minima of the solar cycles shown as vertical lines. Another normalization, which we have to do, is resampling (interpolation) of the cycles. This is because each cycle of the individual station must have same length to be able to present the cycles as a correlation matrix needed for the PCA. This procedure ensures that each cycle has equal weight to the common cycle shape, although their amplitudes are different. The cycles are, however, reverted to their original lengths and amplitudes after the PCA process.

Principal component analysis is a useful tool in many fields of science including chemometrics (Bro and Smilde, 2014), data compression (Kumar, Rai, and Kumar, 2008) and information extraction (Hannachi, Jolliffe, and Stephenson, 2007). PCA finds combinations of variables, that describe major trends in the data. PCA has earlier been applied, e.g., to studies of the geomagnetic field (Bhattacharyya and Okpala, 2015), geomagnetic activity (Holappa, Mursula, and Asikainen, 2014; Takalo, 2021b), ionosphere (Lin, 2012), the solar background magnetic field (Zharkova et al., 2015), variability of the daily cosmic-ray count rates (Okpala and Okeke, 2014), and atmospheric correction to cosmic-ray detectors (Savić et al., 2019). As far as we know, this is the first time to study and compare cosmic-ray cycles using PCA.

To this end, we estimate that the average length of the cycle is 133 months, and use it as a representative cosmic-ray cycle. We first resample the monthly cosmic-ray counts such that all cycles have the same length of 133 time steps (months), i.e about the average length of the Solar Cycles 20-24 (Takalo and Mursula, 2018, 2020; Takalo, 2021b). This effectively elongates or abridges the cycles to the same length. Before applying the PCA method to the resampled cycles we standardize each individual cycle to have zero mean and unit standard deviation. This guarantees that all cycles will have the same weight in the study of their common shape. Then after applying the PCA method to these resampled and standardized cycles, we revert the cycle lengths and amplitudes to their original values.

As said earlier, we standardize each individual cycle of CR flux to have zero mean and unit standard deviation. Standardized data are then collected into the columns of the matrix $X$, which can be decomposed as (Hannachi, Jolliffe, and Stephenson, 2007; Holappa et al., 2014; Takalo and Mursula, 2018)

where $U$ and $V$ are orthogonal matrices, $V^{T}$ a transpose of matrix $V$, and $D$ a diagonal matrix $D=diag\left(\lambda_{1},\lambda_{2},...,\lambda_{n}\right)$ with $\lambda_{i}$ the $i^{th}$ singular value of matrix $X$. The principal component are obtained as the the column vectors of

The column vectors of the matrix $V$ are called empirical orthogonal functions (EOF) and they represent the weights of each principal component in the decomposition of the original normalized data of each cycle $X_{i}$, which can be approximated as

where j denotes the $j^{th}$ principal component (PC). The explained variance of each PC is proportional to square of the corresponding singular value $\lambda_{i}$. Hence the $i^{th}$ PC explains a percentage

of the variance in the data.

The moments are a set of constants that represent some important properties of the distributions. The most common are first moment and second moments, i.e. expectation value (mean value) and variance. Two other important measures are the coefficients of skewness and kurtosis, i.e. third and fourth moments. The skewness measures the degree of asymmetry of the distribution, whereas kurtosis measures the degree of flatness of the distribution. The skewness can be defined as

where $\overline{Y}$ is the mean, s is the standard deviation, i.e. square root of variance, and N is the number of data points. If skewness is positive/negative, the distribution is skewed to the right/left.

The kurtosis is defined as

Kurtosis is sometimes thought to tell about the peakedness of the distribution, but actually it tells about the tails and flatness of the distribution. If kurtosis has high/low value it has heavy/light tails (Krishnamoorthy, 2006). The skewness-kurtosis pair of normal distribution is (0,3). That is why sometimes the coefficient ’excess kurtosis’ is used, i.e. the aforementioned kurtosis-3.

Using Oulu, Hermanus and Climax cosmic-ray data, a PCA analysis was conducted by equalizing the cycles to 133 time steps (months) to get the two main principal components shown in Fig. \iref\thearticlefig:PC1_PC2s. The first and second PC explain 79.0 % and 13.3 % of the total variation of the Oulu C20-C24 data, 77.0% and 13.2% of the variation of Hermanus C20-C24 and 76.7 % and 18.8 % of the Climax C19-C22 data. Hence the two main PCs account together for 92.3 %, 90.2% and 95.5 % of the total variation of CR data for Oulu, Hermanus and Climax, respectively. Note, especially, that the PC2 explains almost one fifth of the Climax CR intensity.

Figure \iref\thearticlefig:EOF1_EOF2 shows the corresponding EOFs of the Oulu, Hermanus and Climax CR intensities. Note the almost sawtooth shape of the EOF2. The EOF2s of odd cycles are positive and EOFs of even cycles are negative, except for EOF2 of Cycle 24 of Oulu and Hermanus which have opposite sign to those for earlier even cycles. Note also that EOF2 for Hermanus Cycle 23 is slightly negative in contrast to other odd cycles. Looking the shape of PC2s of Fig. \iref\thearticlefig:PC1_PC2s, we note that a positive EOF2 means positive phase in the first half of the cycle, and a negative phase in the second half of the cycle. For negative EOF2 the phases are opposite. Note the similarity of the EOF1 and EOF2 in the cycles, which are common for all stations. EOF1s for Climax are little higher, because the PCA is calculated from only four cycles, while Oulu and Hermanus PCA are calculated from five cycles. It is essential, that the cycle 21 has lowest weight to all PC1s, but highest weight to all PC2s. (Later I call PC2 phase positive if it starts with positive phase (as in Fig. \iref\thearticlefig:PC1_PC2sb), and negative if it starts with negative phase.)

Reverting the PC1, PC2 cycles of Oulu, Hermanus and Climax CR back to their original length and concatenating them we get first and second principal component time series for the CR cycles shown in Fig. \iref\thearticlefig:PC_timeseries_firsta and b, respectively (we still use the same normalization as in the Fig. \iref\thearticlefig:Oulu_Climax_ja_muut). It is evident, that the PC1 is dominated by the sunspot cycle related period. Interestingly, the largest size of PC1 is almost equal for Cycles 19 and 22, and smallest for Cycles 20 and 24.

In contrast with PC1, the PC2 shows clear Hale cycle period (see Fig. \iref\thearticlefig:PC_timeseries_firstb). Note that the consecutive even-odd cycle pairs show similar structure for two Hale cycles. Hale cycle is traditionally referred as even-odd cycle pair (Gnevyshev and Ohl, 1948; Wilson, 1988; Makarov, 1994; Cliver, Boriakoff, and Bounar, 1996; Takalo, 2021b). However, we start here with odd Cycle 19, because this is the first cycle with cosmic-ray measurements. Note that cycle pairs 19 – 20 and 21 – 22 have opposite phases, i.e., first cycle is positive and second cycle negative. Note also a very strong peak-to-peak variation of cycle 21 PC2, and also quite strong variation of cycle 22 PC2. After that the situation changes such that Cycle 23 is still positive for Oulu as expected for the odd cycle, but slightly negative for Hermanus. Cycle 24 is again positive for both Oulu and Hermanus CR PC2 time series. There really is change from the patterns of earlier odd-even cycle pairs.

In Figure \iref\thearticlefig:PC_timeseries_resta and b we show the PC1 and PC2 timeseries of the other station we have studied. They seem very similar to the PCs of Oulu CR for Cycles 20 – 24. Note, especially, that the positive phase (seen clearly in magnification) for Cycle 23 is slightly greater for Moscow and Thule than for the other CR stations. This is also the case for Thule in cycle 24, so it seems there is no bias just in cycle 23. Moscow had weird looking cosmic-ray data for the Cycle 24 , which would have distorted the whole PCA for Moscow. That is why we omitted the cycle 24 data for Moscow.

Figure \iref\thearticlefig:CR_and_PC_proxya, b and c show the original CR intensities and PC1+PC2 proxy time series for Oulu and Climax and Hermanus stations, respectively. The main features are produced by these two main PCs quite well, especially for Climax data. The clearest differences are in the Cycle 23 for Oulu and Hermanus data. The late peaks and fast recovery after them are not brought forth by the PC1 and PC2. These features, which are unique for one cycle alone are produced by the higher order PCs. Cycle 23 has large value for EOF3 and when we add the PC3, i.e. plot also PC1+PC2+PC3, the aforementioned features are corrected quite well (green curve in Fig. \iref\thearticlefig:CR_and_PC_proxya. Note that Cycle 22 does not change at all when PC3 is added and green PC1+PC2+PC3 proxy curve is on top the red line for C22. This is because its EOF3 for the C22 is practically zero. The case is similar to Hermanus station.

In order to see if the choice of the starting minima (and the length of the cycle) of the cosmic-ray cycles affects to the PC analysis, we used the dates of the starting minima shown in Table 2. It is also a well-known fact that cosmic-ray flux lags the solar activity cycle from a few months to over a year depending on the station and cycle (Mavromichalaki et al., 1997; Gupta, Mishra, and Mishra, 2005; Rybansky, Kudela, and Minarovjech, 2009; Kane, 2014; Porta, 2018; Iskra et al., 2019; Ross and Chaplin, 2019; Koldobskiy et al., 2022). In the second analysis the PC1 and PC2 explain 78.1 % and 14.4 % of the total variation of the Oulu C20-C24 data, 74.9% and 15.3% of the variation of Hermanus C20-C24 and 74.5 % and 20.4 % of the Climax C19-C22 data. Hence the two main PCs account together almost the same part of the total variation of the CR data as in the fist analysis, i.e. 92.5 %, 90.2% and 94.9 % for Oulu, Hermanus and Climax, respectively. Figures \iref\thearticlefig:PC_timeseries_seconda and b show the corresponding time series of PC1 and PC2, respectively. In the third analysis we used still larger time-lags and these starting dates are also available in the Table 2. In the third analysis the PC1 and PC2 explain 77.3 % and 15.2 % (together 92.5 %) of the total variation of the Oulu C20-C24 data, 73.3% and 16.6% (together (89,9 %) of the variation of Hermanus C20-C24 and 78.3 % and 17.5 % (together 95.8 %) of the Climax C19-C22 data. Figures \iref\thearticlefig:PC_timeseries_secondc and d show the corresponding time series of PC1 and PC2, respectively. Note the similarity of this figure compared to Figs. \iref\thearticlefig:PC_timeseries_firsta and b, and \iref\thearticlefig:PC_timeseries_seconda and b, although the cycle starting dates now differ from 4 to 14 months from the sunspot cycle starting dates. The most striking difference in the last PC2 compared to earlier PC2s is that even Oulu has now slightly negative phase for the odd Cycle 23. This result confirms the statement that Cycle 23 is more similar to earlier even cycles than odd cycles.

Figure \iref\thearticlefig:PC1_PC2_power_spectraa and b show the power spectra of PC1 and PC2 time series of Fig. \iref\thearticlefig:PC_timeseries_firsta, and b, respectively. It is clear that PCA effectively decomposes the solar cycle and Hale cycle related periods from the CR data. According to these spectra the solar cycle period is 10.95 years (on the average) and the Hale cycle about 21.9 years. In the Fig. \iref\thearticlefig:PC1_PC2_power_spectrab, there is another smaller peak at period of one third of the Hale cycle. We believe that this peak is related to the phase change between the succeeding cycles. This causes the maxima (or minima) of the PC2 exist about 2/3 of the solar cycle (marked with two-headed arrows in Fig. \iref\thearticlefig:PC_timeseries_firstb). In the Fig. \iref\thearticlefig:PC1_PC2_power_spectrab there are two spectra of Oulu PC2. The dashed blue spectral peaks are for the whole period C20 – C24, and the solid blue spectral peaks for the case where the last Cycle 24, with ”wrong phase” is omitted from the spectral analysis. We have added power spectra of Thule and Moscow for C20-C23 to the Fig. \iref\thearticlefig:PC1_PC2_power_spectra to confirm that the peaks are the same (it is not necessary put all the power spectra, because Fig. \iref\thearticlefig:PC_timeseries_rest shows that the PC1s and PC2 are very similar to all stations). The 1/3 of Hale cycle peak is even stronger than the whole cycle peak for Oulu, Hermanus, Thule and Moscow. We suppose that this is due to strong peaks for PC2 of the Cycles 20 – 22 for all these stations. However, this peak is somewhat artificial, because it comes from the concatenation of the separate PC2s of the different cycles. The important thing is that we have removed solar cycle related period from the CR data. Note that for Climax the whole Hale cycle peak is higher than the 1/3 Hale cycle peak, because Climax CR data consist of only the distinct Hale cycle pairs 19 – 20 and 21 – 22.

We can now compare our results with the earlier studies. The positive PC2, i.e positive phase in the first half of the cycle, means that the decrease in the CR intensity is slower than if it were negative. In Figs. \iref\thearticlefig:PC_timeseries_first and \iref\thearticlefig:PC_timeseries_rest we see that this is the case for odd cycles. This result is similar to those of, e.g., Mavromichalaki et al. (1997) and Van Allen (2000). For the the even Cycles 20 and 22 the first half of PC2 has negative phase and second half positive phase. This makes decreasing of CR cycle faster and recovering more rapid leading to more peak-like structure of the even cycles. That was the result of the earlier studies too (see (Webber and Lockwood, 1988; Mavromichalaki et al., 1997; Van Allen, 2000; Thomas, Owens, and Lockwood, M., 2014)). We have shown here the essential reason for these behaviors in CR, i.e. they are caused by the Hale cycle related component of the CR data.

We have also conducted the PC analysis for sunspot number 2.0 monthly data (SSN), solar monthly corona index (MCI), and function proportional to the IMF magnetic field intensity multiplied with square of solar-wind velocity Bv². The last one has been shown to play the central role in the process of energy transfer (Ahluwalia, 2000; Verbanac et al., 2011; Takalo, 2021a). Figure \iref\thearticlefig:SSN_CI_Bv2_proxya, b and c show the PC1 and PC2 timeseries for the SSN, MCI and Bv², respectively. Notice that OMNI data, which are used to calculate function Bv² exist only since 1964, i.e. for Solar Cycles 20-24. (The PC2s have been slightly smoothed to better see the phases of the cycles). Note that PC2s are in the same phase for all solar/IMF parameters. We would assume that the solar parameters anti-correlate with cosmic-ray data. Indeed, the PC2s of the cycles 21 and 22 (marked as gray) are exactly in opposite phase with the PC2 of the cosmic-ray data. Note that PC2 is by far strongest just for those cycles in CR data. The solar wind speed and consequently function Bv² seem to be stronger always in the latter half of the cycle, i.e. during the descending phase of the Solar Cycle (see \iref\thearticlefig:SSN_CI_Bv2_proxyc) (Takalo, 2021a). We would then suppose that cosmic-ray has negative phase as preferred state in the second half of the cycle. This is actually the case for cycles 23 and 24 (instead slightly for Hermanus C23 and in one case for Oulu). Although we do not have solar wind data for the cycle 19, the situation could be the same for this cycle also. The cycle 20 is the only one which does not fit in to the aforementioned pattern. Note, however, that the Bv² PC1 is totally different for the Cycle 20 than the other cycles. The PC1 for Cycle 20 is almost constant and the stronger second half rather exists in the PC2 for this cycle. Agarwal and Mishra (2008) already found that the interplanetary magnetic field strength (B) shows only a weak negative correlation (-0.35) with cosmic-rays for the solar cycle 20, whereas it shows a high anti-correlation for the solar cycles 21 – 23 (-0.76,-0.69). They suspect that a significant contribution to modulation from the termination shock during solar cycle 20 could dilute the correlation of cosmic-rays with the interplanetary magnetic field B at one astronomical unit for that cycle. The perturbations of the heliosphere are weaker and less widely spread during solar cycle 20 than during other solar cycles. This might lead to a situation where the heliospheric perturbations are relatively small for cosmic-ray particles allowing these particles to reach the Earth as if it was a minimum solar activity period. Note from figure \iref\thearticlefig:PC_timeseries_firsta that the decrease of the cosmic-ray intensity is smaller for Cycle 20 than for Cycle 23, especially for Hermanus, although sunspot numbers are almost the same for these cycles and Solar corona even stronger for Cycle 20 than for Cycle 23 (see Fig. \iref\thearticlefig:SSN_CI_Bv2_proxy).

We have calculated skewnesses and kurtoses of CR cycles 19-22 for Climax and cycles 20-24 for Hermanus, Kerguelen, McMurdo, Newark, Oulu and Thule. Figure \iref\thearticlefig:Skewness_kurtosis shows the result in skewness-kurtosis (S – K) coordinate system. There are several interesting features. The cycles other than C21 are near the regression line (shown as black dashed line) such that cycles 23 and 24 are in the left lower corner meaning small skewness and low kurtosis. Interestingly Cycle 19 for Climax is among the cycles 23 and 24 in the left lower ellipse. The even cycles C20 and C22 are in the right upper corner with larger skewnesses and higher kurtoses. The correlation coefficient of the regression line is 0.85 and 95% confidence limits are shown as red dashed lines. Note that only even Cycles 20 and 22 for Hermanus are somewhat outside the 95% limits. The Cycles 21 are omitted from this regression and they are located compactly inside a small ellipse with small skewness and moderate kurtosis. We should remember that the Cycle 21 differed from other in PCA such that its weight for PC1 is smallest and PC2 highest of all cycles. Another classification could be that all odd cycles plus cycle 24 have small skewnesses and even cycles 20 and 22 have moderate to large skewnesses. The skewness-kurtosis map also shows that Cycle 21 differs from others such that they are compactly located in a small area in this coordinate system. Also the skewness-kurtosis values for another odd Cycle 23 form a quite compact group, but with somewhat lower kurtoses and slightly larger skewnesses. The interesting thing is that the only Cycle 19 has almost smallest value for skewness and kurtosis of all cycles. It is a pity that we do not have any other results for the whole cycle 19, and even this data was measured with the old IGY-type monitor. There exists another whole Cycle 19 measured at Huancayo, and we have marked it as an ’x’ to the skewness-kurtosis coordinate map. Its skewness is almost zero and it locates near the cycles 21 in the S – K map. The other cycles for Huancayo C20 and C21 are, however, so flat, i.e., have so heavy tails, that their kurtosis is outside the S – K map of Fig. \iref\thearticlefig:Skewness_kurtosis. The reason is probably very high rigidity cutoff for that station, 13.4 Gev, and thus its insensitivity to low energy cosmic-rays (Usoskin et al., 2001, 2017). Another reason to omit Huancayo in PCA is that there exists only three whole cycles for this station.