Abstract
The time variability of the cosmic ray (CR) intensity, recorded by the Climax neutron monitor and covering the period 1953–2004, has been analyzed by the joint application of the wavelet and the empirical mode decomposition (EMD) analyses. Dominant time scales of variability are found at ~11 yr, ~22 yr, ~6 yr and in the range of the quasi-biennial oscillations (QBOs). The combination of the 11 yr cycle and QBOs explains the Gnevychev Gap (GG) phenomenon and many step-like decreases characterizing the CR modulation. The additional scales of variability at ~22 yr and ~6 yr are responsible for other features of the long-term CR trend, such as the intensity flat-topped profile, following the maxima of even-numbered cycles during positive polarity state of the heliosphere (). Comparison with basic time scales of variability derived from the sunspot area (SA) allows the association of the 11 yr cycle and QBOs with solar activity variations, whereas the other two modes with the drift effects govern the CR entrance in the heliosphere.
1. Introduction
It has been known for a long time that galactic CRs are modulated by solar activity (see the pioneering works performed by Forbush [1, 2]). In hundred years of research from the CR discovery [3] the topic of cosmic ray modulation has been widely investigated. During the sixties and seventies, the different forms of modulation (short-, medium-, and long-) were identified by a deep analysis of the available experimental data. Solar and interplanetary structures were recognized as the sources of the 3D variability of the incoming charged particles in the heliosphere and several mathematical techniques were applied on the different datasets during the eighties and nineties to improve the knowledge of the solar-heliospheric domain (see, e.g., [4]). Among them, the Wavelet Transform (WT) was widely recognized as a good tool to investigate the time evolution of periodic or quasi-periodic trends in the experimental data sets [5, 6]. The WT, by decomposing a time series into the time-frequency space, allows to analyze processes containing multiscale features and it is able to determine both the dominant variability modes and how those modes vary in time. The possibility to evaluate a local spectrum distinguishes the use of the wavelet analysis from, for example, Fourier approaches or other mathematical filtering techniques used in the past. Since this technique provides information both in frequency and time domains it allowed to characterize the evolution of CR variability at different scales [7, 8], to identify and characterize the contribution of the different periodicities in the cosmic ray intensity profile for distinct solar cycles and/or different phases of each cycle [9, 10], and to correlate variations, detected in solar cycle indices, with the corresponding fluctuations in CR intensity [11, 12]. Finally the wavelet represented a very useful tool to investigate transitions between different periodicities and to identify a clear difference in the solar activity evolution during odd and even solar activity cycles [9, 10, 13, 14].
More recently, the Empirical Mode Decomposition (EMD) analysis, developed as a signal processing for nonlinear and nonstationary data [15, 16], was also introduced in the study of the solar-heliospheric domain (e.g., [17–21]). The method allows to decompose any complicated dataset into a finite number of components (called Intrinsic Mode Functions: IMF), representing simple oscillatory modes in the time domain with the signal described by a time-dependent amplitude and phase functions.
In this paper we analyze the cosmic ray variations together with sunspot areas recorded for the period from 1953 to 2004. In Section 2 we present results obtained by applying the WT and the EMD to the data set, showing that the EMD modes add more information with respect to the WT only. Section 3 describes two reconstructions of the CR signal, by using different selected EMD modes and discuss the possible association with the different sources of the CR variability, namely solar activity and drift effects. In Section 4 conclusions are drawn.
2. Data Used and Method of Investigation
In order to study the behavior of the CR intensity, monthly mean values derived from records of the Climax neutron monitor (Geographic latitude-longitude: N39.37°-E253.82°; cutoff rigidity: ~3 GV; height: 3400 m a.s.l.; http://ulysses.sr.unh.edu/NeutronMonitor/) have been analyzed for the period 1953–2004. Moreover, monthly averages of sunspot area (SA; http://solarscience.msfc.nasa.gov/greenwch.shtml) have also been considered as a parameter of solar activity responsible for the CR modulation. Search for periodicities in the above data sets is performed by using the WT and the EMD techniques.
As stated in the introduction, the WT analysis offers a significant advantage with respect to the classical Fourier transformation because it allows to localize in time possible periodicities, even those not always present through the whole considered time interval. As mother function we assume the Morlet one, which supplies a better resolution in frequency compared to other available functions. The scales have been chosen as fractional power of 2: The following parameters have been used: , is the frequency resolution, months, is the smallest resolvable scale, and determines the largest scale. Since the set of scales is nonlinear, the associate error is asymmetric. The left and right uncertainties on each wavelet period can be derived as by using (1) and the relation , valid for the Morlet wavelet.
The Wavelet Power Spectrum (WPS), as a function of time and frequency, has been computed in order to characterize the frequencies of the main periodicities. The significance of the power peaks is evaluated against a white noise background [22]. We remark that the wavelet analysis is based on a priori definition of the mother functions, which could not be representative of the eigenfunctions of the investigated phenomenon. For this reason, the EMD, which is an adaptive technique, is more appropriate to deal with nonstationary data [15]. In the EMD framework, a signal can be decomposed into a finite number of intrinsic mode functions (IMFs) as The mode is a zero-mean oscillation, variable in both amplitude and frequency, which can be expressed as , where and represent the instantaneous amplitude and phase of the th EMD mode, respectively. Each IMF has its own characteristic period , defined as the average time difference between local extrema (local maxima and minima). The standard deviation of the differences between local extrema quantifies the uncertainty associated with each period. This kind of decomposition is local, complete, and orthogonal [15, 23] and the residue in (3) describes the mean trend. The statistical significance of each IMF, with respect to white noise, can be checked by applying a specific test based on the following argument [24]. When EMD is applied to a white noise series, the constancy of the product between the energy density of each IMF and its corresponding averaged period can be deduced. This relation can be used to derive the analytical energy-density spread function of each IMF as a function of different confidence levels. This approach represents the analogy of the statistical significance test used to evaluate the significance levels in the Wavelet power spectrum [22]. Thus, by comparing the energy density of the IMFs extracted from the actual data with the theoretical spread function, IMFs containing information at the selected confidence level can be distinguished from purely noisy modes. Finally, the orthogonality of EMD modes allows to reconstruct the signal at a chosen time scale through partial sums in (3) (see also [17–19]).
We take advantage of both WT and EMD for a detailed study of the CR variability. First, we determine the characteristic periodicities of the signal through the wavelet analysis, which provide frequencies equivalent to the Fourier’s ones. Then, we investigate the time variability of all the basic modes, obtained through the EMD, that contribute to the periodicities identified in the WPS. Figure 1 shows the WPS derived from CLI data. The highest values of power are found, as expected, at yr and yr (where the right and left errors are reported as superscript and subscript, resp.). These periodicities are related to the Schwabe cycle and to the polarity inversion of the Sun’s magnetic field (Hale cycle), respectively. Nevertheless, the former cannot be well isolated, because high power is detected in a broad range of periods from 5 to 15 yr; the latter is not reliable, being under the cone of influence, determined by the short length of the data set. Enhanced power is also visible between 1 and 4 yr, but it is significant only during the maximum phases of cycles 21 and 22, being above the confidence level of 80% with respect to white noise.
In order to better characterize the oscillations in the frequency ranges identified through the wavelet analysis, we decomposed the CLI data into a set of IMFs. Ten modes were obtained and displayed in Figure 2 along with their characteristic periods. The energy density of each IMF (, where denotes time averages) is shown in Figure 3 as a function of its corresponding characteristic period, to understand which modes dominate the CR variability. The highest amplitude mode , having yr, can be associated with the Schwabe cycle, whereas , having yr, possibly to the Hale cycle, because its period is comparable with ~22 yr within the error.
The IMF ( yr) is found to be the third important mode, although it was not separated from the ~11 yr periodicity when using the wavelet analysis. Other modes with comparable amplitude are , having periods yr, yr, , respectively, that is, time scales of the so-called quasi-biennial oscillations (QBOs [25–27]). Mode , with yr, presents only one oscillation which cannot be considered as reliable, given the limited time extent of the data set. Finally, modes with yr have the lowest amplitudes. Moreover, Figure 3 shows that all the computed IMFs are significant, their amplitude being located above the spread lines (dashed, dot-dashed, and dotted lines at 90th, 95th, and 99th confidence level, resp.), obtained by performing the significance test described above.
In order to better understand the links between the identified periodicities and the solar activity, we also performed the EMD analysis on the SA data. The derived modes are shown in Figure 4. Four IMFs (, having characteristic periods yr) were found to be not significant at least at 90th confidence level (see Figure 5). Two significant modes () have characteristic periods in the QBO range ( yr and yr, resp.), whereas the IMF with yr has the highest energy (as shown in Figure 5) and it is representative of the Schwabe cycle. Finally, the IMFs (having periods yr and yr, resp.) could be both related to the shape of 11 yr cycle, possibly related to the Gnevyshev-Ohl rule. We underline that no mode at ~6 yr is detected in SA. This indicates that the IMF , detected in CLI data, could not have a solar origin, as will be discussed in Section 3.2.
3. Reconstruction of the CR Signal
3.1. CR Modulation Driven by Solar Activity
An important source of CR modulation is represented by solar activity, which is variable on a wide range of temporal scales. The main modulation is at 11 yr time scale and it is related to the well known Schwabe cycle. As a consequence, the CR intensity shows a pronounced 11 yr variation, with the intensity maximum corresponding almost with the solar minimum and vice versa. In addition, many manifestations of solar magnetism as well as some interplanetary phenomena show quasi-biennial variations. The QBOs have been found in photospheric magnetic field and in many solar activity indices such as the sunspot number and area, the number of flares [26], the 10.7 cm radio emission [12], and the coronal emission ([28], and references therein). In the interplanetary medium, QBOs have been observed in the solar wind speed and the interplanetary magnetic field intensity (e.g., [29]), in the number and flux of interplanetary energetic protons at Earth’s orbit [8, 27], in the density of galactic CRs at 1 AU [9, 14], and in the outer heliosphere [10].
It has been shown that the effect of solar activity on CR flux can be essentially described by superposing QBOs to the ~11 yr mode (e.g., [18, 21]). Thus, in order to investigate the CR modulation driven by the solar activity we perform the reconstruction of both the CR and SA signals, by considering QBOs and the ~11 yr mode. As the wavelet analysis indicates that QBOs are located in the Fourier equivalent period 1–4 yr, we reconstructed the signal at QBOs scales through partial sums of IMFs as in (3) considering only the significant modes () having characteristic periods in this range (, , , resp.). This is also consistent with the criterion used by Laurenza et al., [21] for the selection of clear and unambiguous EMD components contributing to the QBOs. Note in Figure 2 that is generally the dominant mode for the CLI QBOs, having the highest amplitude during most of the time. Nevertheless, mode has a smaller amplitude than during 1988–1992, although they are almost in phase with each other. In addition, the amplitude is higher only during the period 1953–1967 and comparable with during the period 1970–1973. The full QBO signal is depicted in Figure 6 and has a high amplitude during the maximum phase of each solar cycle as expected from results of Bazilevskaya et al. [26], whereas they almost vanish during the sunspot minima.
The importance of QBOs with respect to the ~11 yr, around a solar maximum has been evaluated by computing the ratio between the mean square amplitude of the QBOs and that of the ~11 yr mode. Being the superposition of the modes , representing the contribution of the QBO, the ratio quantifies the importance of QBOs with respect to the ~11 yr mode (being the average computed over a time interval of 5 yr around each maximum of the 11 yr cycle). The QBO contribution is estimated to be higher than 30% with respect to the amplitude of the ~11 yr mode for all the considered cycles (34%, 32%, 38%, 54%, and 31% for cycles 19, 20, 21, 22, and 23, resp., in agreement with the findings of Laurenza et al. [21]).
Then, we compare the superposition of the QBOs and the ~11 yr mode with actual CR data. It is apparent in the upper panel of Figure 7 that the superposition of the ~11 yr and QBOs well reproduces the general trend of CR data, mainly during the maximum and descending phase of sunspot cycles; that is when QBOs reach their highest amplitude. In particular, this superposition accounts for the majority of the step-like decreases, which characterize the CR modulation. Notice that the so-called mini-cycle during 1974 is reproduced. The cosmic ray depression during the mini-cycle was related [30] to a decrease in the radial diffusion coefficient (leading to a reduced efficiency of cosmic ray particles to enter the inner solar system). Under the assumption that scales inversely with as , the depression was attributed to increases in the field magnitude through a diffusion dominated process [31], also consistently with findings of Laurenza et al. [21].
(a)
(b)
The GG phenomenon, namely the peculiar decrease in the temporal trend of various solar-terrestrial indices ([32] and references therein) during the maximum phase of the Schwabe cycles (for early works see [33, 34]), is also observed to be generated by adding QBOs and the ~11 yr mode.
The short term variations of the CR intensity, not explained by the QBOs, could be due to the so called merged interaction regions [35], acting as a barrier to cosmic rays [36]. In fact, solar wind transient flows may interact with one another and with corotating streams to form a large-scale interaction region with a shell-like geometry (global merged interaction region—GMIR) that persists over several solar rotations. These time scales are comparable with the high frequency modes ( yr) we detected in CR data, although their possible connection with GMIR has to be further investigated.
We remark that, around sunspot minima, other low frequency modes would be necessary to explain the variability of the CR flux more related to drift effects in the large scale field of the heliosphere that will be treated in the following subsection.
Finally, we study the relationship between medium and long term CR variations and sunspot area. We describe the temporal variability of SA by superposing the ~11 yr mode and QBOs ( and are the contributing IMFs), which are of the order of 50% with respect to the ~11 yr component during the maximum phase. Figure 8 shows the reconstruction for both CLI and SA, which present almost the same profile and are shifted between each other. In particular, the GG is well detected in CR for all solar cycles and delayed, with respect to SA, of an average value of about 6 months. Similar results were obtained in the past (e.g., [37, 38]) in correlative studies between CR data and solar activity indicators. As SA are a proxy for solar activity causing interplanetary disturbances, which act on CR as propagating diffusion barriers, the found delay can be explained in terms of the propagation time of magnetic perturbations from the Sun towards the edge of CR modulation region. Simple estimations of the propagation time can be obtained for distances of 80 AU () and 100 AU (), where the cosmic ray modulation should be still important [39, 40]. In fact, by assuming the average solar wind speed km/s, the CME mean speed during solar maximum km/s [41] and the velocity of the high speed streams km/s as characteristic velocities of three different solar wind regimes, the corresponding propagation time are months/ months, months/ months, months/ months. Indeed, such values are also consistent with the delay time between transient CR decreases measured at different heliocentric distances, following intense solar-interplanetary events. For instance, a large propagating transient CR decrease was detected at neutron monitor energies at the Earth in October 2003 (the Halloween event; Plainaki et al. [42]) and about 6–8 months later by CR detectors aboard Voyager 2 and Voyager 1 at 73 AU and 92 AU, respectively [43, 44].
3.2. CR Modulation and Drift Effects
It is widely known that the CR intensity curve follows a 22 yr cycle with alternate maxima being flat-topped and peaked (e.g., [45]). This peculiar behavior is described by models of CR modulation [46–52], based on the observed reversal of the Sun’s magnetic field polarity, curvature and gradient drifts in the interplanetary magnetic field. In the drift formalism, during epochs when (i.e., when the Sun’s dipole magnetic field is positive in the northern hemisphere) the approach in the inner heliosphere of positively charged CRs is from over the poles, while during epochs when the preferred direction is from along the heliospheric current sheet (HCS). The average excursion, during a solar rotation, of the HCS from the heliographic equator is referred to as tilt angle (e.g., [53]), which characterize the drift effects on CRs (e.g., [50, 54, 55]). In periods of minimum solar activity, that is, when the CR modulation by solar activity is reduced, the effects related to the polarity state of the heliosphere and to the HCS tilt are supposed to be more important. In order to account for the CR behavior described above, modes have to be considered. As clearly seen in the lower panel of Figure 7, the superposition of these IMFs to the ~11 yr one is necessary to reproduce the long term variation of the actual CR data during the minimum solar activity periods and the maxima of odd numbered cycles. Hence, the combination of both IMFs is fundamental to take into account the full contribution of the drifts to the CR modulation. In fact, should be representative of the heliosphere polarity changes, its characteristic average period being ~22 yr. In addition, the mode seems to be unrelated to solar activity perturbations, because no comparable periodicity was found in the sunspot area, as outlined in Section 2. We suggest that a possible reason for the CR variability at this time scale can be the particle drift related to the variation during the solar cycle of the latitude extent of the HCS. As a matter of fact, the characteristic period yr is a reasonable average time for the latitudinal excursion of the HCS (although the HCS decrease time during the descending phase of the sunspot cycle is usually higher than the rise time, especially for odd numbered cycles, e.g., [56]). Moreover, the mode plays a crucial role in determining the flat-topped maximum of CRs during even sunspot cycles, in agreement with the expectations of drift models that include the changing of HCS inclination during the 11 yr solar cycle (e.g., [50, 55]). Nevertheless, a better understanding of the nature of the mode is needed and will be faced in a future paper.
4. Conclusions
A new approach, through the combined application of two powerful techniques, such as the wavelet and EMD, was performed on CR data from the Climax neutron monitor to investigate the time variability and the relationship with the variations of the solar activity, as parametrized by the sunspot area. Our results provide the first evidence for the different causes producing the CR modulation, identified according to the temporal scales characterizing the CR variability. In particular, we were able to discern the CR variations directly induced by the solar activity and those possibly related to the particle drift in the large scale magnetic field in the heliosphere. Solar activity is found to be responsible for the GG phenomenon and for most of the step-like CR modulation, through the combined effect of the 11 yr cycle and QBOs.
The effect of the drifts has been clearly identified and associated with the CR variations at ~22 yr and ~6 yr. In particular, time variations at ~22 yr are due to the polarity state of the heliosphere, whereas those at ~6 yr could reflect the latitudinal excursion of the HCS. Both contributions explain the CR time profile, as expected from the drift theory and models, during the descending and minimum phase of solar activity of even numbered solar cycles (compare top and bottom panel of Figure 7), that is when and the CRs enter the heliosphere from the poles, while the HCS tilt decays rapidly to low angles, leading to a fast recovery of the CR intensity soon after the solar maximum. In addition, the shape of the CR intensity during the ascending and maximum phase of the odd numbered cycles (again in periods when ) is determined by these variations as well.
cknowledgments
This work was partially supported by the ASI/INAF Contract no. I/022/10/0, by the European Social Fund, European Commission, and by Regione Calabria. Thanks are due also to the Italian PNRA for the use of the RAC-ANT database, prepared at the former IFSI-Roma (now IAPS) under Contract PROGDEF09_37 (code 2009/A3.07).