Scaling Features of Diurnal Variation of Galactic Cosmic Rays

We analyze the scaling properties of the diurnal variation of galactic cosmic rays (GCRs) in Solar Cycle 24 and the solar minima between Solar Cycles 23/24 and 24/25 for 2007 – 2019 based on the count rates of the Oulu, Newark, Hermanus, and Potchefstroom neutron monitors. The scaling features of the GCR diurnal variation are studied by evaluating the Hurst exponent, a quantitative parameter used as an indicator of the state of the randomness of a time series. We estimate the Hurst exponents for GCR diurnal-variation parameters amplitude and phase using structure-function and detrended-fluctuation-analysis methods. Results show that the Hurst exponents for the GCR diurnal variation vary in the range from ≈0.3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\approx0.3$\end{document} to ≈0.9\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\approx0.9$\end{document}, with a general tendency of being systematically above 0.5. It suggests that the GCR diurnal variation reveals a more persistent structure than Brownian motion. However, the time series of GCR diurnal-variation amplitude and phase evolve from a more persistent structure in the solar minimum between Solar Cycles 23/24 in 2007 – 2009 to a more random character in and near the solar maximum 2012 – 2014. This observation seems to be in agreement with the general configuration of the heliosphere through the 11-year solar-activity cycle. Moreover, the temporal profile of the Hurst exponent for GCR diurnal amplitude and phase around the beginning of the solar minimum between Solar Cycles 24/25 (2018 – 2019) differs from the solar minimum between Solar Cycles 23/24 in 2007 – 2009, suggesting a dependence on solar-magnetic polarity. These findings could shed more light on GCR particle transport in the turbulent heliosphere over the solar cycle.

The average solar diurnal GCR variation is shaped due to the modulation of the GCR flux by the solar wind. This has been explained, based on the diffusion-convection theory of the GCR heliospheric propagation (Ahluwalia and Dessler, 1962;Krymsky, 1964;Parker, 1964), as a consequence of an equilibrium established between the radial convection of GCR particles by the solar wind and the radial component of the inward diffusion of the GCR particles, owing to the existence of a radial gradient, along the heliospheric magnetic field. This picture has later been developed to include the effects of particle drifts in the large-scale magnetic field (Jokipii, Levy, and Hubbard, 1977;Kóta and Jokipii, 1983). Anisotropy for the part of the time interval studied here, i.e. 2001 -2014, was discussed by Tezari and Mavromichalaki (2016) for two neutron monitors having various cut-off rigidities located at the same longitude and different latitudes: Athens and Oulu. They have shown that the GCR diurnal anisotropy behaves in a different way during the opposite phases of the solar cycle, depending on the solar magnetic-field polarity, confirming previous results (e.g. Munakata et al., 2014). Thomas et al. (2017) underlined that the differences between solar cycles turned out to be greater when the solar transients visible in GCR were removed. The dependence of anisotropy on the heliospheric magnetic polarity was analyzed by Alania, Bochorishvili, and Iskra (2003), who have shown that within positive sectors anisotropy is greater during the 1976 solar minimum. Kudo and Mori (1990) established that the time intervals with greater values of anisotropy amplitudes from muon telescopes in the declining phase of Solar Cycle (SC) 21 appeared two years before that observed by neutron monitors. Tiwari, Singh, and Agrawal (2012) investigating Solar Cycles 20 -23 showed clear dependencies on the level of the solar-activity cycle, revealing that annual-average anisotropy reaches maximal values during the declining phases for all solar cycles, while during the minimum epochs, anisotropy took the lowest values. Moreover, they presented evidence of substantial differences between the descending epochs of the odd (21 and 23) and even (20 and 22) solar cycles. The GCR anisotropy amplitudes for the even SCs are smaller than for the descending epochs of the odd SCs, while for the ascending epochs there are no significant dissimilarities. Studies of the GCR diurnal anisotropy prove its utility in the recognition of other aspects of heliospheric changes, for instance analyzes of Forbush-decrease features (Gololobov et al., 2017), or a recurrence of the GCR anisotropy as a unique proxy of the level of solar activity, as well as solar-wind states (Mavromichalaki, Papageorgiou, and Gerontidou, 2016;Modzelewska and Alania, 2018).
In this article, we continue research on the temporary changes of GCR diurnal variation alongside the solar-cycle variations. Previously, we showed the usage of diurnal variation for determining the modulation parameters characterizing the heliospheric transport of cosmic rays using NM data (Ahluwalia et al., 2015), muon-telescope observations (Ahluwalia and Modzelewska, 2020), and Forbush ion-chamber data (Ahluwalia and Modzelewska, 2021), earlier discussed also by, e.g., Bieber and Chen (1991), Bieber (1993), andMunakata et al. (2014).
We studied (Modzelewska et al., 2019) in detail the behavior of the GCR diurnal variation over Solar Cycle 24 and the solar minima between Solar Cycles 23/24 and 24/25 (the interval 2007 -2018) and we estimated the contribution of drift effect to GCR modulation for the solar minima between Solar Cycles 23/24 and 24/25. We examined the polarity and rigidity dependence of the recurrent GCR intensity and anisotropy by Advanced Composition Explorer (ACE)/Cosmic Ray Isotope Spectrometer (CRIS), Solar Terrestrial Relations  (Modzelewska and Gil, 2021). The above-mentioned results concern the different average regular patterns of the GCR diurnal variation on different time scales, e.g. the annual or three-year averaging of the GCR diurnal variation used for the calculation of heliospheric-modulation parameters or analysis of its periodic character connected with solar cycle, solar magnetic cycle, and solar rotation. However, the GCR diurnal variation besides the regular character (on time scales of solar rotation, one year, or longer, etc.) has also a highly fluctuating term at short temporal scales (from a few to dozens of days). Moreover, recent studies carried out by Christodoulakis et al. (2019) and devoted to GCR intensity observed by selected NMs underlined the importance of considering smaller scales also. Therefore, the aim of this article is to analyze the scaling properties of the fluctuations of the GCR diurnal variation, focusing on relatively short temporal scales up to ≈ 100 days compared to the ≈ 11-year solar-activity cycle.
In particular, we compute the quantitative parameter called the Hurst exponent, which allows us to consider the state of randomness of the GCR diurnal-variation fluctuations through the SC 24 and the solar minima between Solar Cycles 23/24 and 24/25. In the analysis, data from four NM stations are used, and two independent methods for the determination of Hurst exponent are applied: structure-function scaling and detrended-fluctuation analysis. It is worth stressing that the Hurst exponent has been widely used in various aspects of space-weather studies (e.g. Takalo and Timonen, 1998;Wanliss, 2004;Panchev and Tsekov, 2007;De Michelis et al., 2021;Alberti, Consolini, and De Michelis, 2021), solar activity (Ruzmaikin, Feynman, and Robinson, 1994), or GCR intensity (Christodoulakis et al., 2019), but this is the first time that this parameter is applied for the description for the GCR diurnal variation.
This article is organized as follows: Section 1 is an introduction. In Section 2 the data used in this investigation are briefly presented. In Section 3 we describe the structurefunction and detrended-fluctuation-analysis methods. Section 4 shows the main results of the article. The last section discusses the implications of our findings and summarizes the results.

Data
In our calculations we use the GCR count rates with one-hour time resolution from four neutron monitors with various cut-off rigidities, from ≈ 0.8 GV up to almost 7 GV, located in the Southern and Northern Hemispheres. Details of the stations, such as their location or cut-off rigidities, are gathered in Table 1.
We analyze the daily amplitude and phase of the diurnal (24-hour) variation of the GCR intensity, calculated by means of normalized and detrended (excluding 25-hour trend) hourly data of the GCR intensity using the harmonic-analysis method (e.g. Gubbins, 2004). In our computations of the GCR diurnal variation in 2007 -2019, we follow the procedure described by Modzelewska and Alania (2018) and Modzelewska et al. (2019). We exclude from consideration the diurnal amplitudes > 0.7% as irregular events related to transient disturbances in the heliosphere, generally connected with coronal mass ejections (CMEs) and solar flares. Various methods can be used for padding gaps in the observations (e.g. Munteanu et al., 2016). Here we fill data gaps with random values from empirical distributions of GCR diurnal-variation amplitude and phase.
The daily values of the GCR diurnal-variation amplitude and phase for Newark NM (as an example) in SC 24 and the solar minima between Solar Cycles 23/24 and 24/25 (for 2007 -2019) are presented in the second and third panel of Figure 1. Moreover, as an indicator of the 11-year solar-activity cycle, the upper panels of Figure 1 presents the daily values of the sunspot number (SSN). Figure 1 shows the temporal evolution of parameters of diurnal variation of cosmic rays for the Newark NM at the resolution of one day for SC 24 and solar minima 23/24 and 24/25 in 2007 -2019. Although the dispersion of the data is very high, one could recognize the regular behavior of the parameters of GCR diurnal variation through the solar cycle with an average amplitude ≈ 0.3%. Namely, there is observed a tendency for the 11-year variability of GCR diurnal amplitude to be in phase with the solar-activity cycle (Figure 1 (2010) established that the higher energy NMs contribute more to the 11-year variation of GCR anisotropy due to the diffusion process.
Since the nature of the GCR diurnal variation is highly scattered with relatively small amplitudes (< 0.7%) compared to the solar-cycle effect of modulation of GCR intensity being, e.g. ≈ 10% for Newark NM, we study here the scaling properties of the time series assessing the Hurst exponent, which can be used as an indicator of the randomness of the GCR diurnal-variation fluctuations.

Methodology
We study here the scaling properties of the time series assessing the Hurst exponent [H ]: the parameter used as an indicator of the randomness state of a time series (e.g. Fuss, Weizman, and Tan, 2021). The Hurst exponent, a classical self-similarity parameter, measures three types of behavior in a time series: persistence, randomness, and anti-persistence. More precisely, when the process has a long-range dependency (i.e. correlated) structure, then the Hurst exponent is in the interval 0.5 -1. This situation is typical for processes observed to have more slowly evolving variations (i.e. more persistent), e.g. in monofractal and multifractal time series (e.g. Ihlen, 2012). When the time series has an independent or uncorrelated structure, such as in Brownian motion, then the Hurst exponent is close to 0.5. For time series with anticorrelated structure (anti-persistent), the Hurst exponent is in the interval 0 -0.5.
The concept of the Hurst exponent has its origins in hydrology (Hurst, 1951), and from that time it has been increasingly used in other disciplines: in finance (e.g. Di Matteo, 2007), in biomedical time series (e.g. Ihlen, 2012), or space-weather studies (e.g. Takalo and Timonen, 1998;Wanliss, 2004;De Michelis et al., 2016Alberti, Consolini, and De Michelis, 2021), as well as in solar-activity predictions (e.g. Singh and Bhargawa, 2017). Moreover, using the Hurst exponent (Ruzmaikin, Feynman, and Robinson, 1994;Oliver and Ballester, 1996) revealed the stochastic character of the solar-activity temporal profile. The Hurst-exponent properties were used in the analysis of compound-diffusion properties, i.e. when the particles are closely connected to the magnetic-field lines and the perpendicular transport origins in the random walk (Kóta and Jokipii, 2000). Recently, Christodoulakis et al. (2019) investigated the scaling features of the GCR intensity for 2000 -2017 using NM data. They found that H -values were ≈ 1 for the whole period and consequently GCR intensity time series exhibit positive long-range correlations. Moreover, multifractal behavior was detected at all NMs with a great degree of multifractality only for smaller time scales (≤ 100 days).
There are several techniques for the determination of the H -exponent from experimental data (Esposti, Ferrario, and Signorini, 2008). In this article, we use two approaches for establishing the Hurst exponent: the structure-function and detrended-fluctuation-analysis methods.

Structure Function Scaling
The first technique that we adopt to evaluate the Hurst exponent uses the scaling properties of the structure function (SF). This method is based on the investigation of the scaling properties of a time series directly via the computation of SF being the q-order moments of the distribution of the increments over the time τ (e.g. Takalo and Timonen, 1998;Di Matteo, Aste, and Dacorogna, 2005;De Michelis et al., 2020): The q-order moments of the distribution of the increments are a good quantity to characterize the statistical evolution of a stochastic variable X(t). Hurst analysis examines whether some statistical properties of time series X(t) scale with the time resolution and the observation period [T ]. The Hurst exponent [h(q)] can be defined from the scaling behavior of S(q, τ ) as For q = 1, h(1) describes the scaling behavior of the absolute values of the increments. We are interested in the temporal profile of the Hurst exponent of GCR diurnal variation in Solar Cycle 24 and the solar minima between Solar Cycles 23/24 and 24/25 for 2007 -2019 with resolution of one month. Thus, we consider a moving-scale window of 360 days with the displacement of 30 days. The analyzed data are detrended by removing the linear trend. The scale with T = 360 days is more than ten times larger than the resolution that we want to investigate: 30 days. Here we consider the case for q = 1 (hereafter in this article we designate h(1) = H ) and the scale from 1 to 36 days. This choice permits us to have a reliable estimation of the statistics of the monthly variability of the Hurst exponent [H ] of GCR diurnal variation. The Hurst exponent [H ] is computed through a linear least-squares fitting of Equation 2.

Detrended Fluctuation Analysis
In the next part of the analysis, to determine Hurst exponent for the daily values of the GCR diurnal-variation parameters, we apply the detrended-fluctuation-analysis (DFA) methodology (Peng et al., 1994(Peng et al., , 1995Kantelhardt et al., 2001;Ihlen, 2012). DFA was first proposed by Peng et al. (1994), and from that time the algorithm has found widespread application (e.g. Chen et al., 2002 and the references therein). Moreover, DFA methodology, whose predecessors were the rescaled Hurst interval analysis (Hurst, 1951) and fluctuation analysis (FA) (Peng et al., 1992), has been modified and generalized in subsequent years (e.g. Kantelhardt et al., 2001Kantelhardt et al., , 2002Ihlen, 2012).
In the framework of the basic DFA methodology, the four following steps are executed (Peng et al., 1994;Kantelhardt et al., 2001;Ihlen, 2012). In the first step, the considered time series X of size N (in our case N = 360) is integrated by computing the accumulated departure from the mean of the whole series: where k = 1, . . . , N and X = 1 N N j =1 X(j ). In the second step, the integrated series Y is divided into N s = N s non-overlapping segments v (subseries) of length s. Note that s corresponds to τ mentioned in the context of SF technique (Section 3.1). In the third step, the series Y is locally detrended. More precisely, for a given segment v = 1, . . . , N s of size s, the characteristic size of fluctuation F for integrated and detrended series is calculated by where Y m v (j ) denotes an mth-order polynomial fitted to the Y in segment v. It is worth stressing that various degree polynomial functions can be used (Kantelhardt et al., 2001;Hu et al., 2001). In the analysis presented here, similarly to SF methodology, the first-order (m = 1) polynomial has been applied.
The computation expressed by Equation 4 is repeated over various segment sizes s to provide a relationship between F and s: Finally, a power law is expressed as where the scaling exponent [H ] is the Hurst exponent and is expressed as the slope of a double-logarithmic plot of F (s) as a function of s. One of the important aspects of the procedure matching the exponent H (Equation 6) is the appropriate consideration of the minimal and maximal segment size s. Selected studies considered the maximum box size equal to one-tenth of the signal length (N/10) (Hu et al., 2001) or s > N/4 (e.g. Ma, Yang, andCai, 2005). The limitation of large scale is related to the fact that for very large-scale sizes the function calculated in Equation 4 ceases to be statistically significant as the number of segments N s decreases. Next, for too small scales, systematic variations of the scaling factors may also occur. As many studies showed, the minimum value of the scale parameter should not be less than 10 and generally larger than the order m of DFA (e.g. Phothisonothai and Nakagawa, 2007). In our study, the Hurst-exponent determination based on Equation 6 was performed in a scaling range from 10 to 90 days (≤ N/4). Also, to perform analysis of the diurnal variation of cosmic rays in the years 2007 -2019, the numerical calculations, similarly to SF analysis, were carried out in a 360-day moving window with the displacement of 30 days. It is worth noting that each step of the SF (Takalo and Timonen, 1998;Di Matteo, Aste, and Dacorogna, 2005;De Michelis et al., 2020) and DFA (Peng et al., 1994(Peng et al., , 1995Kantelhardt et al., 2001) procedures presented above has been implemented in MatLab software and checked on data with known Hurst exponents. show a relatively good agreement between the temporal profiles of H obtained by using SF and DFA methods. This can be considered as confirmation of the reliability of the scaling properties of the analyzed data. Of course, one cannot expect one-to-one correspondence between SF and DFA methods, but still they converge, being good tools to study the selfsimilarity problem of the time series.

Results
Detailed analysis of Figure 2 reveals that the Hurst exponents [H ] for GCR diurnal amplitude vary in the range from ≈ 0.3 to ≈ 0.8 for all NMs. One can see that the Hurst exponent is evolving with the 11-year solar-activity cycle with significant variability for    Figure 4 for GCR diurnal amplitude and Figure 5 for the diurnal phase.  Additionally, considering the 37-month smoothing (Figures 4 and 5, bottom panels) for long-term trends, one can see that the level of H is a bit higher for the diurnal phase than for the diurnal amplitude with a tendency for H to be systematically above 0.5 for both. Moreover, one can see the anticorrelation of the Hurst-exponent timelines of both GCR diurnal amplitude and phase with solar activity for 2007 -2014.

Discussion and Conclusions
In this article, we analyzed the scaling properties of fluctuations in the diurnal variation of cosmic rays for Solar Cycle 24 and the solar minima between Solar Cycles 23/24 and 24/25. We applied two methods, structure-function and detrended-fluctuation analysis, to consider the variation of the Hurst exponent as the randomness indicator in GCR diurnal amplitude and phase. The results obtained show that the Hurst exponents [H ] for the GCR diurnal variation parameters vary in the range from ≈ 0.3 to ≈ 0.9 through Solar Cycle 24 and the solar minima between Solar Cycles 23/24 and 24/25. Moreover, the Hurst exponent is evolving with the 11-year solar-activity cycle with significant variability for different GCR diurnal variation parameters. The level of H is a bit higher for the diurnal phase than for the diurnal amplitude, with a tendency for H to be systematically above 0.5 for both. It suggests that the GCR diurnal variation has a more persistent structure than Brownian motion. However, the time series of the GCR diurnal amplitude and phase evolve from a more persistent structure in and near the solar minimum between Solar Cycles 23/24 in 2007 -2009 to a more random character in and near the solar maximum 2012 -2014. Particularly, the results show (Figures 4 and 5, third panels) a transition from a weakly correlated structure (H ≈ 0.63 for the amplitude and ≈ 0.60 for the phase) near the solar minimum between Solar Cycles 23/24 to uncorrelated (H ≈ 0.50 and 0.55 for the amplitude and phase, respectively) near the solar maximum of Solar Cycle 24 of heliospheric dynamics represented by GCR diurnal variation. It is in agreement with the general configuration of the heliosphere, being more regularly structured near solar minimum with the established heliospheric magnetic field and more turbulent near solar maximum. Nevertheless, the temporal profiles of H for GCR diurnal amplitude and phase around the solar minimum between Solar Cycles 24/25 in 2017 -2019 differ in comparison with 2007 -2009, suggesting dependence on solar magnetic polarity. Moreover, during the period 2015 -2019, H -values for GCR diurnal phase differ significantly for SF and DFA methods, with phase causing a different scaling character of the diurnal phase by SF and DFA methods. The above-mentioned scaling features of GCR diurnal variation evolve over 11-year solar-activity cycle, and thus these results could be important for the understanding of the GCR modulation over the solar cycle.