#### Abstract

Cosmogenic isotopes have frequently been employed as proxies of ancient cosmic ray fluxes. On the basis of periodicities of the ^{10}Be time series (using data from both the South and North Poles) and the ^{14}C time series (with data from Intercal-98), we offer evidence of the existence of cosmic ray fluctuations with a periodicity of around 30 years. Results were obtained by using the wavelet transformation spectral technique, signal reconstruction by autoregressive spectral analysis (ARMA), and the Lomb-Scargle periodogram method. This 30-year periodicity seems to be significant in nature because several solar and climatic indexes exhibit the same modulation, which may indicate that the 30-year frequency of cosmic rays is probably a modulator agent for terrestrial phenomena, reflecting the control source, namely, solar activity.

#### 1. Introduction

The importance of cosmic ray variations was pointed out long ago in a vast compendium of relevant research [1]. Here we encounter a periodicity that has not been studied within the framework of cosmic ray variations. Let us begin by emphasizing that 30-years cycles are quite common in nature: the so-called Markowitz wave (Markowitz wobble-MW) is a quasi-harmonic variation of the middle pole of the Earth with a period of 30 year and an amplitude of 0.02′′-0.03′′ [2]. Similar results were obtained in different years by many researchers who were trying to measure the Earth’s magnetic field [3–5], by evaluating the conductivity of the lower mantle and who established 60- and 30-year variations of the geomagnetic field of the lower mantle. The authors in [6], through studies of high-growth anomalies of the secular variation with foci located in South Asia and in the middle of the Indian Ocean, revealed that the increase in the current focus is the initial stage of the 30-year and 60-year variations.

Authors in [4] by evaluating the conductivity of the lower mantle have found evidence of 60 - and 30-year variations in the geomagnetic field. By conducting a qualitative analysis of timelines with data from planetary indices Ap/aa, interplanetary magnetic field intensity and sunspot numbers, as well as cosmic rays data strings taken from 1937 to 2010 [7, 8] have shown quasi-periodic three-cycle trends, which are particularly very well correlated with solar wind, polar coronal holes and the size of solar activity cycle 23 [9, 10].

Authors in [11] have found 32-year variations in temperature by analyzing the spectrum of periodical air-surface temperature fluctuations for 1423 years in Greenland ice cores, and for 1400 and 800 years in the California Arctic pine tree rings. A 30-year variation in storminess data was found by applying the *“Caterpillar”* method in [12]. In fact, the authors identified frequencies at 90 to 100 years, 28 to 32 years, 20 to 22 years, 9 to 13 years, and some others; it is shown that even if the accuracy of this method for determining periods is not high, the fundamental frequency obtained by the *“Caterpillar”* is real.

The application of wavelet analysis to the paleoclimatic proxy data [13] to large-scale atmospheric phenomena (the Atlantic Multidecadal Oscillation and Southern Oscillation Indexes), and hurricane phenomena, led to the discovery of a high coherence with periods of years, between climatic oscillations and cosmic rays (^{10}Be at the North Pole). Furthermore, some properties of hurricanes, such as their total cyclonal energy and the tropical storms appearance along the Atlantic coast of México, together with other properties linked to hurricanes show such a 30-year cycle [14]. Figure 1 illustrates the particular case of category-4 hurricanes

Curiously, not only in nature but in several areas of contemporary human activity, such as business and commerce, 30-year cycle indexes are also found. (e.g., http://www.search.yahoo.com/search?p=30+years+cycle&ei=UTF-8&fr=moz35/). Incidentally, they are often associated with solar activity.

In the present work, we present evidence of the existence of a thirty-year periodicity for the cosmic rays: preliminary results were presented at the 2008COSPAR meeting in Montreal.

#### 2. Data and the Spectral Wavelet Analysis

The spectral analysis of cosmic ray data from neutron monitors has been widely studied through several different methods, for instance [1, 15–18] and references included. Though most of these studies are out of the scope of this work, since they rather concern short-term periodicities, they are, however, very helpful for solar-terrestrial physics.

One of the main problems in determining significant long-term periodicities in the flux of cosmic rays is that the time series of data are relatively very short, as they have been available only for the last five to six decades, when data on cosmic rays (CR) from different stations throughout the world began to be organized and homologated. Because of this restriction, a proxy for cosmic rays has often been employed, one of which in our case allows us in a deterministic way to provide evidence of a 30-year cycle for cosmic rays.

The cosmogenic isotopes Beryllium-10 (^{10}Be) and Carbon-14 (^{14}C) are conventionally considered to be proxies for cosmic rays, in such a way that an adequate spectral analysis may reveal important periodicities. These cosmogenic isotopes are mainly produced by galactic cosmic ray flux modulated by changes in interplanetary and geomagnetic magnetic fields. The analysis of cosmogenic isotopes stored in natural archives, such as ^{10}Be in polar ice cores and ^{14}C in tree rings, provides a means of extending our knowledge of solar variability over much longer periods (e.g., [19, 20]). In addition, the nature of climatic response to solar variability can be assessed over several time scales. It should be remembered that the analysis of the cosmogenic isotopes record is more difficult than the analysis of sunspot numbers. This is due to the fact the ^{14}C and ^{10}Be concentrations reflect the production rate, which is modulated by not only solar activity but also by atmospheric transport and deposition processes [21, 22]

Data on ^{10}Be and ^{14}C can be obtained for periods of thousands of years: we use the INTCAL 98 (http://www.depts.washington.edu/qil/) for the ^{14}C time series and the ^{10}Be time series from [21], which give the concentration found in the Dye-3 ice core (62.5 N, 43.8 W). For the South Pole we used data from [23].

The spectral techniques for analyzing periodicities of cosmophysical phenomena are very varied. The simplest technique for investigating periodicities is the Fourier Transform (FT). Although useful for stationary time series, this method is not appropriate for time series that do not fulfill the steady state condition, as is the case with cosmogenic isotopes.

In order to find the time evolution of the main frequencies of the time series, we apply the wavelet method using the Morlet mother wavelet ([24–26]. Wavelet analysis can be used for analyzing localized variations of power within a given time series at many different frequencies. However, even using this powerful tool, the reconstruction (filtering) of a signal is one of the main problems in the field, as is well known in Signals Theory and signal processing; this follows from the fact that the use of indiscriminate algorithms may lead to findings of spurious nonexistent periodicities in the time series. Furthermore, some real existing frequencies may be masked by more prominent frequencies. To avoid these problems, the Daubechies algorithm [27] has been used. This has proven to be highly efficient for the decomposition of signals in low and high frequencies, with the advantage that it does not create fictitious periodicities in the time series and, in some cases, may be more powerful than the multianalysis wavelet procedure (see Appendix A).

In Figures 2 to 7 we present the results of the wavelet analysis for the time series: the time series studied; the time series themselves are shown in the upper panel. The wavelet Morlet spectrum for the series appears in the mid panel of every figure, and the global wavelet spectrum appears in the right-hand part of the figures, where the dashed line indicates a border of 95% reliability.

It can be seen in Figure 2 that the 30-year periodicity has a confidence of 95%; however, it looks less prominent relative to the 60, 120, and 240 years frequencies, which show confidences far above 95%.

Figure 3 shows the Wavelet of ^{14}C; it can be seen that precisely at the 30-yrs. frequency there appears a small hump relative to the importance of the other periodicities. It would then be natural to ask how to know if such periodicity really does exist with a good reliability.

This turns out to be a very complex problem, to give relevance to the periodicity of 30 years; on the other hand, it is necessary for one side to filter that signal and, on the other hand, to eliminate masking frequencies higher or lower than 30 years. It is precisely in such a situation that the *Daubechies* algorithm [27] turns out to be a very powerful tool, as shown in Figures 4 and 5 where, after the filtering process, the 30-year frequency for ^{10}Be and ^{14}C at the North Pole appears very clearly, both with high reliability, far above 95%.

The way to discern whether the periodicity of 30 years found in the ^{10}Be time series really reflects a comic ray fluctuation or is merely a local phenomenon of the cosmogenic isotope at earth level is to compare the behavior of ^{10}Be at both the North and South Poles. Since concentrations are quite different from one pole to another, it should be expected that their behavior would also be quite different if the existence of such periodicity is a local phenomenon. However, an examination of Figure 6 indicates that the wavelet coherence between the ^{10}Be at both the North and South Poles is very high, >0.9 (red color in the color bar of Figure 1)) at their common frequencies. This is an alternative way to confirm other arguments published in the literature, that this particular cosmogenic isotope is a proxy of cosmic rays.

Cosmic ray variations are mostly caused by time variations in the interplanetary magnetic field (IMF), so periodicities in cosmic ray data must be reflected in IMF data. Unfortunately, confident data from the IMF only date from the beginning of the spacecraft era, not even two 30-year cycles. However, since presumably the main source of such modulation is found in solar activity (SA); to confirm such an hypothesis we carried out a wavelet analysis of SA by using a time series of sunspot data (ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SUNSPOT_NUMBERS/ and also http://sidc.oma.be/sunspot-data/). This analysis allows us once again to find the existence of such a 30-year frequency (Figure 7), where the results are presented without filters.

Since the 30-year periodicity appears as a small peak, (though one with higher than 95% confidence), instead of applying the Daubechies filter [27], to confirm that its existence is real, we have applied collateral methods that confirm such a periodicity for cosmic rays. In the next section the periodicity of 30 years of solar activity is clearly shown, which in the unfiltered spectrum of Figure 7 is just barely perceptible.

#### 3. Autoregressive Spectral Analysis ARMA

To verify the results obtained, Libin and Yudakin have calculated the mutual power spectra and coherence spectra for solar activity and temperature (Figure 8), solar activity and storminess (Figure 9), solar activity and the level of Lake Baykal (Figure 10) and, finally, solar activity and cosmic ray intensity from neutron monitor data and measurements of ^{10}Be (Figure 11) for different measurement periods (see the data in the graphs). Calculations were made using the spectra of autoregressive spectral analysis with simultaneous ARMA (see Appendix B) and filtering the input data with the suppression of 11-year and 22-year variations.

All the above figures show the presence of stable 30-year fluctuations for almost all the processes. Although the coefficients of coherence are not always superior to a 95% confidence interval probability the observed peaks are almost always above 90%.

#### 4. The Standard Lomb-Scargle Periodogram Method

The author in [28] using the Lomb-Scargle technique has built periodograms based on the long-time series data: Figure 12 shows the periodogram, based on data from direct monthly stratospheric measurements of CR provided routinely over a long time period [29, 30]:

It can be seen that an approximately 34-year cycle is present in the data. A similar indication was obtained from the data of a modulation parameter of CR from [31], as shown in Figure 13.

#### 5. Conclusions

We have shown the existence of a 30-year periodicity for cosmic rays. We can argue here that such a frequency is quite probably a modulator for terrestrial phenomena: it seems that in some way cosmic rays modulate climatic phenomena, such as the Atlantic Multidecal Oscillation (AMO) and sea-surface temperature (SST), and these, in turn, modulate hurricane development [31, 32]. Furthermore, a wavelet analysis applied to paleoclimatic proxy data for large scale atmospheric phenomena (Atlantic Multidecadal Oscillation and Southern Oscillation Indexes) has revealed coherence between climatic oscillations and cosmic rays on a 30-yrs cycle [14]. Since this periodicity is present throughout the entire interval under study, the origin of such a periodicity may be associated with the 120-year secular cycle of solar activity whose presence has been demonstrated in [17]. It corresponds to continuous periods of increasing and decreasing activity during maxima and minima of that secular cycle. The 60- and 240-year. solar activity cycles may, then, also be associatedwith the 30-year cycle. This would confirm that *Solar Activity is the source of all cosmophysical modulators of Solar-Terrestrial relationships*, through e a number of intermediaries at short, medium and long frequencies. One of them, in particular, is the 30-years periodicity found in this work, which allows for the analysis of a reasonable long-term variability for Cosmic rays.

#### Appendices

#### A.

*Inverse Wavelet Transform*

We use the inverse wavelet transform to obtain the decomposition of a signal which can be obtained from a time-scale filter [27]. The inverse wavelet transform is defined [25] as
where is the factor for scale averaging, is a constant ( and , for the Morlet wavelet), and removes the energy scaling. We use the inverse wavelet transform to obtain the time series.

#### B.

*Autoregressive Moving-Average Model*

The autoregressive moving-average model (ARMA) is one of the mathematical models used for the analysis and prediction of stationary time series in statistics. The ARMA model is a generalization of two simpler time series models—an autoregressive model (AR) and the moving average model (MA).

The ARMA (, ) model, where and are integers that specify the order of the model is called the next generation process time series :
where is constant, represents white noise, and and are real numbers, the coefficients of the autoregressive, and moving average coefficients, respectively.

Such a model can be interpreted as a linear multiple regression model, in which the explanatory variables are the past values of the dependent variable itself, but as a regression balance—moving averages of the elements of white noise. ARMA-processes are more complex compared to similar processes in a pure form; however the ARMA processes are characterized by fewer parameters, which is one of their advantages.

If we introduce the lag operator , then the ARMA model can be written as follows:

Introducing the shorthand notation for polynomials of the left and right sides, the previous equation can be written as

For the process to be stationary, it is necessary for the roots of the characteristic polynomial of the autoregressive part to lie outside the unit circle in the complex plane (in modulus strictly greater than one). The stationary ARMA process can be represented as an infinite MA process:

For example, the process ARMA (1,0) = AR (1) can be represented as an MA process of infinite order with coefficients in decreasing geometric progression:

Thus, the ARMA processes can be considered to be MA processes of infinite order with certain restrictions on the structure coefficients. There is a small number of parameters to describe the processes they enable rather than a complex structure. All stationary processes can be arbitrarily approximated by an ARMA model of a certain order with considerably fewer parameters than MA models use.

*NonStationary (Integrated) ARMA*

In the presence of unit roots of the p autoregressive polynomial, the process is nonstationary. Roots of less than unity in practice are not considered, since they are processes which exhibit explosive behavior. Accordingly, to test the stationary nature of a time series of basic tests, tests must be run for unit roots. If the tests confirm the presence of unit roots, then we need to analyze the difference between the original time series and a stationary process of the differences of one or two orders (usually the first order is sufficient and sometimes the second) of the ARMA-based model.

Such models are called ARIMA models (integrated ARMA) or Box-Jenkins models. The ARIMA model (, , ), where d is the order of integration (the order of differences in the original time series); and , the order of AR; MA the parts of the ARMA-process differences , the order can be written in the operator form:

The ARIMA process (, , ) is equivalent to the ARMA process (, ) with unit roots.

To construct the ARMA model on a proxy data series of observations, it is necessary to determine the model order (numbers and ), and then the coefficients themselves.

To determine the order of the model an investigation of these characteristics of the time series can be done, seen as its autocorrelation function and partial autocorrelation function.

To determine the coefficients the method of least squares and maximum likelihood method can be used.

*ARMAX Models*

In the classic ARMA model, it can add exogenous factors x. In general, the model involves not only the current values of these factors but also lagged values. Such models are usually denoted ARMAX (, , ), where k-lags come from the exogenous factors. In an operator form, such models can be written as follows (an exogenous factor):
where , , are the order polynomial, respectively, , , of the lag operator.

It should be noted that such models can be interpreted differently, for example, ADL () mode with random errors MA ().

#### Acknowledgments

V. Velasco gives special thanks to the CONACyT (project 180148) for financial support to this work.