Modeling the Time Delay Problem of Galactic Cosmic Ray Flux in Solar Cycles 21 and 23

We present a 2-dimensional time-dependent simulation based on Parker’s transport equation (PTE) describing the propagation of energetic astroparticles – cosmic rays (CR) in the heliosphere. PTE is a second order partial differential equation containing four major processes: diffusion, convection, drift, and adiabatic cooling responsible for modulation of the CR ﬂux in the heliosphere. We implement in the numerical simulation, a few physical parameters as the tilt angle, δ , of the heliospheric current sheet (HCS), module, B , of the interplanetary magnetic ﬁeld (IMF), and variations of drift effect of the CR particles with solar activity (SA). Our approach is the implementation of two independent parameters (proxies), γ and ν . The parameters γ and ν are calculated from various independent sources, γ , from neutron monitors (NM) daily data, and ν , using the IMF hourly data. The solutions of PTE obtained from our numerical model are compared with the variations of the CR ﬂux recorded by NMs. We prove the existence of a varying delay time (DT) between the changes of galactic cosmic ray (GCR) intensity and the parameters characterizing SA. Based on our investigation, we obtained different DTs in Solar Cycles 21 and 23. We conclude that the calculated DTs, after comparison with the observed DTs, are useful parameters for the study of GCR transport in the heliosphere.


Introduction
The problem of delay times between the intensity of galactic cosmic rays, I , and other parameters is relevant regarding different periods of solar activity (SA). Mechanisms of GCR transport in the heliosphere are based on four major processes: diffusion, convection, B M. Siluszyk adiabatic energy changes, and drift caused by the gradient and curvature of the interplanetary magnetic field. The spatial diffusion is caused by random turbulence scattering of the IMF in the heliosphere. All these processes are vital, but their contribution in modulation is different and depends on the SA. The existence of a delay time (DT) been a generally recognized fact for many decades, studied by various scientists though interpreted differently.
Many studies of the DT problem showed that DTs are more pronounced in odd-numbered cycles (Burlaga, McDonald, and Ness, 1993;Parker, 1963) than in even-numbered cycles (Jokipii and Thomas, 1981). Scientists explain this phenomenon by drift effects (Cliver and Ling, 2001;Chowdhury, Kudela, and Dwivedi, 2013;Chowdhury and Kudela, 2018;Mavromichalaki, Belehaki, and Rafio, 1998). Based on the drift theory (Jokipii, 1971) of modulation of GCR in solar cycles with A > 0, where A is the polarity state of the solar magnetic field, a drift stream of GCR caused by the gradient and curvature of the IMF preferentially enters the heliosphere from the polar region and is ejected outward along the equatorial current sheet. The reverse situation occurs in periods when A < 0, the drift stream of GCR enters the heliosphere along of the heliospheric current sheet (HCS) and leaves it near the poles of the Sun. In periods when A < 0, GCR particles are more affected by propagating diffusive barriers associated with SA and the waviness of the HCS, resulting in large DTs. The opposite situation is observed in periods with A > 0. GCR particles are less affected by propagating diffusive barriers which cause small DTs or no DTs.
Several authors have discussed in detail the phenomenon of GCR modulation and the DT problem between the GCR intensity, I (GCR), and various parameters characterizing conditions in the heliosphere. For example Chowdhury, Kudela, and Dwivedi (2013) mentioned that GCR modulation is a very complex phenomenon occurring throughout the heliosphere and depending on several factors. No single solar parameter can account for the GCR intensity variations. The IMF magnitude, B, plays a vital role in the GCR modulation, because the Larmor radius, an important parameter determining particle transport in space, is inversely proportional to the strength of the IMF. An increase in the IMF strength decreases the Larmor radius and the diffusion coefficient, leading to the increase in GCR modulation. It is remarkable, however, that the highest correlation between B and I (GCR) takes place at zero DT in Cycle 23. Chowdhury, Kudela, and Dwivedi (2013) propose that there is no significant DT between B and I (GCR), because the local disturbances, such as CMEs and traveling shocks, which are injected into the inner heliosphere, dominate over the effects of merged interaction regions (MIRs) and global merged interaction regions (GMIRs) operating at large distance in the outer heliosphere in this interval. It is possible that some indices and parameters such as the solar sunspot number (SSN), the solar sunspot area (SSA), and the solar flare index (SFI) represent global effects, whereas others such as B and Ap represent local effects limited to the space near or at the Earth. This is consistent with Usoskin et al. (1998) suggestion that GCR modulation appears clearly correlated only with the global indices because of the complicated transport of GCRs in the heliosphere.
Opposite to the above work on the DT problem, in our last article Iskra et al. (2019) we examined the whole period from 1959 to 2014 which includes two and a half 22-year solar magnetic cycle (SMC), and divided it into five periods corresponding to different signs of the polarity state of the solar magnetic field, A. To examine the DT in a more detailed way each of the five periods with a given sign of the IMF was divided into two sub-periods: ascending and descending periods of SA. Selection of other periods to determine the dependence of DT on the IMF sign showed some inconsistencies with the results of other researchers. We observed a 22-year periodicity of the DTs in the relationship of SSN with I (GCR) and SSN with B, but we did not observe it between the remaining parameters. Based on the results found by Iskra et al. (2019), we can conclude that the structure of IMF turbulence from the minimum to the maximum of SA radically changes. In addition, these structures are different for the periods with A > 0 and A < 0 due to the polarity dependence of the parameters ν y and ν z (Siluszyk, Iskra, and Alania, 2015a;Siluszyk et al., 2015bSiluszyk et al., , 2018. In periods with A < 0, GCR particles are more affected because of the propagation of the appropriate size turbulence-structures of the IMF associated with SA and changing induction of B, resulting in large DTs. The drift of GCR caused by the gradient and curvature of the IMF is an additional factor that strengthens this phenomenon.

Relevant Parameters in PTE
The long-term changes of the CR flux can be described based on the PTE (Parker, 1958). The diffusion flux plays a main role together with convection and adiabatic cooling fluxes in long-term variations of the CR intensity (Dorman, 2006). Nevertheless, this assumption requires the answer to an important question, which of the parameters characterizing SA and the solar wind (SW) can account for the changes in the diffusion of the CR particles.
To answer this question we have to consider all other available arguments. The first one is the parameter γ which characterizes the changes of the rigidity, R, spectrum of CR flux intensity variation given by a formula δI I ∝ R −γ , where I is the intensity of CR. The second parameter is the exponent ν of the power spectral density (PSD) of the IMF fluctuation (PSD ∝ f −ν , where f is the frequency) The exponent ν has been calculated in the frequency interval f = f 2 − f 1 of the resonant frequencies, with f 1 = 1 · 10 −6 Hz, f 2 = 4 · 10 −6 Hz, that are responsible for the scattering of GCR particles with the rigidity range detected by NMs (Alania, Iskra, and Siluszyk, 2003;Siluszyk, Iskra, and Alania, 2014).
Hence, we have two very important physically realistic parameters γ and ν calculated from independent sources. In our article, solutions of PTE obtained in our numerical investigation are compared with data of the CR flux measured by Oulu NM. So, both independent parameters, γ and ν, can be considered like crucial indicators to describe the behavior of CR in space (Alania et al., 2001;Siluszyk, 2008a, 2010). We have found a strong inverse correlation between γ and ν in Solar Cycles 21 (1976-1986) and 23 (1996-2008) (Siluszyk, Wawrzynczak, and Alania, 2011;Siluszyk et al., 2018).
To estimate the role of changes of the diffusion coefficient, K, in long-term changes of the CR flux registered by NMs we consider a parallel diffusion coefficient (DC), k II , having the general form (Jokipii, 1971;Iskra, Siluszyk, and Alania, 2015).
Based on experimental data the diffusion parameter, α, can be exchanged by the parameters γ and/or ν (see Figure 1). The relations between the physical parameters α, γ , and ν are presented in detail in previous articles (Hedgecock, 1975;Alania et al., 2005;Alania, Iskra, and Siluszyk, 2008a;Alania, Iskra, and Siluszyk, 2008b;Alania, Iskra, and Siluszyk, 2010;Siluszyk, Iskra, and Alania, 2015a;Siluszyk et al., 2015b). We call DT to a time interval between the maximum of the CR intensity changes on the one side and any parameters of SA and SW on the other. As mentioned in the introduction, DTs depend on many parameters. DTs are influenced not only by particle drifts, but also by the structure of the magnetic turbulence at various periods of SA, which are determined by the characteristic of diffusion and of coronal mass ejection (CME) activity. As a rule, in this article, we have examined the one to one correspondence between the maximum value Figure 1 The diagram shows the relations between the parameters: α, γ and ν. In this diagram D corresponds to I . of the CR flux and the minimum value of any SA parameters (or vice versa). Usually, to simulate the PTE describing the behavior of energetic CR particles is a difficult problem due to the complexity of electromagnetic processes in the heliosphere. These processes are well reflected in symmetric and asymmetric parts of the 3-dimensional (3D) generalized anisotropic tensor of CR diffusion (Alania, 2002), where K 11 = K II cos 2 δ cos 2 ψ + β cos 2 δ sin 2 ψ + sin 2 δ is the angle between B of the IMF and the B r component of B, δ = arctan(B ϑ /B r ) is the angle between the B of the IMF and B θ component of B in the spherical coordinate system. The ratios β = K ⊥ K and β 1 = K d K are calculated based on the coefficients K , K ⊥ , K d , which are the parallel, perpendicular, and drift diffusion coefficients, respectively. When modeling the variations of CR flux, one can find several problems, i.e. when modeling short term variations of CRs (as the anisotropy, Forbush decreases and the 27-day variations of CR flux) it is essential to consider a 3D time-dependent model of the PTE, while to model 11and 22-year changes it is good enough to perform a 2D time-dependent model of the PTE.
The goal of our article is to develop a 2D model, which in the best case will be adjusted to a real situation in the heliosphere. We compare results of modeling and experimental data of the Oulu NM for two 11-year periods (1976-1988) and (1996-2009) with similar heliospheric structure with respect to the polarity of the Sun's global magnetic field. For this purpose in the next section we describe the experimental data.

Experimental Data and Their Characteristics
Changes of the normalized (to the maximum intensity in 2010) monthly intensity I of the CR registered by Oulu NM and SSN in the period 1964-2016 for the 11-year Cycles 21 and 23 are presented in Figure 2. Figure 2 shows that, in both parameters I (CR) and SSN, long-period changes for odd and even solar cycles are observed, which also depend on the Sun's global magnetic field direction (peak or plateau in the maximum CR intensity). Analysis of Figure 2 shows that I and SSN behave almost similarly in the odd Solar Cycles 21 and 23.
A general difference is that Cycle 23 is almost two years longer than Cycle 21, due to an anomalous elongated minimum. For modeling we have chosen two periods -Solar Cycle 21 (1976-1988) and 23 (1996-2009). To reliably determine DTs between I (CR) and inverted SSN we present the changes in I (CR) versus inverted SSN in Figure 3 for Cycle 21 and in Figure 4 for Cycle 23.
In this article we consider realistic models including the above-mentioned physical approaches. The aim is to find whether the expected changes of density in the CR particles are shifted or not, in relation to the changes of the experimental data of the CR flux.
We present in Figures 5 and 6 the time profiles of the IMF, B(t), and the calculated parameters γ (t) and ν(t). In our previous article  we determined the relation between γ and ν. In general, the structure of the IMF turbulence is different in different SA cycles. Based on the theory of CR particle scattering in space, a natural phenomenon is the existence of an inverse correlation between γ and ν during Solar Cycles 21 and 23. Figure 6 shows that γ and ν are anti-correlated in the periods 1976-1988 and 1997-2009.

Modeling the Long Period CR Variations in 1976-1987 and 1997-2009
Our aim is to compare the results of modeling of the CR transport for the 11-year Solar Cycles 21 and 23. Those cycles were chosen because they have similar changes in CR intensity regarding the polarity of the Sun's global magnetic field (decreasing intensity of  CR flux for positive polarity A > 0 and incrising intensity of CR flux for negative polarity A < 0 epochs of SA). Moreover, in the considered cycles we observe an inverse correlation between the exponent ν of the PSD of the IMF turbulence and the rigidity spectrum exponent γ . To model the 11-year variations of GCRs we use a non-stationary PTE (Parker, 1965;Manuel, Ferreira, and Potgieter, 2014).
where N , R, τ , U , v d are the omnidirectional distribution function, the rigidity of the CR particles, the time, the SW speed, and the drift velocity, respectively. Then, we introduce dimensionless values, density f = N N 0 , time t = τ −τ 0 τ S −τ 0 rigidity R = R i 1 GeV , and distance r = ρ ρ 0 . In the previous equations N 0 = 4πI 0 , where the flux of CR for protons is I 0 = 21.1T −2.8 /(1 + 5.85T −1.22 + 1.18T −2.54 ) in the border of the heliosphere (Webber and Lockwood, 2001;Manuel, Ferreira, and Potgieter, 2014); T is the kinetic energy in GeV  ), e is the charge of the electron. Furthermore ρ, ρ 0 , τ 0 , τ s , τ are the radial distance, the distance to the border of the heliosphere, the year of the beginning of the cycle, the year of the end of the cycle, and the current year, respectively.
We take the distance to the border of the heliosphere as ρ 0 = 100 AU, and the SW speed U = 400 km s −1 . Equation 1 in the 2D spherical coordinate system (r, θ ) using the dimensionless variables f , t , and R, has the form: where coefficients A i (r, θ, R, t) for i = 1, 2, . . . , 8 have the following form: The anisotropic CR diffusion tensor has a form K ij = K (S) ij + K (A) ij and consists of the symmetric and antisymmetric parts K (S) ij and K (A) ij . We consider in the model the drift velocity of the CR particles as, ν D,i = ∂K (A) ij ∂x j (Jokipii, Levy, and Hubbard, 1977;Caballero-Lopez and Moraal, 2004). We take into account the ratios β = K ⊥ K and β 1 = K d K for the CR particles of rigidities R ≥ 10 GV, where K , K ⊥ , K d are the parallel, perpendicular, and drift diffusion coefficients, respectively. We also take into account changes of ωτ 1 in the heliosphere, based on Parker's spiral magnetic field. We suppose that at the Earth's orbit for CR particles with rigidities R = 10 GV, ωτ 1 = 3 and the ratios β and β 1 are where ω is the Larmor frequency rotation of the particle with a charge q and mass m under a regular IMF, B, and τ 1 is the average time between collisions of CR particles with the IMF turbulence. The product ωτ 1 = 3 was estimated based on the data anisotropy of GCR particles in various sectors of IMF. The quasi-linear theory expresses the energy or rigidity dependence as, K(R) ∝ R α , where α is a diffusion parameter, which is satisfied for rigidities R > 10 GV (Rossi and Olbert, 1970;Bieber, 2003;Shalchi, 2009;Shalchi and Schlickeiser, 2004;Iskra, Siluszyk, and Alania, 2015). We can assume that both proxies, γ and ν, can be used adequately to describe the state of the heliosphere (see Figure 5). Hence, based on Equation 2 we construct two different models for both proxies (there are four cases): Model I: We use as a proxy the parameter γ , where the parallel DC can be written as The term describing the energy or rigidity dependence has a form K(R, γ (t)) = R γ (t) .
Model II: We construct the model using the parameter ν, where the parallel DC can be written as The term describing the energy or rigidity dependence is K(R, ν(t)) = R ν(t) . Other components are K 0 = 1.9 · 10 19 cm 2 s −1 , K(r) = 1 + 50r. K(t ) is a function determining the change of DC with time, during the solar magnetic cycle.
In Table 1 we present the functions normalized to the time interval t ∈ [0, 1] used in models describing the PTE for Solar Cycle 21 and 23, respectively.
We have obtained curves based on the experimental data provided by the Neutron Monitor Data Base (NMDB) (http://www.nmdb.eu), OMNIWeb (http://omniweb.gsfc.nasa.gov), Table 1 The functions implemented in the models for Solar Cycles 21 and 23.

Summary
i) We have computed the new 2D time-dependent PTE for long-period changes. In this equation we have used the well-known physical parameters describing the variation of the IMF, B, the tilt angle of the HCS, δ, and the parameters γ and ν for Solar Cycles 21 and 23. ii) The rigidity spectrum exponent, γ , describing a rigidity dependence of the long-period changes of the CR flux, and the exponent ν of the PSD of the IMF turbulence have been used in PTE as proxies. iii) We have shown that, for the first analyzed period (1976( -1987  If material is not included in the article's Creative Commons licence and your intended use is not 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/.

A.1 Calculation of the Exponent γ
Neutron monitor data in NMDB (http://www01.nmdb.eu/) were used to calculate the exponent γ of the rigidity R of the spectrum of GCR intensity variations ( δI (R) I (R) = AR −γ ) in the period of 1960-2002Alania, Iskra, and Siluszyk 2008, 2009, 2010. This was based on the article by Alania, Iskra, and Siluszyk (2003) who calculated the exponent γ using thoroughly selected monthly average data of neutron monitors during the period 1968-2002, including four ascending and four descending phases of solar activity in the A > 0 and the A < 0 epochs. The criterion for data selection adopted in these studies was the continuous operation of neutron monitors with different cut-off rigidities throughout the analyzed period. The magnitude J k i of the monthly averaged GCR intensity variations from the ith neutron monitor was calculated as where N k is the running monthly average count rate (months k = 1, 2, 3, . . . , 12) and N 0 is the monthly average count rate for the year of maximum intensity (usually in the minimum epoch of solar activity). The count rate of the maximum intensity is accepted as the 100% level. The magnitude J k i of the monthly averaged GCR intensity variations measured by the ith neutron monitor with the geomagnetic cut-off rigidity R i and the average atmospheric depth h i is defined by Dorman (1975) as where (δI (R)/I (R)) k is the rigidity spectrum of the GCR intensity variations for the kth month, W i (R, h i ) is the coupling coefficient for the neutron component of GCR by Dorman (1975) and Yasue et al. (1982), R max is the upper limit in rigidity beyond which the GCR intensity variation vanishes. For a power-law rigidity spectrum (δI (R)/I (R)) k = AR −γ k one can write where A k i is the magnitude of the GCR intensity variations rescaled to the heliosphere (free space). From Equation A.2 we obtain The values of A k i should be the same (within the accuracy of the calculations) for any ith neutron monitor if the values of parameters γ k and R max are properly determined. A similarity on the values of A k i for various neutron monitors is an essential argument to affirm that the data from a particular neutron monitor and the methods of calculations for γ k are reliable. To find the temporal variations in the exponent γ k (months k = 1, 2, 3, . . . , 12), a minimization of the expression ϕ = n i A k i − A k 2 (A.4) (where A k = 1 n n i A k i and n is the number of neutron monitors) was provided by Siluszyk et al. (2005) and Alania, Iskra, and Siluszyk (2008, 2009, 2010. The values of the expression Rmax R i R −γ k W i (R, h i ) dR for ranges of R max (from 30 GV up to 200 GV with a step of 10 GV) and γ (from 0 to 2 with a step of 0.05) were calculated based on the method presented in Dorman (1975) and Yasue et al. (1982). The upper limit in rigidity, R max , is taken to be 100 GV. This assumption was regarded as reasonable for the 11-year variation of the GCR intensity by Gushchina et al. (2008). The minimization of (A.4) for the smoothed monthly means (with the interval of 13 months) of the GCR intensity variations was carried out with respect to γ k for the neutron monitors.

A.2 Calculation of the Exponent ν
We collected the data of the B y component of the IMF for the period of 1976-2009 from SPIDR (http://spidr.ngdc.noaa.gov). Then we calculated the power spectrum density (PSD) and its exponent ν on the frequency f by the method of Blackman and Tukey (Lyons, 1996).
First of all, we calculated the autocorrelation function; if {B yi } (i = 1, . . . , N) is a series of daily values of B y of the IMF, the autocorrelation function R r (r = 0, 1, . . . , m) is given as where N = 365 and m = 182. We calculate the discrete Fourier transform of the autocorrelation function, PSD k, (k = 0, 1, . . . , m) from the formula PSD k = 2 t R 0 + 2 m−1 r=1 R r cos πkr m + R m cos πk m .
Here PSD k corresponds to the power spectral density for the frequency f k = k 2m t [Hz]. The series R r and {B yi } are given with a time interval of t = 1 day. The whole frequency range (0, 1 2 t ) is divided into m subintervals of length 1 2m t [Hz]. Therefore, the frequency data points are 0, 1 2m t , 2 2m t , . . . , m 2m t . We approximated the dependence of PSD k on the frequency f k by a power-law function as PSD = Pf −ν . The values of P and ν are derived by using the least-squares method.