On the fractal properties of cosmic rays and Sun dynamics cross-correlations

We investigate the cosmic rays and ten physical parameters of heliosphere behavior considering the scaling features of their temporal evolution. Our analysis start with the cosmic rays measurements by a neutron monitor station located in Moscow, as well sunspot, Ap index, Na/Np, Alfven mach number, DST index, Kp index, Plasma beta (β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\beta $\end{document}), Proton temperature, F10.7\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$F_{10.7}$\end{document} index and Field magnitude |B|\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$|B|$\end{document} (Interplanetary Magnetic Field, IMF), for the period 1976 to early 2020. Each of these datasets was analyzed by using the Multifractal Detrended Fluctuation Analysis and Multifractal Detrended Cross-Correlation Analysis in order to investigate intrinsic properties, like self-similarity and the spectrum of singularities. The main result obtained is that the cosmic rays time series as well the 10 heliospheric parameters exhibit positive long-range correlations with multifractal behavior.


Introduction
At all times and from all directions our Earth is continuously bombed by a type of radiation that we collectively call Cosmic Rays. This radiation is so called because of our understanding of it as a flow of relativistic particles and radiation in a wide range of energy intervals, which crosses the interstellar medium and eventually on its way hits the Earth. The radiation comes from several sources however one of the main ones is our Sun and also from outer space from distant points in our galaxy and super massive objects found in the universe (e.g. galaxies, galaxy clusters, supernova remnants, etc.).
This generic definition allows us to include electromagnetic radiation in all its wavelengths, neutrinos, and, finally, electrically charged particles. Roughly 90 percent of cosmic ray nuclei are hydrogen (protons) and 9 percent are helium (alpha particles). Hydrogen and helium are the most abundant elements in the universe and the origin point for stars, galaxies and other large structures. The study of neutrinos from the sun, or produced as a result of the collisions of protons from outer space with the atmosphere, has been of crucial importance in understanding the mechanisms of energy generation in the interior of stars.
The chemical composition and abundance of cosmic rays provides important information about their origin and propagation processes from the sources to the Earth. The composition is relatively well known up to an energy of the order of 1 TeV, although not its time dependence. Approximately 98% of cosmic particles are atomic nuclei and 2% are electrons. Of the nuclei, 87% are protons, 12% are helium nuclei and the remaining 1% are heavy nuclei (Grupen et al. 2005;Stanev 2010).
One of the most important characteristics of cosmic radiation is its essentially perfect isotropy at large scale (Getmantsev 1963;Pohl and Eichler 2011;Wentzel 1974). This parameter reflects that the rays arrive from all directions with the same frequency. The transport of Galactic Cosmic Rays (GCR) in the heliosphere (understanding as the spatial region under the influence of the solar wind and its magnetic field) is one of the topics of greatest interest in space physics, given its close relationship with the radiation levels in the interplanetary space of the solar system, the study of which is part of what is known as space weather. This transport is modulated by different physical mechanisms, which in turn can be divided into large-scale processes (Boschini et al. 2022;Engelbrecht et al. 2022;Shlyk et al. 2021): re-lated to the heliospheric magnetic field, transient phenomena: those generated by the ejection of solar matter, or with respect to their duration, which range from several solar cycles, such as processes related to the dynamics of the solar dynamo (Ossendrijver 2003;Charbonneau 2014), to a few hours, such as those produced by transients in small-scale activity that propagate along with the solar wind near the Earth. The 11-year variation is due to the solar activity cycle, and the 22-year variation is due to the magnetic solar cycle. A parameter to indirectly measure solar activity is the number of sunspots that the Sun has on its surface and this parameter presents an anticorrelation with the cosmic ray intensity, that is, during solar activity maxima that is related to sunspot maxima, the cosmic ray intensity presents minimums, and vice versa, during periods of quiescent Sun when there are few sunspots, the cosmic ray intensity is maximum. The effects of solar modulation decrease for higher energy's particles.
On the other hand, solar activity increases and decreases with a period of approximately 11 years. The number of sunspots indicates the level of solar activity. The Sun's matter emissions and electromagnetic fields increase during high solar activity, making it difficult for galactic cosmic rays to reach the Earth. The intensity of cosmic rays is lower when solar activity is high. On the contrary, when solar activity is low, the intensity of cosmic rays is high. Therefore sunspots will increase and the intensity of cosmic rays will start to decrease.
The temporal variations of the GCRs flux are modulated by the intensity of the Interplanetary Magnetic Field (IMF), which has a significant 11-year cycle and other higherfrequency variations. It also shows the monthly variation of the 10.7 cm microwave solar flux, which is a measure of the electromagnetic spectrum (solar irradiance), and the monthly relative sunspot number (R z = k(10g + s), with g the number of sunspot groups in the photosphere, s the total number of spots, and k < 1 a scaling constant).
For the period 1986-2003 for example, the correlation between sunspots and 10.7 cm wave solar flux is greater than 0.95 (Singh et al. 2019). In turn, the GCRs flux has a correlation close to −0.8 (Ahluwalia and Ygbuhay 2015; Sierra-Porta 2018) with the above solar measurements, being in counterphase. A minimum of solar activity is associated with a maximum of GCRs flux and vice versa. Consequently, it could be expected that the correlations found between climate variables and the Sun's dynamics may have, in part, the GCRs as an external astronomical forcing if there is any physical connection.
Traditionally in the state of the art of the relationships that can be made between the GCRs flux and solar dynamics parameters, it is common to find many studies associated to one of the parameters, i.e., the number of sunspots. However, it is known that this index is not the only predictor of long term GCRs variations. Some authors have tried to es-tablish new empirical relationships in the study of cosmic ray modulation from various solar and heliospheric parameters (Belov 2000;Lantos 2005;Chowdhury et al. 2013). Some of these relationships include sunspots, of course, but also other indices such as the monthly counts of the grouped solar flares (GSF), the monthly geomagnetic index (Ap), the monthly flare index (FI), the intensity of the interplanetary magnetic field (IMF), heliospheric current sheet (HCS) and coronal mass ejections (CMEs) (Mavromichalaki et al. 2007;Pohl and Eichler 2011).
In addition, it is known that the series of GCR intensities or counts measured from ground-based monitoring stations, as well as sunspots, exhibit multifractal behavior (Movahed et al. 2006;Sierra Porta and Dominguez 2022), i.e., generation of fractal structures at multiple scales of the variation properties of the time series involved. This phenomenon allows us to characterize such time-varying series in terms of the physical properties of the observed phenomena and even to predict some dynamical aspects related to the predictability and memory effect of the series. The use of spectral techniques for the determination of topological characteristics in time series is framed in the area of multifractal analysis, which has found diverse applications in different fields, such as geophysics and economics. In this case, the methodology aims to find a common scaling behavior in the temporal dynamics of the sun's parameters related to cosmic ray intensities, and to study the persistence of correlations between these variables within different time periods (hours, days, months).
In view of the fact that these previously described relationships have been studied and also the corresponding correlations, for example with some relationships between GCRs and clouds (see for example Christodoulakis et al. 2019). Thus it is important to understand not only the spectral and topological properties of the time series involved but also these same properties for the correlations of various connected phenomena. This is especially important because any relationships that may be established between derived parameters of Sun dynamics and cosmic ray intensities would inherit that fractal or memory characteristic, and knowledge of these measurements and considerations may be important in determining the quality of the predictions and their implication for short and long time scales.
The main objective of this manuscript is to examine in more detail these correlations (and from another view point), for the case of GCRs and some heliospheric and solar dynamics parameters. The methodology described in Multifractal Detrended Cross-Correlation Analysis will be used for this study. In the following sections we describe the data sources and the methodology used as well as the analysis of the results and discussions derived from the application and finally, we make some conclusions of the study.
In this manuscript we present for the first time (based on our review of the existing literature) the application of topo-logical methods for series characterization. In particular our interest is to provide information and quantification of the persistence of correlations between solar dynamics parameters and cosmic rays relevant to substantiate and clarify whether the relationships between these are indeed applicable in the long term.

Signal time series data
For this study we used two data sources. The first data source corresponds to the time series of cosmic ray counts or cosmic ray intensities. The time series considered are obtained from the Pushkov Institute of Terrestrial Magnetism, Ionosphere and Radiowave Propagation of the Russian Academy of Sciences (IZMIRAN, Russian, http://cr0. izmiran.ru/common/links.htm) (Belov et al. 2005;Gaidash et al. 2017) which is a scientific institution of the Russian Academy of Sciences, with original data from on-site neutron monitor stations of the Neutron Monitor Database (NMDB, https://www.nmdb.eu/) (Mavromichalaki et al. 2010(Mavromichalaki et al. , 2011 which provides access to Neutron Monitor measurements from stations around the world. The goal of NMDB is to provide easy access to all Neutron Monitor measurements through an easy to use interface. NMDB provides access to real-time as well as historical data.
Depending on the location of the detectors, the intensity of GCRs detected will be affected by the Earth's magnetic field. Detectors located at low latitudes register a lower number of particles, while detectors located near the Earth's magnetic poles register a higher amount. This latitudinal effect is known as the geomagnetic shielding. The number and distribution of GCR observatories helps us to map the arrival of cosmic radiation to our planet.
The data used are hourly resolution neutron counts for three neutron monitor. Moscow neutron monitor is situated in Troitsk city, Moscow Region. It is operated by Pushkov Institute of Terrestrial Magnetism, Ionosphere and radio wave propagation (IZMIRAN) of Russian Academy of Science. The detector 24NM64 with coordinates latitude 55.47 • N and longitude 37.32 • E, is located at an altitude of 200 m.a.s.l with a vertical cutoff rigidity (1965) of 2.43 GV, operating uninterruptedly from 1958 to the present has received several upgrades but with a remarkable data continuity. Oulu (9-NM-64, Chalk River tubes) neutron monitor (latitude 65.05 • N and longitude 25.47 • E, 15 m.a.s.l, vertical cutoff rigidity ∼ 0.8 GV) located in Finland Operated by Sodankyla Geophysical Observatory of the University of Oulu, in operation since April 1964. The last one used is Novosibirsk neutron monitor (24NM64) is operated by Institute og Geology and Geophysics Russian Academy of Science, latitude 54.48 • N and longitude 83.00 • E, is located at an altitude of 163 m.a.s.l with a vertical cutoff rigidity of 2.91 GV, operating uninterruptedly from 1971 to the present.
For this study we are set a data range from 1986 to actual time corresponding with solar cycles 22-24.
The other source of data, corresponding to predictors of cosmic ray intensities, comes from solar dynamics parameters. To investigate these correlations in this case we obtain physical parameter information through high-resolution OMNI (1-min, 5-min): solar wind magnetic field and plasma data in the Earth's Bow Shock Nose (BSN), also geomagnetic activity indices and 5-min energetic proton fluxes (NASA: spacecraft interleaved near-Earth solar wind data, https://omniweb.gsfc.nasa.gov/ow_min.html) (Vokhmyanin et al. 2019;King and Papitashvili 2005).
All of the data to which this interface and its multiple underlying interfaces provide access have in common that they are relevant to heliospheric studies and/or provide solar wind input data for solar wind-magnetosphere coupling studies, and that they reside at GSFC/Space Physics Data Facility. This interface was built in 2008 to provide a more integrated interface to the numerous plasma, magnetic field, and energetic particle data sets relevant to heliospheric studies and residing in Goddard's Space Physics Data Facility. It provides users with underlying interfaces (OMNIWeb, CO-HOWeb, FTPBrowser, CDAWeb, Helioweb, spdf/ftp) that provide diverse functionality for various data sets.
From this dataset we extract 10 different parameters. The sunspot number (R z ); the solar radio flux at 10.7 cm (2800 MHz, F 10.7 -index) is an excellent indicator of solar activity, often referred to as F10.7, is one of the oldest records of solar activity, the radio emissions originate high in the chromosphere and low in the corona of the solar atmosphere, it is well correlated with sunspot number, as well as with a number of ultraviolet (UV) and visible solar irradiance records; the ap-index (Ap) provides a daily average level of geomagnetic activity; the interplanetary magnetic field (IMF) is a vector figure with a three-axis component, two of which (B x and B y ) are oriented parallel to the ecliptic and are not important for auroral activity, so here we use the B z component and the |B| modulus; the Alpha/Proton ratio (Na/Np); the Alfven mach number (AMN) is that ratio specifically with the velocity of Alfven waves in a plasma; the disturbance storm time (DST) index, which is an index of magnetic activity derived from a network of geomagnetic observatories near the equator that measures the intensity of the globally symmetric equatorial electrojet (the "ring current") the Kp index is the global geomagnetic activity index that is based on 3-hourly measurements from ground-based magnetometers around the world, used to characterize the magnitude of geomagnetic storms, it is also an excellent indicator of disturbances in the Earth's magnetic field and is used to decide whether to issue geomagnetic alerts and warnings to users who are affected by these disturbances; plasma beta parameter (β = Np[(4.16T /10 5 ) + 5.34]/B 2 , where T is the proton temperature and B is the magnetic field); and the solar wind proton density (Proton Density, N/cm 3 ).

Methods and strategy
There are many examples in nature in which stochastic processes exhibit cross-correlation, which means that a pair of series are related to each other while retaining a high correlation for a given offset of the series (Yoo and Han 2009;Lynn and Lynn 1973;Welch 1974). The techniques available to measure correlation make use of convolution techniques to measure the optimal offset at which correlation is maximized. This allows us to understand phenomena of simultaneity of the processes as well as the precursor characteristics, in terms of the dynamics of one series in terms of the other. In both one-dimensional and two-dimensional systems this is important in signal theory to search for and locate patterns to understand the proper and articulated dynamics of dynamical systems.
The usual methods for measuring the cross-correlation of two related signal series usually assume seasonal and stationary behavior. In essence, if both series u and v In the real world, however, this property cannot be preserved because the signals are by nature non-stationary.
Recently, new developments for measuring cross-correlation have been studied to overcome the limitations in the series signals discussed above. The methodology, as well as the associated algorithm, is collectively called Multifractal Detrended Cross-Correlation Analysis (MFDCCA) (Podobnik and Stanley 2008;Zhou et al. 2008), which can be understood as a generalization for the cross-correlation of two series of the Multifractal Detrented Fluctuation Analysis (MFDFA) method (Kantelhardt et al. 2002). The algorithm also allows to measure in terms of fractality analysis, the persistence of correlations between two series. Its simplicity is remarkable as well as its application in many fields of science, for example in financial stock markets (Preis et al. 2013;Bordino et al. 2012;Junior and Franca 2012;Moat et al. 2013), geophysical process (Marinho et al. 2013;Shadkhoo and Jafari 2009), environmental pollution problems (Zhang et al. 2015;Tzanis et al. 2020), energy availability (Zhuang et al. 2014;Pal et al. 2014), medical science (Wang and Zhao 2012;Ghosh et al. 2014), among many others applications.
The mechanism is quite simple. It needs two time series to be calculated u, and v. The idea consist on divide both signal series into N − n overlapping time intervals, with N and n the length of series and partition intervals, respectively, each one containing n + 1 values. The covariance is then computed, where s = 1, 2, . . . , N − n andû s , v s are regression fit over each interval s with a polynomial or order to decide. Finally, the fluctuation function is calculated summing over all the overlapping time inter- . Recently, a cross-correlation index ρ(n) has been introduced in order to make the interpretation of detrented cross-correlation analysis clearer (Zebende 2011), ρ(n) = F 2 u,v (n)/ (F u,u (n)F v,v (n)). Similar to the standard crosscorrelation coefficient, ρ(n) = 1 indicates a perfect crosscorrelation, while ρ(n) = −1 means a perfect anti-crosscorrelation. If ρ(n) = 0, there is no cross-correlation between the two time series at the scale.
The qth-order detrended covariance is calculated as generalization of F u,v (n) as: for q = 0, but for q = 0 ln F u,v (n, s) . (2) Finally, F 2 u,v (n) is expected to scale as n H , where H is the Hurst exponent, and its value can indicate if a process is persistent or anti-persistent: (a) H < 0.5 indicated the process is anti-persistency, i.e. the correlations under study are anti-persistent and tends to decrease (increase) after a previous increase (decrease), an anti-persistent correlations appears very noisy; (b) H = 0.5 imply a uncorrelated process; (c) 0.5 < H < 1 correspond to a persistent, correlations has been increasing (decreasing) for a period T , then it is expected to continue to increase (decrease) for a similar period, persistent processes show long-range correlations and exhibit relatively little noise; and (d) H > 1 nonstationary process, stronger long-range correlations are present.
If H (q) is approximately constant for all values of q, the time series is said to be monofractal, that is, it exhibits the same scaling at all scales. If H (q) varies significantly, the time series is multifractal, with different scaling at different scales. The different scaling values are better described by where α = H (q) + q dH (q) dq . The spectrum is concave down, and the wider it is, the more multifractal the time series.
Essentially, f (α) is the dimension of points with pointwise dimension α. Another way to relate each of the singularity measures of the series more compactly is through a Legendre transformation For the application of the method we use the algorithm developed in python for this purpose and described in Bianchi (2020).

Individual series analysis
The series studied for both GCRs and solar parameters have a length of 403224 sample points in hourly frequency and starting from 1976-01-01 to 2021-12-30. In the case of the GCRs series in case of Moscow neutron monitor, there are 1759 missing values, less than 1% of data sample. For other neutron monitor series the characteristics of data are similar to the first one. In this section the following analysis is taking only into consideration Moscow monitor. In the case of the heliospheric parameters there is 62053 missing values in average with a standard deviation of 42419 values, it's represent less than 1% of data. All series have also been normalized and scaled. For this estimator, each characteristic is translated individually so that it is in the range between zero and one: X scaled = (X − X min )/(X max − X min ) + X min . Figure 1 shows the correlation coefficients between GCRs and the heliospheric indices considered. All the solar dynamics signals and parameters were chosen from the maximum correlations found for these series in correspondence with GCRs. The correlations are sufficient to consider the application of the MFDCCA method. However, additionally, we also find sufficiently strong correlations between them, in particular, sunspots are strongly related to F 10.7 index, Interplanetary Magnetic Field (Field Mag. Avg. |B|), Na/Np and Alfven mach number; also Ap index, DST index, Kp index and Protron temperature are closely correlated. Figure 2 shows the comparative dynamics of the cosmic ray series and its heliospheric predictor counterpart. As can also be seen in Fig. 1, the parameters exhibiting positive correlation are Alfven mach number, DST index and Plasma beta (beta), while those exhibiting anti-correlation are R z , F 10.7 index, Ap index, Field Mag. Avg. |B|, Na/Np, Kp index and Proton temperature.  To make a rigorous analysis of the correlations we use a cross-correlation study between each of the solar parameters with the cosmic ray count series, from which we calculate the corresponding time lag to account for the time interval in which the maximum correlation occurs. For the case of the number of sunspots the maximum correlation is −0.89 with a time lag of about 2 months. This has been studied in multiple references (see for example Sierra-Porta 2018; Kane 2014). The variability of cosmic ray counts measured by ground-based detectors are highly influenced by solar activity but not at a given time, but rather over a period that can generally be longer (see for example Belov 2000;Nagashima and Morishita 1979). In other words, it takes a long time, in most cases, for the sun to produce changes in the interstellar medium to produce significant changes in the cosmic ray rate. One explanation of this effect according to the Nagashima and Morishita (1979) point of view is: "When the Sun's polar magnetic field is parallel to the galactic magnetic field, they could easily connect with each other, so that the galactic CRs could more easily enter the heliosphere along the magnetic lines of force, compared to those in the antiparallel state of the magnetic fields...... The transition between the two states has a delay in polarity reversal, depending on the rigidity of the observed CRs; such as about 1 year for a neutron monitor and about 3.5 years for low rigidity components.". More recently, however, it is known that the 22-year variability of the correlation lag is explained by the action of heliospheric drifts (for a current review see Koldobskiy et al. 2022;Wang et al. 2022). In Table 1 we show the result of cross-correlation analysis.
In order to remove seasonal variation and non-seasonalities for the application of the method and algorithm described by MFDCCA, we specifically take the studied series and make N = 30 partitions of the series with linearly spaced window sizes between 6720 and 100806, corresponding with L d /60 and L d /4, respectively, with L d the data length. This implies that the minimum and maximum windows are between 9 months and 7 years, respectively. Additionally, for the qth-order detrended covariance we test q = [−6.0, 6.0] spaced by 0.5 and for or results we use fifthorder polynomials for fitter in each interval. From this we then compute the variability measure of the fluctuations for each interval and for each q in order to investigate the empirical relationship with the measures of the scaled window widths as n H , where H is the power-law scaling exponent that can be considered as the important parameter to measure self-similarity and that can be associated with long-range fractal properties in the studied series. Figure 3 shows the corresponding fluctuation function F (n) for various time scales in the series partitions for each of the solar dynamics indices as well as for the cosmic ray count series at the Moscow station. The fluctuation function for the observed range of scales exhibits persistent longrange behavior with Hurst exponent greater than 0.5 for the case of the cosmic ray time series signals for all orders of covariance q (the exponents in Eq. (1)). This has been found previously in Sierra Porta and Dominguez (2022) and furthermore this behavior has also been linked to the cutoff rigidity of each of the stations, which means that for various time periods the series exhibits increases (or decreases) followed by decreases (or increases) over a similar period of time. This is closely linked to the quasi-seasonal characteristic of the series mediated largely by the dynamics of the heliosphere. For the case of Field Magnitude Average |B|, Na/Np, Alfven mach number, DST index, Kp index and Proton temperature, the series exhibit anti-persistence over the whole range of q values, while in the case of sunspots (R z ), F 10.7 index, Ap index and Plasma beta (β), the series exhibit persistence for values of q < −1 and anti-persistence for q > −1.
The next step is to study the singularity spectrum for all the series considered using the MFDFA technique, which is usually based on the scaling and Legendre transform associated with the q-th order moments that depend on the signal lengths, as mentioned in the method description section. In Fig. 3, the spectrum of the fluctuation functions for the most constant part of the fluctuations and time scales has been shown. However, this is generally not always constant. When the spectrum of the fluctuations is not constant over some range of scales the series is said to exhibit fractal behavior. Figure 4 shows the singularity or multifractality spectrum associated with each of the individual series using Eq. (3). Figure 4 shows a corresponding vertical line with the maximum value of f (α) for some associated q and a strip that determines the corresponding α(= α max − α min ) between the extremes of the spectrum and also accounts for the fluctuation asymmetry for the singularity spectrum. When the curves are rightward asymmetric, i.e., the curve fluctuates more to the right towards smaller values of f (α), the fluctuations are said to be sensitive to small scales, conversely, when the curve fluctuates to the left towards smaller values of f (α) then the fluctuations dominate at large scales. These develop suggests that the fluctuations are multifractal for all the series studied. Indeed, for the cosmic ray case this has been confirmed in Sierra Porta and Dominguez (2022) even for several cosmic ray counting stations at different cutoff rigidity.
In none of the series studied the spectrum fluctuates similarly on both sides of the maximum value, this is confirmed as shown in Fig. 3.

The cross-correlations between GCRs and solar dynamics parameters according with MFDCCA
For the application of MFDCCA we have used the same parameters as for the individual series, this means that we have divided our series signals into N = 30 intervals of similar widths as above, furthermore, in the results shown below we performed the multifractality analysis for the correlations of GCRs compared to the heliosphere parameters and indices. We use here three dataset of neutron monitors described previously.
The main result of the analysis described in the previous section is that the cosmic ray series for the Moscow monitoring station exhibit positive long-range autocorrelation with pronounced multifractal behavior. This may also be true for some other parameters of solar dynamics such as Na/Np, Alfven mach number and DST index.
Here we now attempt to correlate the behavior of cosmic ray dynamics with heliospheric parameters that have been shown to have an impact as precursors or predictors of these dynamics. The Spearman correlation coefficients have been shown in the preceding section, so there is evidence of relationship and causality between each of the parameters. Evidently sunspot number (r = −0.82) is a very good predictor but so are for example Field Mag. Avg. |B| (Interplanetary Magnetic Field, IMF, r = −0.84), Na/Np (r = −0.81), Kp index (r = −0.71) and Beta Plasma (β) (r = −0.71), all exhibiting anti-correlation. The cross-correlation measurement result yields lag times of 68, 4, 0, 1 and 28 days, respectively.
To continue with the analysis we measure the fractal spectrum of each correlation. Figure  In all cases for the studied indices and all neutron monitor stations the slopes for the fit of the fluctuation functions are slightly higher for the case of F 10.7 Index radio emissions which has proven very valuable in specifying and forecasting space weather and also correlated with cosmic ray intensity. Additionally, for the Moscow and Novosibirsk stations the results are very similar and slightly different (lower in intensity) than for the Oulu monitor station. An interesting result that follows from this observation has to do possibly with the dependence that could be exhibited for different geomagnetic cutoff rigidities at the geographic points of each station. (1) and H (q) represents the Hurts exponent or slope in the linear fit for each q A reasonable hypothesis for the moment that emerges is that the slopes in the fluctuation functions may have dependence on the geomagnetic cutoff rigidity.
In all cases the correlations exhibit characteristics of nonstationary processes or strong long-range correlations, all Hurts exponents occur for values greater than 1 (H (q) > 1) for all orders of covariance q. This can be confirmed in the continuous spectrum of H (q) in Fig. 6 in the left panels for Moscow station and Fig. 9 and 10 for Oulu and Novosibirsk, respectively, in Sects. A.1 and A.2.
Additionally, the singularity spectrum does exhibit significant differences between all correlations (see Fig. 6 right panels). It is possible to see that none of the spectra fluctuate symmetrically with respect to the maximum of f (α). This implies long-range correlations for sunspot, Ap index, Na/Np, Alfven mach number, DST index, Kp index, Plasma beta (β) and Proton temperature, while short-range multifractal correlations for F 10.7 index and Field magnitude |B|.
The above indicates that cosmic ray counts are not only connected to all heliospheric parameters but also that systematically a growth (or decrease) of GCRs is associated with a decrease (or growth) in the dynamics of some indices for negative correlations and for positive correlations growth (or decrease) of GCRs occurs after growth (or decrease) of heliospheric parameters for the case of positive correlations. However, in almost all cases this behavior has strong multifractal characteristics of non-stationary processes with longrange correlations.
To complete our analysis, the following Table 2 presents the values of f max (α == α max ) and f max (α == α max ). As can be seen, for each neutron monitoring station, the values follow similar trends, however we can notice a slight tendency to higher values for the Oulu monitor than for the other two. The Moscow and Novosibirsk monitors have similar cutoff rigidity parameters but Oulu has a significantly lower cutoff rigidity. This could confirm the hypothesis of a trend of the multifractal intensity measurement of the fluctuation spectrum with the geographical position of the monitors. In all cases, the maximum value of f (α) corresponds to values very close to q = 0, the measure of the singularity of the multifractal spectrum is not symmetric but in all monitors the trends are maintained.

Conclusions
We investigate the intrinsic self-similarity and the spectrum of singularities of the cosmic ray time series in one of the monitoring stations of the Russian monitor network (Moscow), for the period 1976-2020, as well as some indices of the heliosphere that determine high correlations with cosmic rays, that is, they can and have been used as precursors or predictors of cosmic ray activity measured on the ground. The techniques used were MFDFA and MFDCCA. The main finding from the results maintains that the temporal evolution of cosmic rays presents long-range positive autocorrelations and multifractal behavior, this also occurs for some indices of solar dynamics studied.
After applying the technique, the optimal range of multifractality for small time scales has been detected in the order of ln(n) ≈ [8.5, 11.5], which is expected since we know that both cosmic rays as some of the heliosphere indices have quasi-seasonal characteristics with periods associated with cycles of repetition of general behavior, in the case  of sunspot for example of about 11 years. Furthermore, the spectrum of fluctuations and their corresponding Hurst exponent vary significantly in large windows spanning several local periods with small and large fluctuations (i.e., negative and positive q, respectively), and thus will average their differences in magnitude. However, common features of multifractality were revealed for both positive and negative q values. Next, we tried to compare the time series of GCRs with the time series of the physical parameters and properties of the heliosphere, detecting statistically significant correlation and multifractal behavior in all correlations. It is important to point out that many of the investigations related to the possible association of GCRs and heliophysical parameters have improved the understanding of the phenomena that occur in the heliosphere and their relationships are known in terms of models that cause causality, however, until now there had been no These associations and their persistence characteristic have been studied from the topological point of view.
The results obtained could be used to improve the existing models of regressions between the series studied, including some parameters derived from the intrinsic multifractality of each pair of correlations between the series, in particular, this could be interesting for the inclusion of these values in models of machine learning that could improve the effectiveness and predictability of the estimates. For example, the results mentioned above could give more information about the difference in the evolution of solar activity during odd and even solar activity cycles (Durney 2000;Usoskin et al. 2001).
As a future work arising from and inspired by the results of this analysis, the authors are developing for now a much more detailed study to consider the dependence of the geomagnetic cutoff rigidity of neutron monitor stations at various latitudes.

Final comments. Supplementary data
Although this study was conducted without any type of financing from any entity and/or organization, the authors would like to express their gratitude to the Universidad Tec-nológica de Bolivar (UTB) and the Research Department for their constant support and providing a computer to perform the calculations in this study. Additionally, we would like to thank the Dirección de Investigación of Universidad Tecnológica de Bolivar for their support in revising the style of the manuscript for an improved presentation.
For a review or for replication and reproduction purposes of our results, a version of the dataset (Sierra Porta 2022) used for this study and the results obtained can be found in above reference.  The singularity spectrum f (α) for the correlations versus singularity strength α for the all examined times series permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.