© ESO 2022

1. Introduction

Recurrent high speed streams (HSSs) in the solar wind that originate from coronal holes cause various important effects in the heliosphere (e.g., Nolte et al. 1976; Schwenn 2006; Tsurutani et al. 2006; Vršnak et al. 2007; Rotter et al. 2012; Hofmeister et al. 2018, and references therein); for a review, see Cranmer (2009). The interaction of HSSs with the preceding slow solar wind creates the so-called corotating interaction regions (CIRs; e.g., Gosling & Pizzo 1999; Schwenn 2006; Tsurutani et al. 2006, and references therein); for a review, see Richardson (2018). These HSS and CIR structures cause recurrent moderate geomagnetic storms and substorms, where the Disturbance Storm-Time index (Dst) rarely goes below −100 nT (e.g., Tsurutani et al. 2006; Vršnak et al. 2007, 2017; Zhang et al. 2007; Verbanac et al. 2011a,b, 2013), as well as temporal depressions in the galactic cosmic ray (GCR) count rate, which are generally known as recurrent Forbush decreases (RFDs; for a review, see, e.g., Iucci et al. 1979; Kunow et al. 1995; Richardson 2004, 2018). Generally, Forbush decreases (FDs; Forbush 1937) are temporal, well-defined decreases in the GCR count rate that are caused by interplanetary coronal mass ejections (ICMEs) and CIR structures (e.g., Lockwood 1971; Cane 2000; Richardson 2004, 2018; Dumbović et al. 2011, 2012, and references therein); for reviews, see Cane (2000) and Richardson (2004). We note here that some authors exclusively use the term Forbush decrease for GCR depletions caused by ICMEs, and those that are caused by CIRs, they call “recurrent GCR depressions”. RFDs are usually relatively weak, and their amplitudes are considerably lower than those of “standard” FDs caused by ICMEs (e.g., Richardson 2004; Čalogović et al. 2009; Dumbović et al. 2011, 2012; Badruddin & Kumar 2016; Melkumyan et al. 2019, and references therein).

In recent years, a number of studies related to RFDs were published. For example, Kumar (2014), Badruddin & Kumar (2016), Ghanbari et al. (2019), and Guo et al. (2021) studied the behavior of RFDs by applying the superposed epoch analysis approach. We also refer to the observation-based studies by Modzelewska et al. (2020) and Shen et al. (2020), as well as to the numerical studies by Alania et al. (2011), Guo & Florinski (2016), Kopp et al. (2017), and Luo et al. (2020).

The interaction of the high-speed stream with the preceding slow solar wind creates the CIR structure by compressing the plasma in the region in which they interact (Gosling & Pizzo 1999; Jian et al. 2006; Richardson 2018). The compression region includes parts of the slow solar wind and parts of the HSS, divided by the so-called stream interface. The compressed interaction region is characterized by a continuously increasing flow speed, consistent with the continuity equation. The compression region is followed by the rarefaction region, in which the flow velocity gradually decreases, which is again consistent with the continuity equation. The whole CIR and HSS structure (especially the compression region) is characterized by significantly increased magnetic field fluctuations that are mainly related to the Alfvénic-wave perturbations.

Transport of GCRs in the heliosphere is usually described in terms of the so-called Parker GCR transport equation (Parker 1965). It considers the GCR particle density at a given time and position, which is determined by the effect of the convection and diffusion (i.e., two terms from the standard Fokker-Planck equation) and by a term added by Parker, describing the effect of adiabatic cooling of particles within an expanding volume (or heating in the case of compression). Later on, the Parker GCR transport equation was employed in many studies concerning various aspects of the heliospheric propagation and behavior of GCRs (e.g., Gleeson & Axford 1968; Perko & Fisk 1983; Jokipii et al. 1993; Fisk et al. 1998; Potgieter 1998; Wibberenz et al. 1998; Wibberenz & Cane 2000; Morales-Olivares & Caballero-Lopez 2010; Gieseler & Heber 2016; Zimbardo et al. 2017; Richardson 2018, and very many references therein).

The diffusion rate is determined by the diffusion coefficient, which depends on the mean-free path of particles, where the scattering is dominantly caused by the magnetic field fluctuations (Parker 1965). While the diffusion coefficients are calculated, the corresponding total fluctuating magnetic energy is usually computed using the root mean square variation in the vector of the interplanetary magnetic field (e.g., Zank et al. 1998). The particle drift effect can also be included in the transport equation (e.g., Fisk et al. 1998, and references therein), as well as the effect of mirroring on steep magnetic field gradients and shocks (e.g., Kirin et al. 2020, and references therein). The diffusion regulates the inward transport of particles, whereas the solar wind flow regulates the outward transport of GCRs. The adiabatic term describes changes in the GCR energy spectrum.

This paper aims to explain the behavior of recurrent Forbush decreases caused by CIR and HSS structures in the solar wind at heliospheric distances about 1 au, based on a convection-diffusion approach. Global 3D time-dependent simulations obtained by solving the GCR transport equation have shown that the smaller diffusion coefficient in the regions of stronger magnetic field related to CIRs results in a rapid local decrease in cosmic ray flux (Kota & Jokipii 1991). We also note that using the convection-diffusion approach to explain the CIR-related GCR modulation is not new (e.g., Chih & Lee 1986; Wibberenz et al. 1998). Richardson et al. (1996) considered the steady-state diffusion-convection model, which includes adiabatic deceleration and longitudinal variations in the solar wind speed. In this model, the effect of the convection due to the outward flow of the solar wind is locally balanced by the inward diffusion of GCRs. The model that we present in this paper is to a certain degree similar to the more sophisticated 3D numerical model by Alania et al. (2011). They regarded a steady-state heliosphere with CIR represented as the heliolongitudinal variation in solar wind speed and applied a full Parker transport equation. They found that the variation in GCR intensity is inversely correlated with the modulation parameter, which is proportional to the product of the solar wind velocity and the strength of the interplanetary magnetic field. Moreover, their calculations have shown that the heliospheric current sheet (HCS) structure in the model does not significantly affect the results. A similar result was also found by Guo & Florinski (2016) using the 3D coupled MHD-CR transport model. They found that GCR variations depend on the ratio of diffusion coefficients in the fast and slow solar winds. Finally, we compare our results to the most recent dedicated 3D coupled MHD-CR transport models by Kopp et al. (2017) and Luo et al. (2020). Their simulations show that enhanced convection, reduced diffusion, and drifts due the magnetic field enhancement may contribute to the recurrent FD. However, the simulations cannot distinguish the relative contribution of different mechanisms.

In this respect we note that an application of the modeling results to observations and the corresponding quantitative analysis is still lacking. Hereafter, the GCR count rate measured in situ by spacecraft (s/c) is denoted as CR, whereas for the analog in the model option, the symbol Φ is applied.

For the purpose of our analysis, we use some observational examples related to a long-living CIR structure that was recorded over 27 solar rotations in the period from June 2007 to May 2009. It was analyzed in detail by Dumbović et al. (2022, hereafter, Paper I; see also Kühl et al. 2013). In Sect. 2 the observational aspect is described in detail, and in Sect. 3 the analytical model is proposed. A comparison of model results with observations is given in Sects. 4 and 5. The outcome of the analysis is discussed in Sect. 6, and conclusions are drawn in Sect. 7.

2. CIR structure and the related GCR flux modulation

In Fig. 1, two examples of well-defined CIRs are presented. The events were recorded in the periods between the day of the year DOY = 184 and 191, 2007 (Carrington rotation 2059; Fig. 1a), and DOY = 58 and 68, 2008 (Carrington rotation 2067; Fig. 1b). Both CIR episodes were related to a long-living equatorial coronal hole that produced a persistent CIR and HSS that lasted for 27 rotations. The presented examples were recorded during its second and eleventh rotation (hereafter, rot-2 and rot-11, respectively.

Fig. 1. 1 au measurements of CIRs recorded from DOY 183 to 191, 2007 (rot-2; left column) and from DOY 58 to 68, 2008 (rot-11; right column). Panels a and e: flow speed, V. Panels b and f: magnetic field strength, B (red), and fluctuations (root mean square), δB (blue). Panels c and g: proton number density, N (red), and temperature, T (blue). Panels d and h: EPHIN cosmic ray count-rate change, CR, expressed in percent of the pre-event count rate. The dashed rectangle in both panels marks the period of the flow speed increase.

In Figs. 1a–c and e–g, the 1 h resolution data from the OMNI database¹ (King & Papitashvili 2005) are displayed, and Figs. 1d and h show the 1 h resolution data from the Electron Proton Helium Instrument (EPHIN; Müller-Mellin et al. 1995) on board the Solar and Heliospheric Observatory (SOHO), measured by its F-detector, which records hydrogen and helium nuclei above 53 MeV (Müller-Mellin et al. 1995). However, it is important to note that although the lower cutoff of SOHO-EPHIN-F is about 53 MeV, the response function of SOHO-EPHIN is such that it increases with energy, that is, it is less sensitive to low-energy particles (Dumbović et al. 2020).

The rot-2 event began on DOY = 184.3, 2007, starting from a quiet-Sun wind speed of V ≈ 340 km s⁻¹ (Fig. 1a). The HSS speed attained a maximum of ≈610 km s⁻¹ at DOY = 185.3, after which it gradually decreased to only ≈305 km s⁻¹ on DOY = 190.8. The phase of the rising speed is marked in panels a–d of Fig. 1 by the dashed blue rectangle. The period of increasing flow speed corresponds to the enhanced magnetic field, B, and plasma density, N (red graphs in Figs. 1b and c, respectively) caused by the interaction of the HSS with the ambient slow solar wind. This period is also characterized by enhanced magnetic field fluctuations, δB, and rising temperature, T (blue graphs in Figs. 1b and c, respectively). The change of the EPHIN cosmic ray count rate, CR, expressed in percentiles of the pre-event value, is presented in Fig. 1d. The CR(t) graph very closely corresponds to the inverted V(t) graph (for more examples, see Paper I). The CR data show a clear “over-recovery” phenomenon (e.g., Jämsén et al. 2007; Dumbović et al. 2012) after the main CR dip.

The rot-11 CIR began on DOY = 58.5, 2008, starting from a speed of V ≈ 350 km s⁻¹ (Fig. 1e). The HSS speed attained a maximum of ≈780 km s⁻¹ at DOY = 61.2, after which it gradually decreased to ≈320 km s⁻¹ on DOY = 68. The phase of the rising speed is marked in panels e–h of Fig. 1 by the dashed blue rectangle. The period of increasing flow speed again corresponds to the enhanced magnetic field, B, and plasma density, N (red graphs in Figs. 1f and g, respectively), as well as to the enhanced magnetic field fluctuations, δB, and rising temperature, T (blue graphs in Figs. 1f and g, respectively). The change in the EPHIN cosmic ray count rate is presented in Fig. 1h. Again, the CR(t) graph very closely corresponds to the inverted V(t) graph. The data presented in Figs. 1f–g show a distinct secondary peak in the B(t), δB(t), and N(t) graphs between DOY ≈ 65 and 66, corresponding to a weak increase in flow speed that is recognizable in the V(t) graph (Fig. 1e), and is also notable as a weak new rise in the T(t) graph (Fig. 1g). At the same time, CR(t) shows a distinct secondary decrease, with a minimum at DOY = 65.7 (Fig. 1h).

In Fig. 2a, a simplified generic temporal behavior of CIR structure and the associated RFD is sketched schematically (based on the detailed statistical study presented in Paper I). The abscissa shows the mean value of the relative time differences, Δt, of the onsets, maxima, and ends of V, N, B, and CR relative to the onset of the flow-speed increase. It is expressed in days (thus, Δt = 0 for the V onset). The y-axis in the y > 0 range represents the mean CIR-disturbance amplitudes, normalized with respect to the pre-CIR values, based on the mean values of the amplitudes of N, B, and V (to highlight the relative amplitude of V, it is multiplied by a factor of two). The [x_i, y_i] points obtained for the mentioned key parameters, that is, onsets, maxima, and so on, are finally connected by straight lines (this is just a schematic drawing). The graph includes the rarefaction region at the rear of HSS. In the y < 0 range, the mean RFD amplitude is shown as the percentage of the pre-CIR level.

Fig. 2. “Generic” temporal behavior of the flow speed, V (blue), magnetic field, B (red), density, N, (green), and the cosmic ray count-rate change expressed in percent, CR (black), based on the study of a sample of CIRs presented in Paper I. Panel a: simplified, generic, presentation based on the mean values of the amplitudes and timing of the B, N, and CR onsets, peaks, and ends, relative to the onset of V (to see better the relative amplitude of V, it is multiplied by a factor of two). The peak times are marked by vertical dashed lines. Panel b: SEA results; the variations relative to the pre-CIR values are normalized to the peak values, ΔV/Vmax, ΔB/Bmax, ΔN/Nmax, and ΔCR/CRmax. The time Δt = 0 is set at the onset of the V increase.

The plasma density starts to increase before the onset of the flow-speed rise on average, whereas the magnetic field in this representation starts to increase approximately simultaneously with the flow-speed increase. After the compression period, a rarefaction phase occurs, characterized by weakening of the magnetic field and plasma density decrease (e.g., Hundhausen 1973; Gosling 1996; Gosling & Pizzo 1999; Richardson 2018, and references therein). The CR decrease occurs during the period of the speed increase, corresponding to the phase of the enhanced magnetic field. Before the CR decrease, as is quite often observed, a weak enhancement in the count rate is seen in the period when only the density is increased, like in the super-epoch analysis (SEA) representation (Fig. 2b). The CR recovery phase closely follows the decrease in the flow speed. Regarding the SEA approach, it should be mentioned that Badruddin & Kumar (2016) analysed the differences between characteristics of CIRs that developed a frontal shock and those without the shock. We emphasize that in our sample, CIRs did not show a shock signatures, and our SEA results are quite similar to those presented by Badruddin & Kumar (2016).

Figure 2b shows an analogous, but more detailed, statistics-based generic time profile of V, B, N, and CR, where the displayed results are based on the SEA, performed in Paper I, covering the complete set of 27 rotations. Unlike the schematic presentation shown in Fig. 2a, the graph displayed in Fig. 2b reveals that not only does N starts to increase before V, but so does B (but within the standard deviation of Δt). Moreover, the graph V(t) reveals a slight depletion before the main rise of V. The graph CR(t) shows the pre-CIR slight increase in CR in more detail (e.g., Lockwood 1971), which is hereafter denoted as “nose”. As we show below, the nose is most likely a consequence of the pre-event phase of the decreased V and the steep magnetic field gradient related to the compression caused by the interaction of the fast and slow solar wind (Kirin et al. 2020).

The characteristics of the CIR structure and the associated RFD behavior presented in Figs. 1 and 2 are the basis for constructing a model for the CIR-related RFDs. It is also used for a direct comparison of the observations and model results.

3. Model

3.1. Physical background

In the following, we present a simple model for CIR and HSS effects on the cosmic ray count rate in the heliocentric distance range of r = 1 au (the global situation is sketched and the coordinate system explained in Figs. 3a and b). The model describes a single rotation, that is, rotations are considered independent (differences from one to another rotation are due only to different quiet-Sun values of Φ₀, V₀, and B₀, and the temporal behavior of Φ(t), V(t), and B(t) across the CIR and HSS at 1 au).

Fig. 3. Panel a: sketch of the global situation (upper panel: rest frame; lower panel: frame of reference rotating with the Sun); for details, see the main text. Panel b: schematic presentation of the geometry and CIR characteristics (fixed reference frame), as applied in the model (x0 and xe represent the beginning and end of CIR and HSS at 1 au, respectively; for details, see the main text). Panel c: dependence of the radial CR-gradient, Gr, for three specific s/c positions (x1, x2, x3) as a function of heliocentric distance, r. The heliocentric distance of the CIR frontal edge for the three specific s/c positions, xi, are denoted r1, r2, r3. Panel d: corresponding radial distance dependence of the cosmic ray count rate, Φxi(r), between CIR and HSS outer boundary and s/c. The red line marks the behavior in a quiet-Sun situation at the CIR outer boundary. Panel e: sketch of the resulting Forbush decrease profile, Φ(t), as measured by the s/c.

The model is based on the Fokker-Planck diffusion equation,

$$\begin{array}{c}\hfill \frac{\partial \mathrm{\Phi }}{\partial t}+\mathrm{\nabla }·(\mathrm{\Phi }\mathit{V})-\mathrm{\nabla }·({\kappa }_{\mathit{ij}}\mathrm{\nabla }\mathrm{\Phi })=0,\end{array}$$(1)

which was considered by Parker (1965) as a staring point in developing a model for the heliospheric propagation of GCRs. In Eq. (1), Φ represents the GCR particle number density at a given time and position, V denotes the convective velocity, and κ_ij represents the diffusion coefficient tensor. Equation (1) does not consider the change in particle energy spectrum due to the adiabatic heating/cooling of particles (for a complete GCR transport equation, including the adiabatic effect, see Parker 1965). In other words, we consider in the following only the change in total particle number density, regardless of the shape of the energy spectrum.

We note that although the lower cutoff of SOHO-EPHIN is about 53 MeV, the response function of SOHO-EPHIN-F is such that it increases with energy, that is, it is less sensitive to low-energy particles. On the other hand, the spectrum of GCRs decreases with energy. As a result, the spectrum of particles detected by EPHIN therefore looks like a bell-like curve with a peak somewhere around 1 GeV (depending on IP conditions, see, e.g., Dumbović et al. 2020). Bearing in mind that at 1 au, the proton spectrum recorded by the EPHIN-F peaks somewhat below but close to 1 GeV, small changes in proton energy spectrum therefore do not significantly affect the spectral distribution in the energy range ≈1 GeV measured by s/c.

In the following, we consider a CIR structure in a stationary state, that is, we assume that its structure and characteristics do not change significantly over the time period required for its passage over the spacecraft. In this case, the time-derivative ∂Φ/∂t = 0 at a given position within the CIR and HSS, so that Eq. (1) can be written as

$$\begin{array}{c}\hfill \mathrm{\nabla }·(\mathrm{\Phi }\mathit{V})=\mathrm{\nabla }·({\kappa }_{\mathit{ij}}\mathrm{\nabla }\mathrm{\Phi }),\end{array}$$(2)

that is,

$$\begin{array}{c}\hfill \mathrm{\Phi }\mathit{V}={\kappa }_{\mathit{ij}}\mathrm{\nabla }\mathrm{\Phi },\end{array}$$(3)

which means that the effect of the convection due to the solar wind outward flow is locally balanced by the inward diffusion of GCRs. The CIR-related decrease in Φ is a result of the suppressed inward diffusion due to the enhanced magnetic field and its fluctuations, and it is also due to the enhanced outward convection. The diffusion coefficient parallel to the field lines κ_∥ is far larger at 1 au than the perpendicular one κ_⊥, that is, the diffusion perpendicular to the magnetic field is negligible compared to the diffusion throughout the field (e.g., Parker 1965; Fisk et al. 1998; Giacalone 1998; Richardson 2018). We emphasize that the value of κ_∥ is quite ambiguous as it cannot be determined directly from observations. Its estimate therefore depends on a number of assumptions. Consequently, quite different values are found in various studies; they vary from 10¹⁷ to > 10¹⁹ m² s⁻¹. For example, Parker (1965) estimated it to be 10¹⁷ − 10¹⁸ m² s⁻¹, ≈3 × 10¹⁹ m² s⁻¹ was reported by Giacalone (1998), 10¹⁸ m² s⁻¹ was adopted by Gleeson & Axford (1968) and Schwenn & Marsch (1991), and Wibberenz et al. (1998) considered values 1.4 − 2.8 × 10¹⁸ m² s⁻¹. In the following, we use as a starting point for the quiet-Sun solar wind the radial diffusion coefficient of κ_r = 1.25 × 10¹⁸ m² s⁻¹, which is obtained from the CR radial gradient at 1 au of G_r = 4.8%/au (Gleeson & Axford 1968) and for the solar wind speed of V = 400 km s⁻¹, and using κ_r = V/G_r (e.g., Gleeson & Axford 1968; Schwenn & Marsch 1991, see also Eq. (6)).

In constructing the model, we assumed that the diffusion coefficient is related to the magnetic field magnitude B, based on a number of previous studies (e.g., Potgieter 2013; Richardson 2018, and references therein). In particular, we considered the simplest option, κ ∝ 1/B (see, e.g., Potgieter 1998; Wibberenz et al. 1998; Wibberenz & Cane 2000). In addition, because the mean free path of particles depends on the magnetic field fluctuations (e.g., Wibberenz et al. 1998; Wibberenz & Cane 2000; Michałek & Ostrowski 2001; Shalchi & Schlickeiser 2004; Chhiber et al. 2017, and references therein), as does the diffusion coefficient, we also took the magnetic field fluctuations δB into account, considering the simplest possible option κ ∝ 1/δB. In this respect, we note that according to the model proposed by Qin (2002), which is based on standard quasi-linear theory, the parallel mean free path should be proportional to B^5/3/δB². As we show below, this is not consistent with observations, which show that the diffusion coefficient is roughly inversely proportional to both B and δB. Under these two approximations, the perturbed value of the diffusion coefficient could therefore be estimated as ${\kappa }_B^{\prime }={\kappa }_0$B₀/B, or equivalently, κ_δB′=κ₀δB₀/δB, where κ₀, B₀, and δB₀ represent the values in the unperturbed solar wind.

However, the relation is most likely oversimplified, that is, it is more complex than the simplest options we assumed (for other options, see, e.g., Wibberenz et al. 1998; Wibberenz & Cane 2000; Michałek & Ostrowski 2001; Shalchi & Schlickeiser 2004; Chhiber et al. 2017, and references therein). For example, the mean free path, thus also the diffusion coefficient, depends on the particle energy, meaning that from the observational point of view, it depends on the energy channel of measurements (which is not of concern for our study because we considered only one channel). The relation κ(δB) also depends on the spectrum of fluctuations. As a consequence, it probably depends on the heliocentric distance and on the solar cycle phase (in our case, this is not important as we considered a relatively short time period relative to the solar cycle duration, and we focused on heliocentric distances of about 1 au). To compensate for all these ambiguities, we introduced a free parameter k that is expected to be different from event to event due to differences in fluctuation time-spatial characteristics, and it would also depend on the energy channel of the measurements.

Thus, we introduce the scaling coefficients k_b and k_δb, which we incorporated in b^* = k_bb and δb^* = k_δbδb, where we used the abbreviations b = B/B₀ and δb = δB/δB₀. In other words, b^* and δb^* were used in the model calculations, where k_b and k_δb were determined by fitting the model results to the observations by minimizing the rms differences between the modeled and observed data. With this notation, the effective local radial gradient can be expressed as $G_r^{*}=V/{\kappa }^{*}$, where V is the local solar wind velocity, and κ^* is defined by the unperturbed value of κ₀ and the effect of the perturbed magnetic field, as determined by the parameters b^* or δb^*, that is, κ^* = κ₀/b^* or κ^* = κ₀/δb^*, which means G^* = (V/κ₀)b^* = G₀b^* and equivalently, G^* = (V/κ₀)δb^* = G₀δb^*.

We now consider as the first step of the following considerations the situation in the reference frame that rotates with the Sun, that is, the frame whitin which the CIR and HSS are at rest, whereas the s/c passes through it. In this reference frame, the solar wind flow velocity, which is radial in the rest frame of reference, becomes aligned with the magnetic field lines because the velocity now has an additional component, corresponding to the solar rotational velocity (see Fig. 3a, and, e.g., Figs. 1, 4, and 5 in Hundhausen 1972). The advantage of the frame rotating with the Sun is that the diffusion and the wind flow are coaligned, both following the magnetic field lines, so that the diffusion is antiparallel with the convective flow. Because κ_⊥ ≪ κ_∥ means that the distribution of Φ along a given field line, Φ(ξ), where ξ represents the curvilinear coordinate along the field line, is independent of the distribution of Φ along neighboring field lines. Thus, considering the situation along a particular field line, Eq. (3) can be written as

$$\begin{array}{c}\hfill \mathrm{\Phi }(\xi )V_{\xi }={\kappa }_{\Vert }\frac{\partial \mathrm{\Phi }(\xi )}{\partial \xi },\end{array}$$(4)

that is,

$$\begin{array}{c}\hfill \frac{\partial (ln\mathrm{\Phi }(\xi ))}{\partial \xi }=G_{\xi },\end{array}$$(5)

where we have introduced the abbreviation for the GCR gradient along a given field line,

$$\begin{array}{c}\hfill G_{\xi }=\frac{V_{\xi }}{{\kappa }_{\Vert }}·\end{array}$$(6)

Equation (6) defines in firm physical terms a local gradient of Φ(ξ), G_ξ, along a given field line. If G_ξ is approximately constant over a distance range from ξ₁ to ξ₂ along a particular field line, or in other words, if we take the mean value of G_ξ over this segment of the field line, ${\overline{G}}_{\xi }$, the solution of Eq. (5) reads

$$\begin{array}{c}\hfill \mathrm{\Phi }({\xi }_2)=\mathrm{\Phi }({\xi }_1)e^{{\overline{G}}_{\xi }({\xi }_2-{\xi }_1)}.\end{array}$$(7)

For the vicinity of 1 au, we can write

$$\begin{array}{c}\hfill \mathrm{\Phi }(\xi )={\mathrm{\Phi }}_{1\mathrm{au}}{e}^{{\overline{G}}_{\xi }\xi },\end{array}$$(8)

where Φ_1 au represents the value of Φ at s/c located at 1 au, and we set ξ = 0 at the intersection of a given field line with the 1 au s/c trajectory.

This means that when we have the 1 au s/c measurements of Φ_1au (or analogously, at any other heliocentric distance at which a particular s/c is located), as well as the 1 au measurements of V and B (providing an estimate of G_ξ), Eq. (8) provides an estimate of Φ at any position [x, r] within or out of CIR and HSS. Moreover, if measurements of Φ(t) and V from two s/c located at two different heliocentric distances, r₁ and r₂, were available, we could easily estimate the value of ${\overline{G}}_{\xi }$ within the corresponding distance range ξ₁ − ξ₂ along a particular field line, assuming the nominal shape of the corresponding Parker spiral. This would also provide an estimate of the mean value ${\overline{\kappa }}_{\Vert }$ within the distance range ξ₁ − ξ₂ by using Eq. (6). By also employing the measurements of B and δB, we can also find an empirical expression relating κ_∥ to B or δB.

Finally, we emphasize that Eq. (8) clearly shows that after the s/c exits the HSS, in an ideal case in which the flow speed and the magnetic field attain pre-CIR values, that is, the value of G_ξ returns to the initial state, the cosmic ray particle density should also return to the pre-event level. If the flow speed and the magnetic field values differ from those measured before the CIR, the post-event cosmic ray count rate level should be different from that recorded in the pre-event state. If values of V and B are lower than in the pre-event state, the so-called over-recovery would be observed, in which the cosmic ray count rate after the HSS is higher than prior to it. If values of V and B in the post-HSS phase remain higher than in the pre-event phase, the count-rate will stay lower than it was before the CIR.

3.2. “Operative” approximation

Although the consideration presented above describes the basic behavior of Φ throughout the CIR and HSS structure, it is inconvenient for quantitative analysis and comparison with the measured values of CR at 1 au because the values of Φ_1 au at different 1 au s/c positions that are basically governed by Eq. (5) cannot be related directly to the pre-event values Φ₀. This is primarily because in the absence of measurements by two (or more) s/c overtaken by the same CIR structure at different heliocentric distances, the value of Φ at a certain referent distance is not available (for uncertainties of Φ at large heliocentric distances, see, e.g., Fig. 3 in McKibben et al. 1982).

In the following, we therefore approach the problem in an alternative way that provides an approximative, but directly operative quantitative analysis of Φ_1 au(t). More precisely, the proposed procedure described below enables estimating the ratio of the value of Φ at 1 au along the x-axis (for the definition of the coordinate system we used, see Fig. 3; in Fig. 3b we show the situation in the fixed reference frame) relative to the pre-event value Φ₀, that is, the values of Φ_1 au(x)/Φ₀, which can then be readily converted into the ratio Φ_1 au(t)/Φ₀, as measured by s/c during the CIR and HSS passage over it. To do this, we converted the dependence Φ_x(ξ) into an equivalent that defines the dependence of Φ on the radial coordinate, r, at a given x-coordinate, Φ_x(r). In other words, the gradient G_ξ(x) as estimated for 1 au was converted into the radial gradient, ${\overline{G}}_r(x)$.

The procedure of obtaining the ratio Φ_1 au(x)/Φ₀ was started by calculating the values of G_ξ along the x-axis G_ξ(x, 1 au). Then, taking into account the nominal orientation of the local Parker spiral (geometrically, the Archimedian spiral), we calculated the component of the gradient G_ξ(x, 1 au) in r-direction to obtain local values, G_r(x, 1 au). The Archimedian spiral is defined in polar coordinates [r, θ] as

$$\begin{array}{c}\hfill r=\frac{v_r}{\omega }\theta ,\end{array}$$(9)

where v_r is a constant radial speed and ω is a constant rotational speed. The winding angle, representing the tilt angle of the field lines relative to the radial direction, reads

$$\begin{array}{c}\hfill tan\psi =\frac{B_{\varphi }}{B_r}=\frac{\omega r}{v_r}·\end{array}$$(10)

The winding angle relates G_ξ and G_r as

$$\begin{array}{c}\hfill G_r=G_{\xi }tan\psi .\end{array}$$(11)

This situation is sketched in Fig. 3b (fixed reference frame) for three different s/c positions, x_i, denoted as x₁, x₂, and x₃, where the corresponding values G_ξ1, G_ξ2, and G_ξ3 are also sketched together with the related components G_r1, G_r2, and G_r3. The CIR and HSS structure is shaded in gray, and the frontal edge of the CIR is marked by the bold red line. The magnetic field lines are sketched by the thin blue lines. The heliocentric distance of the CIR boundary at coordinates x₁, x₂, and x₃ are denoted as r₁, r₂, and r₃, respectively. Ahead of the CIR boundary lies a quiet-state slow solar wind (shaded light gray), characterized by the GCR radial gradient G_r0 and the field-aligned gradient G_ξ0. At the bottom of Fig. 3b, the typical behavior of the 1 au solar wind flow speed, V, the magnetic field magnitude, B, and the cosmic rate count rate, CR, are sketched schematically.

In the next step, we estimated for each s/c position x within CIR and HSS the values of G_r(x, r) along the r-direction within the CIR and HSS structure, that is, from r = 1 au up to the outer CIR boundary located at r = r_cir(x), meaning that we converted the measurements obtained by s/c at 1 au into what is expected in the entire CIR and HSS structure within the range from its 1 au onset to its 1 au end. To illustrate this mapping procedure, we considered a moment t = t_i after t = 0, set at the CIR onset as measured by the s/c. At this moment, the plasma element that was at 1 au at t = 0 now already lies at r_cir = 1 au + V_r0t_i, where V_r0 is the solar wind velocity measured at t = 0, so that r_cir represents the radial distance of the frontal CIR boundary at t = t_i. Analogously, we can track the trajectory of any plasma element that was at 1 au in the period from t = 0 to t = t_i. For example, the element that was at 1 au at t_j < t_i, showing the radial speed V_j, at t_i will be at the distance r_j(t_i) = 1 au + V_j (t_i − t_j). When we take into account that the solar wind speed around 1 au does not change with distance, and when we also assume that in the considered distance range the magnetic field does not change much, the physical conditions within a given plasma element remain mainly unchanged over its trajectory within the t_i − t_j interval. Thus, the values of G_r estimated from s/c measurements (as previously described) for any t_j in the period from t = 0 to t_i can be attributed to the corresponding distance r_j between r_cir(x_i) and 1 au. Based on this mapping of G_r(r) at each x_i, it is straightforward to estimate the mean value of G_r from 1 au to r_cir for a given s/c position x_i, that is, for the corresponding moment t_i.

As a digression, we note here that the increasing V_r with the increasing time t_i, that is, increasing distance x_i within the CIR and HSS structure, also means that faster plasma elements catch up the preceding slower elements. Thus, over a larger heliospheric distance range, the previous approximations cannot be applied because the plasma and magnetic field within a plasma element are compressed due to [∂V_r/∂r]_x < 0. However, the compression effect is partly compensated for by the 1/r² lateral expansion of the structure.

The steepening of the frontal CIR structure leads to one more important effect, and that is the forward-shock formation (Gosling 1996; Gosling & Pizzo 1999). As the spatial V_r(r) profile causes a gradual steepening of the frontal part of the CIR, it will at a certain distance r_s inevitably transform into a shock. Denoting the radial speed of the frontal CIR element as V_r0, the maximum 1 au HSS radial speed as $V_r^{\text{max}}$, their difference as $V_r^{\text{max}}-V_{r0}\equiv \mathrm{\Delta }V$, the time elapsed from the 1 au CIR onset to the moment at which the HSS attains $V_r^{\text{max}}$ as Δt_max, and the corresponding distance between the CIR outer boundary and 1 au as Δr, we can estimate the time needed for the shock formation as Δt = Δr/ΔV = V_r0Δt_max/ΔV.

We take as typical examples the measurements presented in Figs. 1a and e, where the frontal radial speed was V_r0 ≈ 350 km s⁻¹ and the maximum speeds within the HSS were $V_r^{\text{max}}\approx 600$ and 800 km s⁻¹, achieved within ≈1 and 3 d, respectively. The corresponding times needed for the full completion of the shock after the frontal CIR edge passed over the 1 au s/c are then equal to 1.4 and 2.3 d, respectively, and the corresponding heliocentric distances are ≈1.5 and 2.1 au. This estimate is broadly consistent with observations as well as with numerical results (e.g., Gosling 1996; Gosling & Pizzo 1999).

We return to the previous considerations of the radial GCR gradient. The obtained values of G_r(x_i, 1 au), whose behavior is schematically sketched in Fig. 3c for three different s/c positions, x₁, x₂, and x₃, provide the estimate of the cosmic ray count rate variation along r for a given s/c position x_i. Beyond the outer CIR boundary, located at r_cir(x_i) the value of G_r(x_i, r_cir(x_i)) is taken to be constant (e.g., McDonald et al. 1981), having the value G_0r (Fig. 3b), which can be chosen from some of the values reported in a number of previous studies, for instance, 4.8%/au, as reported by Gleeson & Axford (1968); 3.5%/au, as estimated by McDonald et al. (1981); 1−10%/au, and 3−7%/au over different phases of the solar cycle, as reported by Fisk et al. (1998) and Fujii & McDonald (1993), respectively; 2.7%/au, as found by de Simone et al. (2011) and Gieseler & Heber (2016); 2−6.6%/au, as determined for different energy ranges by Marquardt & Heber (2019). These values depend not only to the phase of the solar cycle and the heliospheric distances considered, but also on the proton energies taken into account, consistent with the predictions by Parker (1965).

This enables an estimate of the values of Φ(x, 1 au), starting from the quiet-state value of Φ(x, r_cir) at the outer boundary of CIR and HSS as defined by the quiet-state radial gradient G_0r. The behavior of Φ_xi(r) at a particular x-coordinate, x_i, is now described by the equation for the r-direction,

$$\begin{array}{c}\hfill \frac{1}{{\mathrm{\Phi }}_{\mathit{xi}}(r)}\frac{\partial {\mathrm{\Phi }}_{\mathit{xi}}(r)}{\partial r}=G_r(x_i,r).\end{array}$$(12)

Equation (12) defines the behavior not only at 1 au, but rather, at any heliocentric distance r within, as well as beyond 1 au, as given by the solution of Eq. (12),

$$\begin{array}{c}\hfill ln{\mathrm{\Phi }}_{\mathit{xi}}(r)={\displaystyle \int }{G}_r(x_i,r)\mathrm{d}r.\end{array}$$(13)

To obtain the values of Φ at 1 au for a given x_i coordinate, Eq. (13) has to be integrated step by step, starting from r_cir(x_i), where the cosmic ray count rate has the nominal quiet-state value for that heliospheric distance, Φ₀(r_cir(x_i)), down to r = 1 au to derive the needed value Φ_1au(x_i).

The heliocentric distance

$$\begin{array}{c}\hfill r_{\mathrm{cir}}(x_i)=1\mathrm{au}+\mathrm{\Delta }{r}_{\mathrm{cir}}(x_i)\end{array}$$(14)

is calculated by substituting into Eq. (14) the value Δr_cir(x_i) = V_0rΔt_i, where V_0r denotes the pre-event radial solar wind speed, and Δt_i represents the time measured from the onset of CIR at 1 au, t₀, that is, Δt_i = t_i − t₀.

The obtained values of G_r(x, r) within the CIR and HSS structure define the behavior of G_r(r) for any coordinate x_i. This is schematically sketched in Fig. 3c for three different x-coordinates, x₁, x₂, and x₃, corresponding to segments at the beginning portions of the CIR, the place where G_r(1 au) is maximum, and a portion within the rear parts of HSS, respectively. Thus, Eq. (13) can now be evaluated to obtain Φ_xi(r), and thus Φ_xi(1 au), and consequently, Φ_1 au(t_i). This is schematically depicted in Figs. 3d and e, respectively.

On the one hand, bearing in mind all the approximations that we used, and on the other hand, the inherent noise in measurements and the process of smoothing, we simplify below the calculation of Φ_xi(1 au), i.e., Φ_1 au(t_i), where instead of integrating Eq. (13) step by step along the r-direction at a given x-coordinate, x_i, we estimate the mean value of G_r, along r at x_i, to get ${\overline{G}}_r(x_i)$, and use it to obtain the value of Φ_1 au(x_i) as

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{1\mathrm{au}}(x_i)=\mathrm{\Phi }(r_{\mathrm{cir}}(x_i)){\mathrm{e}}^{-{\overline{G}}_r(x_i)\mathrm{\Delta }r(x_i)},\end{array}$$(15)

where Δr(x_i) = r_cir(x_i)−1au represents the radial extent of CIR and HSS beyond 1 au at a given s/c position x_i, as defined below Eq. (14).

Having the values of Φ_1 au(x_i), i.e., Φ_1 au(t_i), the values of Φ_1 au(t_i) can be straightforwardly converted into the change of Φ_1 au relative to the pre-CIR value Φ₀, expressed in percentiles,

$$\begin{array}{c}\hfill \mathrm{\Delta }{\mathrm{\Phi }}_{\%}(t)=100\frac{{\mathrm{\Phi }}_{1\mathrm{au}}(t)-{\mathrm{\Phi }}_0}{{\mathrm{\Phi }}_0}=100(\frac{{\mathrm{\Phi }}_{1\mathrm{au}}(t)}{{\mathrm{\Phi }}_0}-1),\end{array}$$(16)

which can be directly compared with the observations in the form as presented in Figs. 1d and h, or the statistics-based generic form shown in Figs. 2a and b.

Finally, we note that the described model is to a certain degree similar to the propagating diffusive barrier approach (PDB, Wibberenz et al. 1998; Wibberenz & Cane 2000), which describes the effects of merged interaction regions (MIRs) at large heliocentric distances. However, we note that in PDB, which are valid for large heliocentric distances at which the magnetic field lines are almost concentric circles, the radially outward convection is balanced by the diffusion perpendicular to the magnetic field.

4. Comparison with observations

In the following, the model proposed in the previous section is applied to the observations, using the results of the study presented in Paper I, where the characteristics of a long-living CIR that was recorded over 27 solar rotations in the period from June 2007 to May 2009, are analyzed in detail. First, we consider in Sect. 4.1 the generic scheme shown in Fig. 2a, as well as the results of the SEA displayed in Fig. 2b. Then, in Sect. 4.2 we consider two specific examples, recorded from DOY 183 to 191, 2007 (rot-2) and from DOY 58 to 68, 2008 (rot-11). They are presented in Fig. 1.

As specified in Sect. 3, the main model input parameters are the solar wind radial flow speed, V_r and the magnetic field strength, B, or its fluctuations, δB. We also note that we employed the previously mentioned free parameter k when we relate the magnetic field strength (or its fluctuations) to the radial gradient of the GCR count rate, G_r, because the dependence of the diffusion coefficient on the magnetic field strength is not well known and the background heliospheric GCR radial gradient, G_r0, as inferred in a number of previous studies, spans quite a broad range of values. The GCR gradient should be higher in stronger fields, that is, roughly G_r ∝ B due to κ_∥ ∝ 1/B (Potgieter 1998; Wibberenz et al. 1998; Wibberenz & Cane 2000). We emphasize that the free parameter k does not affect the overall behavior of the modeled cosmic ray count rate within CIR and HSS. It is only used to adjust the overall amplitude of the model outcome to the observed amplitudes of the cosmic ray count rate variations. In other words, when the value of k is determined from fitting the model results to the observations, we obtain a better insight into the relation between the diffusion coefficient and the magnetic field strength, or its fluctuations.

4.1. Generic SEA and “schematic” time profiles

In the left column of Fig. 4 we show a comparison of the outcome of SEA (see Fig. 2b), with the model result for scaling factors k_B = 0.8 and k_δB = 0.85, estimated by matching the model with the observations. The values of coefficients k were obtained by minimizing the rms differences between the modeled and observed data. The model output is based on a smoothed V_r(t) curve and smoothed curves b(t) = B(t)/B₀ and δb(t) = δB(t)/δB₀. The smoothing is based on segmental linear fitting (hereafter, SLF), where linear fitting is applied to each segment between two neighboring knees in the measured time curves of a given parameter (for details, see Appendix A). We chose the SLF option instead of directly comparing the 1 h observed data with the model results simply to reduce the noise in the observed-minus-model data, especially in the cases of B and δB, which are characterized by quite high noise. In Fig. 4c the SEA SLF-smoothed V_r(t) curve (dashed black line) is compared with the inverted curve of the SEA SLF-smoothed cosmic ray count rate variation, CR(t) (blue lines), to illustrate a close similarity of the V_r(t) and −CR(t) behavior. In Fig. 4e the modeled SEA-based Φ(t) curves are directly compared to the SLF-smoothed SEA CR(t) curve. The red rectangle emphasizes the period of rising solar wind speed, ∂V_r/∂t > 0.

Fig. 4. Comparison of the model results and observations: SEA, based on the study in Paper I, shown in Fig. 2a (left) and the schematic generic behavior, based on Fig. 2b (right). Panel a: magnetic field, B (gray), and its fluctuations δB (green), normalized with respect to the pre-event values, b = B/B0 and δb = δB/δB0, multiplied by the free parameter k (written in the inset), based on SEA. Panel b: magnetic field, B, normalized with respect to the pre-event values, b = B/B0, multiplied by the free parameter k (written in the inset), based on the schematic behavior. Panels c and d: solar-wind flow speed, Vr (dashed black line) and the inverted graph of the cosmic ray count rate, CR (blue line). Panels e and f: observed CR (blue line) compared to the model results (dashed gray line) based on the behavior of b; in panel e: model results for the δb option are presented as the dashed green line. The red rectangles mark the period of the rising solar wind speed.

In Figs. 4b, d, and f, a comparison of the generic schematic time profiles (Fig. 2a) and the model B option results for k = 1.5 is presented, analogously to that concerning the generic SEA-option shown in Figs. 4a, c and e, respectively. In both options, model results and measured data match well. Moreover, the nose in the CR time profile corresponds to a slight pre-CIR depletion in the velocity curve.

4.2. Two specific examples: rot-2 and rot-11

In Fig. 5 we present a comparison of the model results and observations for rot-2 (left column) and rot-11 (right column). These two examples are chosen due to their well-defined and “clean” CIR and HSS structure, that is, there was no effect of other transient disturbances, such as ICMEs or another CIR that overlapped those we studied. Furthermore, they are different to a certain degree, not only by the amplitudes and perturbation duration, but also because rot-2 shows the over-recovery of CR at the end of HSS (Fig. 5e, see also Fig. 1d), whereas rot-11 shows a distinct secondary peak in B and δB at the rear of HSS with a corresponding secondary decrease of CR (Fig. 5f, see also Fig. 1h).

Fig. 5. Comparison of the model results and observations for two CIR events. Left: rot-2 (DOY = 184 and 191, 2007). Right: rot-11 (DOY = 58–68, 2008). Panels a and b: magnetic field and its fluctuations normalized with respect to the pre-event values (b and δb, gray and green; respectively). Panels c and d: solar-wind flow speed (dashed black) and the inverted CR(t) graph (blue). Panels e and f: observed CR(t) values (blue) compared to the model results based on the behavior of b (dashed gray) and δb (dashed green).

Figures 5a and b show the SLF-smoothed behavior of the normalized magnetic field strength, b = B/B₀, and its fluctuations, δb = δB/δB₀, for rot-2 and rot-11, respectively. The displayed values are scaled by the factor k we used in the calculations (k_b = 0.65 and 0.9 for b, and k_δb = 0.3 and 0.75 for δb, respectively). In Figs. 5c and d, the SLF-smoothed behavior of V_r is drawn together with the SLF-smoothed inverted CR curve to emphasize the similarity of the two curves. In Figs. 5e and f, the observed behavior of CR (blue) is compared with the model results based on b(t) (dashed gray) and δb(t) (dashed green). The graphs are very similar to the modeled and observed behavior, including the over-recovery in the case of rot-2, as well as the secondary dip in rot-11.

5. Analysis of the model results

In Fig. 6 the model results (y-axes) are plotted versus the observed results (x-axes) for the rot-2, rot-11, SEA, and generic cases. The slopes of the linear least-squares fits, y = c₁x + c₀, are found to be c₁ ≈ 1 on average, and the corresponding y-axis intercept is c₀ ≈ 0 (a perfect match corresponds to c₁ = 1, c₀ = 0). The correlation coefficients are high; they range from cc = 0.78 to 0.94. The results displayed in Fig. 6 show that the correlation coefficients are better for the b-option than for the δb-option, which is expected because of the various effects that have to be taken into account (high noise in the δb-data, a strong dependence on the time resolution, the unknown spatial scales and timescales of the fluctuations that affect CR, etc.). However, we note that the clear correlation between the model results and the observations for the δb-option is consistent with the hypothesis that the magnetic filed fluctuations govern the GCR diffusion term of the GCR transport equation, as expected from the theoretical point of view (Parker 1965). We note that δB ∝ B (see Sect. 6), which basically indicates that the fluctuations are of comparable normalized values, meaning that the turbulence attains larger amplitudes in stronger fields. The relative values remain comparable, however, δB/B ≈ 0.1.

Fig. 6. Quantitative comparison of the model values (y-axes) with the observed values (x-axes). Panels a–c: data drawn in blue and red represent the model results based on the normalized magnetic field strength b and its normalized fluctuations δb, respectively. Panel d: outcome for generic schematic profile based on the b model option is displayed. In the insets, linear least-squares fits are written, together with the corresponding correlation coefficients, cc.

In Fig. 5e the CR over-recovery is evident in both the observed and modeled values for rot-2, which is also reflected in Fig. 6a (dots in the top right quadrant; x > 0, y > 0, respectively). Moreover, the model reproduces the secondary dip of CR in the case of rot-11 (Fig. 5f) very well.

In Fig. 7 we present the distributions of the observed-minus-calculated CR values, O–C, for rot-2, rot-11, and SEA for model results based on the b- and δb-option. In all cases, the O–C values of the CR instantaneous amplitudes (expressed in percent; see Fig. 5) peak in the range from −0.5 to +0.5%. The best results are obtained for the rot-2/b-option and SEA/b-option (Figs. 7b and e, respectively).

Fig. 7. Distributions of observed-minus-calculated (O − C) values: Panel a: rot-2 modeled using the normalized magnetic field values b. Panel b: rot-2 modeled using the normalized magnetic field fluctuations δb. Panels c and d: the same as in panels a and b, but for rot-11. Panels e and f: the same as in panels a and b, but now for the SEA results.

In Table 1 we summarize the basic results of our comparison of the model results and observations. In Col. 1 all the considered options are specified. In Cols. 2 and 3 the average and median of the observed-minus-calculated values are shown. In Cols. 4–6 the basic parameters of the linear least-squares fit are presented (the slope, c₁, the y-axis intercept, c₀, and the correlation coefficient, cc, respectively). In the bottom row, the mean values for each column are displayed. The average and median values show a good overall match between the calculated and the observed values (the mean values presented in the last row are 0.05 ± 0.20% and 0.10 ± 0.22%, respectively). Generally, the slopes of the calculated-versus-observed linear least-squares correlations are grouped around c₁ = 1, with the mean value 0.99 ± 0.16. The y-axis intercept, c₀, ranges from −0.44 to +0.46, with the mean value of −0.04 ± 0.36, that is, c₀ ≈ 0. The correlation coefficients range from cc = 0.78 to 0.94, with a mean value of 0.85 ± 0.06.

6. Discussion

First of all, we emphasize that our results for the overall generic data, as well as for the two specific examples that were analyzed in more detail, are related to only one long-lasting CIR. Thus, the results we obtained should not be a priori generalized to all CIR events. This can be seen already from the fact that the scaling factor k is quite different in the case of rot-2 and rot-11, as well as in comparison with the SEA and schematic generic profiles. Still, the value of k for the b-options is approximately the value of k ≈ 1 on average (from the four inferred values of k, see Table 1, the mean value is ${\overline{k}}_b=0.96\pm 0.37$). The values of k for the δb-option seem to be somewhat lower (k_δb = 0.30, 0.75, and 0.85 for rot-2, rot-11, and SEA, respectively), with a mean value of $\overline{k}=0.63\pm 0.29$, but they are still within the same order of magnitude.

The dependence δB(B) for rot-2 and rot-11, presented in Fig. 8, shows that δB and B are correlated (the correlation coefficients are cc = 0.59 and 0.78, respectively). Although the number of data points is relatively small (n = 15 and 27, respectively), the F-test confidence is high (P = 98.5% and P ≫ 99%, respectively). When the data for both data samples are put together, the correlation coefficient adds up to cc = 0.70, with an F-test confidence P ≫ 99%.

Fig. 8. Dependence of the level of the magnetic field fluctuations, δB, on the related magnetic field strength, B, for rot-2 (red) and rot-11 (blue), together with the corresponding linear least-squares fit. The fit embracing all data points is shown by the dashed black line.

The linear least-squares fit for rot-2 results in δB = (0.06 ± 0.02) B + (0.13 ± 0.12), for rot-11 in δB = (0.13 ± 0.02) B − (0.21 ± 0.11), and for the joined data, the relation becomes δB = (0.098 ± 0.015) B − (0.06 ± 0.08). The slopes of c₁ ≈ 0.1 and the y-axis intercepts of c₀ ≈ 0 confirm the close correlation of these two magnetic field parameters, where the field fluctuations δB achieve about 10% of the instantaneous magnetic field strength B on average. The tight correlation of δB and B explains why the model b-option and δb-option both work well and why the scaling factor k in both options has a quite similar value.

At this point, it is appropriate to discuss the effect of the magnetic field fluctuations δB on the GCR particle propagation. The OMNI-reported fluctuations δB represent the root mean square of the average magnetic field magnitude within the time interval defined by the time resolution², in this case, 1 h. The measurements show that during the increased-B phase, the level of the fluctuations most often varyies from 0.2 to 2 nT, for example. The typical duration of these episodes of increased fluctuations is δt = 1 − 5 h. The length-scale equivalent, δx, can be calculated as δx = v_rotδt, where v_rot represents the velocity of the CIR rotation at 1 au. In this way, we obtain δx ≈ 0.01 − 0.05 au, that is, δx ≈ 1.5 − 7 × 10⁶ km, and this is somewhat larger than a typical 1 GeV proton gyroradius in the 5−10 nT environment. This means that these fluctuation clumps taken as an entity should not affect the propagation of these particles too much (Parker 1965), and thus should not affect the diffusion coefficient significantly. The considered fluctuation clumps are a result of many fluctuations on still smaller spatial scales that would also affect the propagation of 1 GeV protons, that is, the related diffusion coefficient and thus the GCR count rate, and this would need an additional, more careful comparison of the CR behavior with δB on different timescales.

Finally, we comment on the pre-CIR density increase and the associated V depletion (Figs. 1a and c, respectively; see also Fig. 2b). It can be speculated that a possible explanation might be related to the plasma pressure gradient within the frontal part of the CIR. This pressure gradient should cause a component of the flow speed along the field lines out of the frontal boundary of CIR that is compressed and thus carries additional plasma ahead of the CIR, that is, causes a weak pre-CIR density increase. As this flow component has a radial component that opposes the overall pre-CIR solar wind flow, this effect should also lead to the formation of the depletion in the radial speed time profile. In the situation shown in Fig. 1e, where there is no significant pre-CIR V depletion, there is no pre-CIR density increase either.

This scenario might be also related to the formation of the pre-CIR nose in the CR time profile. Namely, Fig. 2b shows that the nose occurs during the period of V depletion (see also Fig. 1d). According to the previously proposed model, a lower V, that is, a reduced convective effect, should lead to an increase in CR. The situation presented in Figs. 1e–h is somewhat different because the nose is still observed in CR, although there is no pre-CIR V depletion. On the other hand, inspecting the B and δB graphs in Fig. 1f, we find that in the period of the nose-related CR increase, there is a depletion in both B and δB, and according to the model, this should enable a more efficient diffusion and thus the CR increase.

Because the proposed explanation for this weak increase in GCR flux before the main FD decrease is quite speculative and cannot be fully confirmed by our analysis, we should also consider another possible explanation that is related to the mirror effect at a steep gradient of the magnetic field (including the situation in which a forward shock has already formed) in the region in which the fast solar wind compresses the plasma in the interaction with slow solar wind. In this respect, we note that in the situation in which cylindrical geometry cannot be applied and the adiabatic approximation is not valid, the mirror effect is quite complex and cannot easily be included in the hydrodynamical model such as employed here, that is, the numerical approach has to be applied (see Kirin et al. 2020). The results presented by Kirin et al. (2020) show that a certain population of particles tends to stay bound to the shock for a certain period of time, indicating that this effect may explain the weak increase in GCR flux before the main FD decrease. Related to the compression in the frontal parts of CIRs, we also emphasize the so-called snowplow effect (e.g., Guo et al. 2021), which is observed as the plasma density increase that starts even before the solar wind speed starts to rise, and which may also cause the GCR pre-CIR increase (Guo et al. 2021).

Finally, we note that although our model is based on 2D approach, it uses a realistic CIR structure in the equatorial plane as a structural background for a quite simple analytical model for the GCR modulation, which is based on the simplified Parker transport equation. Our model is similar to the steady-state diffusion-convection model by Richardson et al. (1996), which in addition considers the adiabatic deceleration and longitudinal variations in the solar wind speed. In both models, the effect of the convection due to the solar wind outward flow is locally balanced by the inward diffusion of GCRs. Richardson et al. (1996) concluded that the main cause of the CIR-related GCR depression is the increase in solar wind speed, that is, the enhanced outward-convection effects, whereas the enhanced magnetic field in the frontal parts of CIR as well as the turbulence-related drifts play a less important role. In contrast, we have found that the magnetic field and turbulence enhancement play an important role in the GCR modulation. The results of our model are on the other hand very similar to the much more sophisticated 3D numerical model by Alania et al. (2011). They considered a steady-state heliosphere in which the CIR is represented as heliolongitudinal variation of solar wind speed and applied the full Parker transport equation. They found that the variation in GCR intensity is inversely correlated with the modulation parameter, which is proportional to the product of the solar wind velocity V and the strength of the interplanetary magnetic field B. This is exactly the case in our model. Moreover, their calculations have shown that changing the HCS structure in the model does not significantly affect the results (for a detailed study of the effect of HCSs on GCRs, see, e.g., Thomas et al. 2014).

Results similar to those obtained by Alania et al. (2011) were also found by Guo & Florinski (2016) using the 3D coupled MHD-CR transport model. They also found that that GCR variations depend on the ratio of the diffusion coefficients in the fast and slow solar winds, as is the case in our model. Finally, we compared our results to the most recent dedicated 3D coupled MHD-CR transport models by Kopp et al. (2017) and Luo et al. (2020). Their simulations show that enhanced convection, reduced diffusion, and drifts due the magnetic field enhancement may contribute to the recurrent FD. However, the simulations cannot distinguish the relative contribution of different mechanisms. Our results indicate that convection and diffusion are enough to describe the basic properties of recurrent FDs. On the other hand, we do not exclude the possibility that drifts are needed to explain the fine structures within the recurrent FD. We also note that in all of the models above, including our model, the transport parameters need to be tuned to successfully reproduce the observations. In the quantitative aspect, our simple analytical model therefore performs equally well as the much more sophisticated and physically more realistic models.

7. Conclusion

Our analysis explains the physical background of RFDs caused by CIR and HSS structures in terms of the convection-diffusion approach. It is demonstrated that the enhanced convection-effect caused by the increased velocity of HSS and the reduced diffusion-effect caused by the enhanced magnetic field fluctuations within the CIR and HSS structure are sufficient to explain the observations. This approach was verified by applying a relatively simple analytical model to the observations.

The main guiding observational fact that leads in this direction is the close similarity of the V(t) curve and the inverted CR(t) curve (Figs. 4 and 5; for more examples, see Paper I). However, the enhanced convection effect cannot explain the CR(t) behavior without an additional ingredient. After we inspected the in situ solar wind data, we realized that this has to be the increase in magnetic field strength B within the CIR structure, or a physically more appropriate option, the enhanced magnetic field fluctuations δB that reduce the CR diffusion.

The combined effect of the increased solar wind speed and the enhanced magnetic field fluctuations practically entirely explain the overall RFD behavior in the 1 au range, related to the CIR and HSS structure at these distances. The two effects taken together can also explain the phenomenon of the over-recovery that is sometimes observed after the RFD events. It is likely that the over-recovery occurs (based on the two employed examples) in the events in which the solar wind speed and the magnetic field fluctuations are both lower than in the pre-event state.

Acknowledgments

This work has been supported by Croatian Science Foundation under the project 7549 “Millimeter and submillimeter observations of the solar chromosphere with ALMA”. B.V. and M.D. also acknowledge support by the Croatian Science Foundation under the project IP-2020-02-9893 (ICOHOSS). B.H. acknowledges the financial support via projects HE 3279/15-1 funded by the Deutsche Forschungsgemeinschaft (DFG).

References

Appendix A: Segmental linear fitting

In Fig A.1 the SLF procedure is illustrated. It is applied to the EPHIN GCR count rate, solar wind velocity time profile, and interplanetary magnetic field time profile. The procedure is as follows. First, the moments in a given curve are identified, where the slope of the curve (neglecting the noise) significantly changes (nodes). Then, the parts of the curve between nodes are fit by a linear fit (red lines).

Fig. A.1. Three examples of the SLF (rot-11). Panel a): EPHIN GCR count rate. Panel b): solar wind velocity time profile. Panel c): interplanetary magnetic field time profile.

All Tables

All Figures

Fig. 1. 1 au measurements of CIRs recorded from DOY 183 to 191, 2007 (rot-2; left column) and from DOY 58 to 68, 2008 (rot-11; right column). Panels a and e: flow speed, V. Panels b and f: magnetic field strength, B (red), and fluctuations (root mean square), δB (blue). Panels c and g: proton number density, N (red), and temperature, T (blue). Panels d and h: EPHIN cosmic ray count-rate change, CR, expressed in percent of the pre-event count rate. The dashed rectangle in both panels marks the period of the flow speed increase.

In the text

Fig. 2. “Generic” temporal behavior of the flow speed, V (blue), magnetic field, B (red), density, N, (green), and the cosmic ray count-rate change expressed in percent, CR (black), based on the study of a sample of CIRs presented in Paper I. Panel a: simplified, generic, presentation based on the mean values of the amplitudes and timing of the B, N, and CR onsets, peaks, and ends, relative to the onset of V (to see better the relative amplitude of V, it is multiplied by a factor of two). The peak times are marked by vertical dashed lines. Panel b: SEA results; the variations relative to the pre-CIR values are normalized to the peak values, ΔV/Vmax, ΔB/Bmax, ΔN/Nmax, and ΔCR/CRmax. The time Δt = 0 is set at the onset of the V increase.

In the text

Fig. 3. Panel a: sketch of the global situation (upper panel: rest frame; lower panel: frame of reference rotating with the Sun); for details, see the main text. Panel b: schematic presentation of the geometry and CIR characteristics (fixed reference frame), as applied in the model (x0 and xe represent the beginning and end of CIR and HSS at 1 au, respectively; for details, see the main text). Panel c: dependence of the radial CR-gradient, Gr, for three specific s/c positions (x1, x2, x3) as a function of heliocentric distance, r. The heliocentric distance of the CIR frontal edge for the three specific s/c positions, xi, are denoted r1, r2, r3. Panel d: corresponding radial distance dependence of the cosmic ray count rate, Φxi(r), between CIR and HSS outer boundary and s/c. The red line marks the behavior in a quiet-Sun situation at the CIR outer boundary. Panel e: sketch of the resulting Forbush decrease profile, Φ(t), as measured by the s/c.

In the text

Fig. 4. Comparison of the model results and observations: SEA, based on the study in Paper I, shown in Fig. 2a (left) and the schematic generic behavior, based on Fig. 2b (right). Panel a: magnetic field, B (gray), and its fluctuations δB (green), normalized with respect to the pre-event values, b = B/B0 and δb = δB/δB0, multiplied by the free parameter k (written in the inset), based on SEA. Panel b: magnetic field, B, normalized with respect to the pre-event values, b = B/B0, multiplied by the free parameter k (written in the inset), based on the schematic behavior. Panels c and d: solar-wind flow speed, Vr (dashed black line) and the inverted graph of the cosmic ray count rate, CR (blue line). Panels e and f: observed CR (blue line) compared to the model results (dashed gray line) based on the behavior of b; in panel e: model results for the δb option are presented as the dashed green line. The red rectangles mark the period of the rising solar wind speed.

In the text

Fig. 5. Comparison of the model results and observations for two CIR events. Left: rot-2 (DOY = 184 and 191, 2007). Right: rot-11 (DOY = 58–68, 2008). Panels a and b: magnetic field and its fluctuations normalized with respect to the pre-event values (b and δb, gray and green; respectively). Panels c and d: solar-wind flow speed (dashed black) and the inverted CR(t) graph (blue). Panels e and f: observed CR(t) values (blue) compared to the model results based on the behavior of b (dashed gray) and δb (dashed green).

In the text

Fig. 6. Quantitative comparison of the model values (y-axes) with the observed values (x-axes). Panels a–c: data drawn in blue and red represent the model results based on the normalized magnetic field strength b and its normalized fluctuations δb, respectively. Panel d: outcome for generic schematic profile based on the b model option is displayed. In the insets, linear least-squares fits are written, together with the corresponding correlation coefficients, cc.

In the text

	Fig. 7. Distributions of observed-minus-calculated (O − C) values: Panel a: rot-2 modeled using the normalized magnetic field values b. Panel b: rot-2 modeled using the normalized magnetic field fluctuations δb. Panels c and d: the same as in panels a and b, but for rot-11. Panels e and f: the same as in panels a and b, but now for the SEA results.
In the text

	Fig. 8. Dependence of the level of the magnetic field fluctuations, δB, on the related magnetic field strength, B, for rot-2 (red) and rot-11 (blue), together with the corresponding linear least-squares fit. The fit embracing all data points is shown by the dashed black line.
In the text

	Fig. A.1. Three examples of the SLF (rot-11). Panel a): EPHIN GCR count rate. Panel b): solar wind velocity time profile. Panel c): interplanetary magnetic field time profile.
In the text