Short-Term Periodicities Observed in Neutron Monitor Counting Rates

Neutron monitor counting rates and solar wind velocity, interplanetary magnetic field, sunspot number and total solar irradiance measurements from 2013 to 2018 corresponding to the end of solar maximum and the decreasing phase of the Solar Cycle 24 have been used. The main objective is to check whether the periodicities observed in the cosmic rays are affected by the magnetic rigidity or the height at which the neutron monitors are placed. A Global Neutron Monitor (GNM) has been defined as representative of the neutron monitor global network. This GNM is constructed by averaging the counting rates of a set of selected neutron monitors. The selection process is based on the combination of three new data quality criteria, which are applied to neutron monitors in the Neutron Monitor Data Base giving a final pool of 22 selected neutron monitors. Morlet wavelet analysis is applied to the GNM and the selected solar activity parameters to find the common periodicities. Short-term periodicities of 13.5, 27, 48, 92, 132 and 298 days have been observed in cosmic ray intensity. A clear inverse relationship between rigidity and spectral power has been obtained for the 13.5-, 48-, 92-, 132-day periods. A not so clear but still observed direct relationship between the height of the neutron monitors and the spectral power for the 48-, 92-, 132-day periods has been also found. The periodicity of 92 days is the one which shows the highest dependence with rigidity cutoff and height. As far as we know, this is the first time that these dependencies are reported. We think that these observations could be explained by assuming some cosmic ray intensity energy dependence in such periodicities and a competitive effect between rigidity and height.


Introduction
The primary cosmic rays (CR) are electrically charged high energy particles, mostly originating in violent phenomena of our galaxy (such as supernova explosions, pulsars with very strong magnetic fields) and, to a lesser extent, solar or extragalactic phenomena, which continuously affect the terrestrial atmosphere with energies between 10 6 and 10 20 eV/nucleon. When these primary CR interact with particles present in the Earth's atmosphere, other particles called secondary cosmic rays, such as protons, muons, neutrons and mesons, are created. Some of these secondary particles are registered by ground-based neutron monitors (NMs) and muon telescopes, respectively.
The year 1964 was designated by cosmic ray investigators as the international quiet year of the sun (IQSY) and Hugh Carmichel developed a neutron monitor with a statistical accuracy of 0.1% for hourly data called "NM64" or also known as "supermonitor" (Shea and Smart, 2000). This "NM64" is used as standard to obtain data which to compare with other stations. From that date to the present day, other similar instruments have been developed along different locations on the terrestrial sphere. Most of them report their measurements through the Neutron Monitor Data Base (NMDB). To study the solar activity through NM requires the use of multiple stations located in different geomagnetic locations because their position with respect to the Earth's magnetic field determines the minimum energy of CR to generate counts on a given NM. This need justifies the foundation of the NMDB. The Cosmic Ray Intensity (CRI) is anticorrelated to the solar activity, measured by the sunspot number, with a certain delay caused by irregularities in the Interplanetary Magnetic Field (Forbush, 1958;Usoskin et al., 1998;Singh, Singh, and Badruddin, 2008). The historical maximum value of the CRI measured by the NM is being registered during this year 2020 and has not yet reached the maximum value, previously, the maximum value of the CRI was recorded at the end of Solar Cycle (SC) 23, during the solar minimum between 2007-2009. It was an unusual minimum followed by a weak solar maximum in the SC 24 (2012)(2013)(2014) compared to the last previous cycles (Ahluwalia, 2014).
The cosmic ray flux presents a quasi-periodic modulation of 27 days. This phenomenon is attributed to the solar rotation and is related to the coronal holes and co-rotating interaction regions, i.e. interactions between fast and slow solar streams (Grieder, 2001). Ahluwalia and Ygbuhay (2015) compare the modulation of cosmic rays along the SCs 20-24 and they conclude that this modulation is significantly weaker in the SC 24, although the variation of CRI with the magnetic rigidity is similar in all the studied SCs.
This period is not related to the solar electromagnetic radiation according to Poblet and Azpilicueta (2018). Nevertheless, temporal analysis indicates that the geometric configuration changes of the Sun's magnetic field and the Earth's magnetic field during a rotation could be an important factor in addition to any temporary solar wind structure in generating the 27-day period.
The CR suffer a modulation process when they pass through the heliosphere. In this process, CR flux is maximum coinciding with the minimum solar activity, i.e. low sunspot number and in the same way, minimum CR flux is given when the maximum solar activity occurs. This maximum in solar activity is reflected by a high number of sunspot and solar flares coupled with an inversion in the magnetic polarity of the Sun. These variations are repeated cyclically with a period of 11 years (Gnevyshev, 1967).
Finally, the well-known Rieger-type periodicity was firstly observed in the gamma-ray flares around the SC 21 by Rieger et al. (1985). This periodicity was also observed by Kudela et al. (2002) in the CRI. This period is not stable and only appears around solar maximum during the polarity change.
This paper examines the short-term periodicities, i.e. periodicities between 2-512 days, during the period 2013-2018 using the wavelet analysis. The facts that several different phenomena can contribute to different periodicities and that these are not constant in time could make Fourier analysis unsuitable (Schreiber, 1998) and therefore wavelet analysis has been used to observe the intrinsic temporal variation in the periods. This paper is organized by the following way: in Section 2 are described the data used and the wavelet analysis method; in Section 3 we present the obtained results; in Section 4 one finds the discussion and finally, in Section 5 the summary and the conclusions are described.

Data and Analysis Method
A set of neutron monitors with available data from 2013 to 2018 has been selected. These data are one hour pressure corrected counting rates and were collected from the NMDB web page. A check of data availability for the data provider stations in NMDB was performed initially. Only stations with less than 5% of missing values and/or outliers, i.e. clearly wrong data, were selected for the study. This gave 26 of 53 stations in NMDB and four of these 26 were taken out from the pool of selected stations because their wavelet spectrum showed an anomalous behavior when comparing to the rest of stations. This will be explained in detail in Section 2.1. Thus, the final selection of stations comprises 22 NMs at different vertical cutoff rigidity and height in meters above sea level (m.a.s.l.). The 22 stations and their characteristics are listed in Table 1.
Selected dataset covers the period from 1 January 2013 to 31 December 2018, i.e. the descending phase of the solar cycle 24. It was downloaded from http://www01.nmdb.eu/, the main website of the NMDB, and http://cr0.izmiran.ru/common/links.htm. Hourly averaged solar wind speed and interplanetary magnetic field data measured by ACE (Advanced Composition Explorer) from https://cdaweb.sci.gsfc.nasa.gov/index.html/ and daily observations of sunspot number (SSN) and the total solar irradiance (TSI) from https://www.ngdc.noaa. gov/stp/SOLAR/ and http://lasp.colorado.edu/home/sorce/data/tsi-data/ respectively, as complementary data and to establish a comparison of neutron monitor behavior and solar conditions.
Once the data from neutron monitor are selected, it is necessary to eliminate the outliers in the dataset before applying the wavelet analysis. To carry out the data cleaning is necessary to define the concept of outlier. An outlier or atypical value is an observation that is numerically distant from the rest of the data. These outliers can be removed or replaced by other values. An experimental point is consider as outlier if the measurement is out of the interval defined byx ± kσ x , σ x being the standard deviation of the data,x the mean value and k is an integer number whose value is 1, 2 or 3, depending on the study. This outlier detection criterion is only applicable to normal distributions, so that it must be checked if the data used in this study follow a normal distribution. Three normality tests counting rates have been applied to the neutron monitor: the Kolmogorov-Smirnov test, the Jarqua-Bera test and the Anderson-Darling test. The three ones provide the same conclusion, the neutron monitor counting rates do not follow a normal distribution during the period 2013-2018. With this result in mind, a different method to detect outliers into the dataset was chosen. The analysis of Box and Whiskers is a robust method to detect outliers which only depends on the data and not on the distribution of the data. The points out of the range [Q 1 − 1.5 (Q 3 − Q 1 ) , Q 3 + 1.5 (Q 3 − Q 1 )], Q 1 and Q 3 being the first and third quartile respectively, were considered outliers. The set of outliers and missing values were substituted by linear interpolation. The percentage of outliers with respect to the total set of points for a given neutron monitor is collected in Table 1 in the P.C.P. (Percentage of Corrected Points) column. The stations with P.C.P. > 5% were left out of the analysis as it has been explained above.
Once the outliers are substituted in the data from the selected stations, the wavelet analysis is applied to their counting rates to find the periodicities in the data. This analysis decomposes a signal into a sum of wavelets which come from a "mother" wavelet function, in this case the Morlet function. The continuous normalized Wavelet Transform applied to a discrete sequence x n with equal temporary spacing δt and n = 0, . . . , N − 1 is where x n are the different points of the signal and * is the conjugate complex of the "mother" wavelet function and the coefficients s and n are relative to the scale and the translation of the wavelet respectively. Concretely, we used the "mother" Morlet function with β = 1, given by: where ω 0 is the nondimensional frequency and η = s · t . The Spectral Power P is given by the square of the Wavelet Transform which is a complex number, P = |W n (s)| 2 . The Wavelet Power Spectrum (WPS) shows the temporal distribution of the Spectral Power for each frequency and on the other hand, the Global Wavelet Spectrum (GWS) is the average of the Spectral Power (P = |W n (s)| 2 ) at each frequency s: The significant peaks are above the 95% confidence level in the GWS and they are represented as a black line in the WPS figures in this paper. The black lines depict p-values less than 0.05 which means that the probability of the null hypothesis (there is no periodicity in the temporal series) is very low. The cone of influence in WPS figures is shown as a shaded area which corresponding to the zone influenced by the edge effects.
To determine significance levels above which |W n (s)| 2 represents real characteristics of the time series, a background spectrum is taken. Many geophysical and astrophysical time series can be modeled as either white noise or red noise. A simple model for red noise is the univariate lag-1 autoregressive process. The Fourier spectrum of this simple model is given by where k = 0, .., m/2 is the frequency index and .
x n is a time series in the discretized time n, N andx are the number of points and the average value of this time series, respectively. If a peak in the GWS is significantly above this background spectrum P k , then it can be assumed to be a true feature with a certain percent confidence, (Torrence and Compo, 1998). We determine the 95% confidence level (significant at 5%), to multiplying the background spectrum P k by the 95th percentile value for chi-squared distribution χ 2 function with two degrees of freedom. The calculated periodicities in this study show this confidence level and they are collected in Tables 3 and 4. The R package "WaveletComp" developed by Angi Röesch and Harald Schmidbauer has been used to perform wavelet analysis of time series (see the handbook by Roesch and Schmidbauer (2018)).

A Global Neutron Monitor
We define as a global neutron monitor (GNM) an ideal neutron monitor whose power spectrum is obtained using the averaged of the counting rates, in typified units, of all the selected neutron monitors in a particular study. Criteria based on three quality indices have been used to define the final list of stations to build the global neutron monitor.
These three quality indices are referenced as Q j , Q x and Q u ; they, by distinct methods, compare the different stations and give us a quantitative idea of which stations have a distant behavior of the general trend. For the calculation of the averaged counting rates we only took into account the NM whose value is greater than 0.5 in the three indices. This threshold comes up because, generally speaking, statistical studies consider that a value equal to or greater than ± 0.5 of the cross-correlation function is considered such as a moderate linear association positive or negative depending of the sign. Since the coefficients Q x and Q u encompass the cross-correlation function the same criterion is taken into account. This criterion is also applied to the coefficient Q j and it is considered that a value of 0.5 or higher indicates that the difference between two spectra is moderate or low respectively.

Quality Index Q j
We have defined a quality factor Q j to evaluate how much a single power spectrum differs from the one that is obtained for the reference station. The proposed quality index of each neutron monitor j is given by the following equation: where the second term to the right is the mean squared error: GWS j (s i ) are the points of the GWS obtained for a certain neutron monitor j , GWS(s i ) are the points of the GWS computed by averaging the GWS points of all the initial set of selected neutron monitors, and m = 801 is the number of discretized frequencies. Although this index is not normalized, when Q j → 1 implies that the station j behaves like the reference station. We consider the following criterion: if Q j < 0.5 the behavior of the station is very different from the average, if Q j ∈ [0.5, 0.65) the spectrum for this station is reasonable with respect to the average, for Q j ∈ [0.65, 0.85) the spectrum is good and finally if Q j ∈ [0.85, 1] is excellent. Four of the initial 26 neutron monitors, FSMT (Fort Smith), JUNG1 (Jungfraujoch-1), MWSN (Mawson) and TXBY (Tixie Bay), showed a Q j < 0.5. These stations were removed from the selected list and GWS(s i ) was recalculated taking into account only the final 22 selected NMs. These 22 NMs were used to construct the GNM spectra. Figure 1a shows the GWS obtained with the 22 neutron monitors with Q j > 0.5 (red line) and the GWS of two NMs with Q j < 0.5: JUNG1 (blue line) and TXBY (green line). The values of Q j for the 22 selected stations are collected in Table 2. Figure 1b shows the best and the worst station according to the Q j index computed with the initial set of selected NMs. Blue one represents NEWK (Q j =NEWK = 0.975), green line is YKTK (Q j =YKTK = 0.563) and red line is the average station.  Correlation between the measurements of the different NMs is expected since their variations are mainly due to the solar activity affecting Earth globally. Based on this assumption, a second index Q x is defined to evaluate the degree of correlation between stations as follows: is the normalized cross-correlation function between two time series x and y, and is the sample cross-correlation or cross-covariance, N being the number of points of the time series, x n+h the counting rate of a particular NM, y n the counting rate of a second NM to be compared with the first one,x andȳ are the average value of each time series, respectively, and h is the delay or temporary advance (h = 0, ±1, .., ±100 hours) between the two set of measurements. As a way of example, we explain below how Q x is computed for CALM station. First, the normalized cross-correlation function ρ xy between CALM and NEWK neutron monitors versus h is computed (Figure 2). The maximum value of ρ xy = 0.868 is found when h = −4 (red line). This h = −4 hours is compatible with the 0.21 h of delay between CALM and the Greenwich meridian and the 5.05 delay between NEWK and the Greenwich meridian taking into account the accuracy ±1 h of the computation. It is accepted that a cross-correlation coefficient ρ xy is statistical significant, at 95% confidence interval, if where N − h is the number of data points used in the calculation of the ρ xy where N = 52584, the number of data points of the time series, and h = 0, ±1, .., ±100 is the lag. In this case, ρ xy is statistically significant since ρ xy > 0.0087. Finally, the cross-correlation function between CALM and the rest of stations is computed keeping the maximum value of ρ xy (h) and the index Q x is calculated using equation 5. The 22 NMs gave Q x > 0.5 (fourth column in Table 2). Q x expresses the average degree of correlation between station x and the rest of stations so that Q x → 1 implies that the station x is strongly correlated with the remaining 21 NMs. We consider the following criterion: if Q x → −1 there is anticorrelation between station x and the remaining stations, if Q x → 0 the station x is not correlated with any other NM, if Q x < 0.5 the correlation is very weak or non-existent, if Q x ∈ [0.5, 0.65) the counting rates for this station x is reasonable correlated with respect to the other NMs, for Q x ∈ [0.65, 0.85) the correlation is strong and finally if Q x ∈ [0.85, 1], the station x and the other stations are very strongly correlated.

Quality Index Q u
The Q x index is based on the correlation between temporal series of neutron monitor counting rates, this Q u index is based on the correlation between the GWS of the NMs. The procedure to build Q u is similar to the used to compute Q x . Thus, Q u represents the degree of correlation between the GWS of a station and the others. Equation 5 is rewritten as follows: where u and v represents the GWS of each neutron monitor, M = 22 the number of NMs and ρ uv is the sample cross-correlation function between two GWS u and v, and m = 801 being the number of discretized frequencies, GWS u (s i+h ) is the GWS value for the frequency s i+h of a particular station u, GWS v (s i ) is the GWS value for the frequency s i of a second station to be compared with the first one, GWS u and GWS v are the average value of GWS, averaging for all study frequencies, of the monitor u and v respectively, and h is the delay (h = 0, ±1, .., ±100). Figure 3 shows the cross-correlation function ρ uv between CALM and SOPO for an interval of 1h temporal steps from −100 h to +100 h. The correlation maximum is on h = 1 (red line on the figure) what is expected for two NMs located close to the Greenwich meridian. The statistical significance is given by the equation The value of h that maximizes the statistical significance is 0.074, then for ρ uv > 0.074 the cross-correlation function is considered significant.
We consider the same criterion that with the last index: if Q u → −1 there is anticorrelation between the GWS of the station u and the GWS of the remaining stations, if Q u → 0  the spectrum of the NM u is not correlated with the GWS of the other NMs, if Q u < 0.5 the correlation is very weak or non-existent, if Q u ∈ [0.5, 0.65) the GWS for this station u is reasonably correlated with respect to the GWS of the other NMs, for Q u ∈ [0.65, 0.85) the correlation is strong and finally if Q u ∈ [0.85, 1], the GWS of station u and the GWS for another stations are very strongly correlated. All the stations have Q u > 0.5 (see Table 2, fifth column).

The Global Neutron Monitor Counting Rate
The counting rate at a given neutron monitor is a function of rigidity, altitude, the building where the monitor is located among other things. Nevertheless, we propose a virtual neutron monitor whose counting rate can be characterized by a counting rate that can be representative of the response on neutron monitor counting rates to solar activity. To build the global neutron monitor (GNM) counting rate, we calculate typified units or reduced centered units for the 22 neutron monitor in Tables 1 and 2. The positions on the Earth, latitude and longi-tude, of the selected NMs are depicted in Figure 4. If x n is the counting rate in a discretized time n in any neutron monitor, the typified unit in that time is given by wherex and σ x are the mean and standard deviation respectively. This reduced centered variable have null mean (z = 0) and standard deviation equal one (σ z = 1). Once z n is computed for a neutron monitor, the time is shifted according to the h obtained in equation 5. This shifts the counting rate to Greenwich meridian time. The sixth column in Table 2 collect the delays of each station with respect to CALM which has a similar geographical longitude to the Greenwich meridian. This procedure is repeated for the 22 NMs. Finally, the GNM counting rate in typified units is computed by averaging the 22 typified and shifted counting rates. This counting rate is shown in Figure 5a. The WPS and the GWS of the GNM are shown in Figure 5b and 5c, respectively. The results are similar to those obtained for the 22 NMs separately. The GNM counting rate will be used in Section 3 to calculate the CR periodicities and the degree of correlation of CR with different solar parameters.

Neutron Monitors
Wavelet analysis applied to the neutron monitor counting rates presents, in general, similar results for all the NMs. The WPS and GWS computed for the 22 selected NMs are gathered in the appendix. The 1-day period, associated with the terrestrial rotation, was observed in all stations but it has not been included because it is out of scope of this study. The GNM shows a representative behavior of all stations ( Figure 5) and so, it can be used to describe the characteristics of the periodicities observed in cosmic rays counting rates. The detected periodicities are listed in Table 3.
The most prominent period observed in the GNM is the ≈ 27-day one and it appears throughout the complete time series except from mid-2013 to mid-2014 ( Figure 5). This period, related to the synodic solar rotation, was observed by the 22 stations (see the figures in the appendix).
A second peak, probably related with solar rotation, is centered in ≈ 13.5 days. The 13.5-day period (the second harmonic of the fundamental solar rotation) has been associated with both solar active longitudes as well as tilted dipole structure (Vipindas, Gopinath, and Girish, 2016). The presence of the heliospheric current sheet (HCS) splits the heliosphere in two and the Earth orbit cross the HCS twice every solar rotation when the tilt angle is large enough. This higher tilt angle is expected out of the solar minimum. This period presents higher power between 2014 and 2017 ( Figure 5b) and it is especially strong in the second half of 2016 during the decreasing phase of the solar cycle. This periodicity is much weaker and almost undetectable from 2017 to the end of the studied time interval in coincidence with the solar minimum.
The nearly annual period is the second highest peak (   Table 3 Periodicities in GNM counting rates, sunspot number (SSN ), total solar irradiance (T SI ), interplanetary magnetic field B and solar wind speed V and its components. The second line expresses the notation used to refer to each obtained periodicity in the different ranges.
The three periods were identified between 2014-2016 with similar spectral power in all the NMs and in coincidence with a lower CRI. El-Borie, Aly, and El-Taher (2011) performed a study based in measurements of three NMs (Kiel, Rome and Huancayo/Haleakala) during three solar minima (1975-1976, 1985-1986 and 1994-1995) and three solar maxima (1979-1980, 1989-1990 and 2000-2001). He reported a set of periodicities coincident with the ones that we have observed. The 27day and 13.5-day period were observed during the three solar maxima and the three solar minima obtained in all cases. A 9-day period was detected along to the three solar minima and the 1979-1988 solar maximum. This period is not clearly detected in this work. This can be due to time interval covered. From the solar maximum to the solar minimum and in that stretch the 9-day peak is weaker.
The 48-day period was observed by the three NMs during the 1975-1976 solar minimum. This could be consistent with our observations as it is only observed in the descending phase of the solar cycle ( Figure 5).
The 92-day period was registered in Kiel and Rome during the 1989-1990 solar maximum; nevertheless, Huancayo/Haleakala did not register this period possibly due to its high cutoff rigidity (R c = 13.3GV ) implying that the period could be related to variations in the cosmic ray flux at energies below 13 GeV/nucleon. The period between 2013 and 2018 analyzed in this paper comprises part of the solar maximum, the decreasing phase of the cycle and, probably, the beginning of the solar minimum. These two periodicities were clearly observed at the beginning of the descending phase what could explain the El-Borie, Aly, and El-Taher (2011) observations.
We must to point out the 154 day period known as Rieger period (Rieger et al., 1985) has not been found in our analysis. Rieger-type periodicity usually occurs near maxima sunspot cycles and it was firstly reported in cosmic ray by Kudela et al. (2002). In a recent paper, Saad Farid (2019) reported the observation of a Rieger-type periodicity in the cosmic ray flux in 2012, i.e. the beginning of the solar maximum. 2012 is not included in the present work but 2013 it is, and the Rieger-type periodicity was not observed at the end of the solar maximum before the solar polarity inversion.
One of the open questions regarding the observation of periodicities by neutron monitors is the role of rigidity and altitude above sea level on the relative strength of a particular periodicity. Figure 6 shows the spectral power of the 13.5-, 27-, 48-, 92-, 132-and 298-day peaks vs. magnetic rigidity R c . A different behavior between southern NMs (green circles in Figure 6 The spectral power of the 13.5-, 27-, 48-, 92-, 132-and 298-day period vs. magnetic cutoff rigidity R c of each neutron monitor. In green, the neutron monitor located in the southern hemisphere and in blue the NMs located in the northern hemisphere. The black line is the regression line obtained by linear adjust for the northern hemisphere NMs. R 2 is the coefficient of determination obtained for the northern hemisphere stations. The grey band displays the error range of fit line. Figure 6) and northern NMs (blue circles in Figure 6) has been detected. Southern stations do not show dependence but this is not truth for northern stations. A clear dependence between the spectral power and rigidity is not observed for the 27-and 298-day periods. Nevertheless, a decreasing relationship is observed for 13.5-, 48-, 92-and 132-day periods, being stronger for the 92-day one (R 2 = 0.5422). As far as we know, this is the first time that this kind of dependence has been reported. The result obtained for the 92-day period could also explain the observations by El-Borie, Aly, and El-Taher (2011) commented on above, although in that case the altitude of the station (Haleakala 3030 m.a.s.l. and Huancayo 3400 m.a.s.l.) could play a competitive role against the station rigidity as is explained below. The dependence of the spectral power with the altitude of the NM site has been also studied. A direct result appears when spectral power is represented against altitude (Figure 7). NMs at sea level (below 200 m.a.s.l.) do not show a relationship with height (green circles in Figure 7). Nevertheless a direct relationship emerges for stations above 200 m.a.s.l. (blue circles in Figure 7) for the 48-, 92-and 132-day periods. This seems to suggest rigidity plays the main role in the spectral power for low altitude stations, but this role change as the station height enhances. For stations at sea level or low altitude (less than 200 m), the effect of cutoff rigidity is the dominant factor. This fact would explain the strong dispersion of the spectral power values for these stations (neutron monitors with altitude < 200 m). From approximately 200 m the altitude begins to play a significant role. The highest correlation is observed for the 92-day period with a R 2 = 0.3621. The same, the highest coefficient of determination, is observed in the relationship between spectral power and rigidity in the 92day period. The spectral power decrease with rigidity and the spectral power enhancement with altitude could be both explained if the cause forcing this period is more effective on CR with less energy. The contribution of CR with lower energy to the neutron monitor counting rate is lower for higher rigidity and higher for higher altitudes. This implies stronger periodic variations in the CRI and so that higher spectra power. Regarding the dependence with altitude, it has to be taken into account that rigidity and height have opposite competing effects on the weight of low energy cosmic ray on the neutron monitor counting rates. The NMs considered in this study are located at heights between 260 and 3570 m.a.s.l. (blue circles in Figure 7) and in a range of rigidities between 0 and 8.53 GV. According to our results it could be said that for this selection on the NMs the height effect is dominating the rigidity effect. Nevertheless, out of this selection this balance can change. This change could explain because El-Borie, Aly, and El-Taher (2011) did not detect the 92-day period in Huancayo/Haleakala data. The neutron monitors were at a height of about 3000 m.a.s.l. but at a rigidity of 12.9 GV.

Solar Parameters
The complete study of periodicities in the cosmic ray flux have to include the context that could be driving the observed periodicities. This context is what is analyzed under the epigraph solar parameters. We have selected four magnitudes that have a clear influence on the cosmic ray propagation throughout the heliosphere. These four magnitudes can be split into two categories. Those related to the solar wind conditions and those related to the solar cycle status. Regarding solar wind conditions, solar wind velocity and interplanetary magnetic field play a central role on the cosmic ray propagation. Many papers have point out this role, for instance, Blanco et al. (2013) analyzing non-periodic events such as Forbush decreases, concluded the central role of the solar wind speed, magnetic field strength and the deceleration of interplanetary coronal mass ejections on the depth of the Forbush decreases. Regarding the solar cycle status, the sun spot number (SSN) and the total solar irradiance (TSI) are the chosen magnitudes. SSN is a clear proxy of the solar activity. Moreover, it is well known the anticorrelation between the SSN and the cosmic ray intensity. A recent paper regarding this relationship is Ross and Chaplin (2019). Similar comments, regarding solar activity and solar cycle status can be done about the TSI. Kopp (2016) revised the observed temporal variations of the TSI.

Solar Wind Velocity
The solar wind is the medium through which cosmic rays propagate in the heliosphere and the solar wind speed defines the solar wind in slow solar wind and fast solar wind. The solar wind regime, i.e. slow or fast, establishes the transport conditions for the cosmic rays. It is therefore to be expected that the periodicities observed in the speed of the solar wind might appear in the intensity of the cosmic rays.
As it is already commented in Section 2, solar wind velocity data were recorded by ACE spacecraft at L-1 position. The Geocentric Solar Ecliptic (GSE) coordinates has been used and the speed V = V 2 x + V 2 y + V 2 z , the V x component and the perpendicular component of solar wind velocity, V perp = V 2 y + V 2 z , have been analyzed accordingly with the methodology described Section 2. The result of this analysis is shown in Figure 8.
The WPS and GWS are similar in V and V x what is expected because the solar wind is strongly radial. The GWS shows clearly marked periodicities in 27, 13.5 and 7 days. These periodicities can be related to solar rotation and the presence of the HCS which is usually embedded in the slow solar wind. Other two weaker periodicities, 92 and 298 days are observed (bottom panel Figure 8). The detected periodicities are presented in Tables 3 and 4. The evolution in time of these periodicities is revised with the WPS. 27-and 13.5-day periodicities are present along the complete period. Otherwise, 298-day periodicity is also present along the whole period although its power is weaker than the 27-and 13.5-day ones. Something different happens with the 92-day period. It is stronger near the solar maximum and weakens as it approaches the solar minimum.
The behavior of V perp is completely different from those observed for V and V x and, as we explain below, it is strongly influenced by the spacecraft orbital parameters. It shows a prominent peak at 176 days (GWS in Figure 8), but it is not due to solar phenomena. This outstanding periodicity is related with the orbital period of ACE of approximately 177 days and appears throughout the entire interval in the WPS (Figure 8). Other detected periodicities in V perp are 92, 59 and 27 days.

Interplanetary Magnetic Field
Because of cosmic rays are charged particles, the interplanetary magnetic field governs their transport through the heliosphere. This is clear attending the dependence of the shape of the cosmic ray intensity curve with the polarity of the solar cycle (Potgieter, 2008). The effect that the field has on the propagation of cosmic rays depends on their energy which determines the length of their gyroradious. The value of the gyroradius determines the spatial size that a given magnetic structure must have in order to effectively affect the transport of the cosmic rays. With these thoughts in mind, it is reasonable to expect that periodic variations in the magnetic field can be induce periodic variation in cosmic ray.
As has said in Section 2, interplanetary magnetic field data from ACE spacecraft magnetograph has been used. The same procedure describe in that section has been applied to magnetic field measurements to search the main periodicities of the data in the magnetic intensity B = B 2 x + B 2 y + B 2 z , the B z component and the radial component in the ecliptic plane B r = B 2 x + B 2 y . The GWS shows a different behavior for each of the analyzed magnetic magnitudes (Figure 9). The radial component shows three clear periodicities of 9, 13.5 and 27 days. An additional weaker periodicity of 383 days is also observed. On the opposite, B z component presents no clear signatures of periodicities, nonetheless the 27-day and 59-day periods were still detected although with a not very strong spectral power. Finally, the magnetic field intensity presented a mixed behavior with periodicities of 7, 13.5, 27, 48, 92 and 298 days. These periodicities are included in Tables 3 and 4. Regarding the persistence of a periodicity along the studied period, the WPS shows that the periodicities of 7, 13.5, 27 days detected in the magnetic field intensity are very irregu- lar, disappearing several times briefly. The same is observed for the B z component. Something completely different is observed in B r , the 13.5-and the 27-day periodicities remain throughout the entire period. We have to point out that the 27-day periodicity turned weaker during the first months of 2014, the central months of 2017 and the end of 2018. During these periods, the 13.5 periodicity turned stronger.

SunSpot Number and Total Solar Irradiance
The SSN is a widely used proxy for the solar activity. Sunspots are grouped around active regions on the Sun surface and their number grows up when the solar cycle moves toward the solar maximum. Their spatial distribution on the Sun surface is not homogeneous and their number and location on the solar disc changes with the solar cycle. Sunspot lifetime is ranged from hours to weeks, so their spatial distribution can change from one solar rotation to the next one.
The TSI is determined by the balance between the presence of sunspots and faculae on the solar disc. On solar-cycle and solar-rotation timescales, the majority of the TSI-amplitude fluctuations can be explain by this balance and sunspots and faculae are manifestations of the magnetic field on the solar surface (Lean, 2010).
The GWS of the SSN presents a strong periodicity of 27 days. Other detected but less strong periodicities are 48, 92 and 298 days (Figure 10).
In previous papers, Joshi, Pant, and Manoharan (2009)  The persistence of these periodicities from 2013 to 2018 is shown by the WPS (first panel in Figure 10). It is clear that all of them persist along the complete period except for 2018 when the sunspots practically disappeared from the solar disc.
Although the TSI has a dependence on the SSN, the main observed periodicities are different from those found for the SNN. These periodicities are 9, 13.5, 27, 48, 92, 190 and 298 days (bottom panel in Figure 11 and Tables 3 and 4). It must be remembered that the TSI emerges from the competition between sunspots and faculae. As in the case of SSN, the TSI periodicities are present from 2013 to 2017, disappearing in 2018 The WPS of TSI shows a blue spot in 2014 that is due to missing data that has been corrected for by linear interpolation (third panel in Figure 10).

Discussion
The variations in the flux of CR observed by ground level instruments such as neutron monitors correspond to changes in the propagation conditions of CR in the heliosphere. Recurrent variations in CR flux should reflect recurrent variation in the CR transport conditions. Transport models describe the modulation of CR in terms of the solar wind speed, diffusion tensor and the heliospheric magnetic field. For a revision as regards modulation of CR, see Potgieter (2013). The chosen parameters selected in this paper, heliospheric magnetic field, solar wind velocity, sunspot number and total solar irradiance should have an effect on the propagation conditions of CR in the heliosphere. Similar periodicities in such magnitudes could explain the observed periodicities in cosmic ray counting rates. On the other hand, the use of NMs at different cutoff rigidities allows one to evaluate the behavior at different cosmic ray energy thresholds. Tables 3 and 4 show the periodicities found in this study. The error intervals of the periodicities have been calculated fitting each peak to a Gaussian function and calculating its standard deviation. Six periodicities with associated periods of 13.5, 27, 48, 92, 132 and 298 days, have been found in GNM counting rates during the interval 2013-2018.
After evaluating the relationship between the power of the spectral power for the observed periodicities with rigidity and the height at which the neutron monitor is located, two main results arise. A rigidity dependence has been observed in the power of the spectral power in some of the periodicities (Figure 6). These were the 13.5-, 48-, 92-and 132day periods. This dependence decreases with the neutron monitor vertical cutoff rigidity showing coefficients of determination of R 2 13.5d = 0.362, R 2 48d = 0.287, R 2 92d = 0.542 and R 2 132d = 0.133, respectively. This dependence that could not be significant in the 132-day period, is not observed in the 27-and 298-day periods. These three periodicities, 13.5, 48 and 92-day periods, show higher power when the cutoff rigidity is lower, i.e. the period is clearer in NMs with lower energy threshold for primary CR. This result could be explained if the transport conditions are changing with these periods in a way that CR with less energy are more effective affected than CR with higher energies.
The second main observation inferred from our work is that there also exists a dependence with respect the height at which neutron monitor is located. This dependence is clear at heights above 200 m.a.s.l. and, as the rigidity dependence, only is observed for certain periodicities. These were the 48, 92 and 132-day periodicities. This dependence increases with the altitude showing coefficients of determination of R 2 48d = 0.112, R 2 92d = 0.362 and R 2 132d = 0.273, respectively (see Figure 7). The study of variation of CRI with altitude is known as altitude effect. The secondary CR cascade starting at the upper atmosphere gradually gets intensified at lower altitudes and have a maximum radiation at a height about 15 -20 km depending on the latitude (Bazilevskaya, Krainev, and Makhmutov, 2000). In other words, the CRI increases with the altitude and reaches a maximum at a height of about 20 km, known as the Regener-Pfotzer maximum. Below this height there is a gradually fall in intensity due to absorption and decay processes (Grieder, 2001). The experimental results are similar at different places of the earth. The higher altitudes, the higher intensity of secondary CR.
The dependence with rigidity and height shows an opposite behavior regarding the power of the spectral power on GNM counting rate and both effects are observed for the 48-, 92and 132-day periods. These effects are clearer for the 92-day periodicity. As it has been already commented in Section 3.1, we think that this behavior could be explained considering a cosmic ray energy dependence. The contribution of CR with lower energy to the neutron monitor counting rate is lower for higher rigidity and, at a same rigidity, higher for higher altitudes. This could imply stronger periodic variations in the CRI and so that higher spectra power. Although the entrance of cosmic ray to Earth atmosphere is ruled by the rigidity of cosmic rays, rigidity and height have opposite competing effects on the weight of low energy cosmic ray on the neutron monitor counting rates. The range of NM heights is between 260 and 3570 m.a.s.l. (blue circles in Figure 7) and the range of rigidities is between 0 and 8.53 GV. According to our results it could be concluded that for our selection on the NMs the height effect is dominating the rigidity effect. Nevertheless, out of this selection this balance could change. This change could explain because El-Borie, Aly, and El-Taher (2011) did not detect the 92-day period in Huancayo/Haleakala data. The neutron monitors were at a height of about 3000 m.a.s.l. but at a rigidity of 12.9 GV.
The periodicities detected in CRI are also observed in the selected solar parameters (Table 3). We discuss now the periods in solar parameter and their role in forcing CRI periodicities.
The 48-day period is also observed in the module of interplanetary magnetic field, SSN and TSI while the 92-day period is registered in the same parameters in addition with the solar wind velocity. On the other hand, previous work has reported the relationship between solar wind speed and the fractional coronal hole area (Tulasi Ram, Liu, and Su, 2010) and between the TSI and the presence of coronal holes (Woods, 2010). Maybe, this is pointing to a possible relationship between solar wind structures and the observed flux of cosmic ray which is reflected by the observed 48 and 92-day periodicities. Close to the 92-day period observed in cosmic rays, a periodicity of 88 days has been found in the solar wind velocity. This period is only observed in the solar wind velocity and could be a harmonic of the ACE orbital period of 177 days. Nevertheless after analyzing the ACE orbit, no such period of 88 days is detected, therefore it should be concluded that 88 days is a true period in the solar wind velocity. This could reinforce the hypothesis of solar wind structures producing energy depending periodic variations in the cosmic ray flux.
The 13.5-and 27-day periods, related with the solar synodic rotation, are observed in almost all the analyzed parameters. The most prominent peak in the power of the wavelet analysis for CR counting rate is in 27 days ( Figure 5). This period is attributed to the solar rotation commonly, and is related to the coronal holes and co-rotating interaction regions, i.e. interactions between fast and slow solar streams. On the other hand, changes of the geometric configuration in the Sun's magnetic field and the Earth's magnetic field along a solar rotation could be also, an important factor on the 27-day period in cosmic ray flux (Grieder, 2001;Poblet and Azpilicueta, 2018).
The 132-day period was only detected in CR. This could imply that this periodicity is induced by a process not related with solar wind velocity, heliospheric magnetic field, SSN or TSI. However, the relationship of this periodicity with a solar phenomenon cannot be ruled out. It has been reported the observation of the period of 134.5 days in the sunspot data of the northern hemisphere and a similar period (133 day Finally, the 298-day period is detected in the magnetic field intensity, TSI, SSN and in the component V perp perpendicular to the ecliptic plane in addition to GNM counting rates. The 132-and 298-day period are very close to subharmonics of the solar synodic rotation of ≈ 27 days, the fifth (≈ 135 days) and eleventh (≈ 298 days) subharmonic. Therefore, the 132-and 298-day period could related to the fifth and eleventh subharmonic of the solar synodic rotation of ≈ 27 days.
Having the same period by two magnitudes does not mean a causal relationship or a common origin producing such period in both magnitudes. Cross correlation function has been computed to check if there would be common origins and/or causal effects among cosmic rays and solar magnitudes. Cross correlation among cosmic rays and different solar parameters has been studied by several authors in previous work. Chowdhury and Kudela (2018) analyzing 27-day averaged data from solar cycles 21 to 24 obtained a normalized crosscorrelation coefficient of ≈ 0.74 − 0.82 and ≈ 0.97 between CRI and SSN and CRI and interplanetary magnetic field respectively. Singh and Mishra (2018) determined the values −0.87, −0.75 and −0.60 for the cross-correlation coefficients between CRI and SSN during SCs 22, 23 and 24, respectively, i.e. a decrease in the degree of correlation is observed as solar cycles pass. Recently Oloketuyi et al. (2020) obtained correlation coefficients of −0.72 and −0.73 between CRI and SSN in the SCs 23 and 24, respectively, which is consistent with our work.
Cross-correlation function among GNM counting rates and solar parameters has been obtained. We calculate the cross-correlation function with different lag times being between −100 and +100 hours for solar wind velocity and interplanetary magnetic field and between −10 and +10 days for SSN and TSI. The maximum value of this function for each case is  Table 5. From this Table 5 is clear the anti-correlation between SSN, TSI and magnetic field with a time lag of +9 days, −2 days and −54 hours, respectively. However, GNM counting rate and solar wind velocity do not show a correlation. The value obtained in this work for the SSN is consistent with previous studies but B reflects a slight anticorrelation while Chowdhury and Kudela (2018) obtains a very high anti-correlation.

Conclusions
A new representative measurement of a set of NMs has been defined by imposing three quality indices for the selection of the NMs, we named this the Global Neutron Monitor (GNM). It has been defined for a set of 22 stations listed in Table 1. This GNM has been calculated by averaging the counting rates in typified units of these 22 stations. The Morlet Wavelet Analysis is applied to the GNM counting rate and to the set of 22 NMs in the period 2013-2018, corresponding to the end of solar maximum and the decreasing phase of the Solar Cycle 24. Before applying this spectral analysis, an outlier detection is carried out on the 22 NMs. To detect outliers, we performed the analysis of Box and Whiskers, a robust method which only depends on the data and not on the distribution of the data. The set of outliers and missing values were corrected by linear interpolation. The 13.5-, 27-, 48-, 92-, 132 and 298-day periodicities have been found. The higher spectral power corresponds with 27-day period. An inverse dependence between the spectral power and rigidity has been found for the 13.5-, 48-, 92-and 132-day periods. A second direct dependence between the spectral power and the height of the neutron monitor for the 48-, 92-and 132-day periods. The 92-day period showed the highest correlation in both cases. To the best of our knowledge, this is the first time that these dependencies are reported. We think that these observations could be explained by assuming some CRI energy dependence in such periodicities and a competitive effect between rigidity and height.
The same analysis has been applied to sunspot number, total solar irradiance, interplanetary magnetic field and solar wind velocity in order to find common possible origins/relationships with the CRI periodicities. The 13.5-, 27-, 48-, 92-and 298-day periodicities have been also found in different solar parameters. Nevertheless, the 132-day period has not been found in solar parameters. Additional periodicities have been found in solar parameters (Table 4).
The cross correlation has been performed to evaluate dependencies between CRI and solar parameters along the complete studied period. A strong anticorrelation (< −0.7) is found with SSN and TSI. A slight anticorrelation with the interplanetary magnetic field. No correlation is observed between CRI and solar wind velocity. The whole studied period (2013)(2014)(2015)(2016)(2017)(2018) comprises the end of the last solar maximum, the descending phase and the solar minimum. Maybe, this period including different phases in the solar cycle evolution could explain the poor cross-correlation between CRI and magnetic field or solar wind velocity.

Figure 17
The WPS and GWS of Thule and Yakutsk NM.