Guo2014 Cosmic Ray Modulation

Journal of Geophysical Research: Space Physics RESEARCH ARTICLE 10.1002/2013JA019546 Key Points: • Corotating interactio...

1 downloads 144 Views 3MB Size
Journal of Geophysical Research: Space Physics RESEARCH ARTICLE 10.1002/2013JA019546 Key Points: • Corotating interaction regions are simulated during the recent solar minimum • The 27 day periodic variation of cosmic rays are reproduced near Earth • Mechanisms of cosmic ray depression are discussed

Correspondence to: X. Guo, [email protected]

Citation: Guo, X., and V. Florinski (2014), Corotating interaction regions and the 27 day variation of galactic cosmic rays intensity at 1 AU during the cycle 23/24 solar minimum, J. Geophys. Res. Space Physics, 119, 2411–2429, doi:10.1002/2013JA019546.

Received 17 OCT 2013 Accepted 6 APR 2014 Accepted article online 8 APR 2014 Published online 28 APR 2014

Corotating interaction regions and the 27 day variation of galactic cosmic rays intensity at 1 AU during the cycle 23/24 solar minimum X. Guo1 and V. Florinski1,2 1 Center for Space Plasma and Aeronomic Research, University of Alabama, Huntsville, Alabama, USA, 2 Department of

Space Science, University of Alabama, Huntsville, Alabama, USA

Abstract We investigate the formation and evolution of corotating interaction regions (CIRs) in the solar wind and their effects on galactic cosmic rays (GCR) during the recent solar cycle 23/24 solar minimum. The output from a three-dimensional MHD model serves as background for kinetic time-dependent simulations of GCR transport based on the Parker equation. The results show that the CIR forward/reverse shock pairs or compression/rarefaction regions play important roles in the transport of GCR particles and directly control the observed 27 day periodic intensity variations. We find that stream interfaces (SIs) in CIRs and the heliospheric current sheet (HCS) are both closely associated with the GCR depression onset, in agreement with the observations at 1 AU. The HCS is more important when its tilt angle becomes small during the declining phase of the solar minimum, while the passages of SIs control the onset of GCR depressions for larger HCS tilt angles. The mechanism of GCR intensity variation near 1 AU can be explained through an interplay between the effects of particle drift and diffusion. The simulated plasma background and GCR intensity are compared with the observations from spacecraft and a neutron monitor on the ground, to find good qualitative agreement. Evidently, CIRs had a substantial modulational effect on GCR during the recent solar minimum. 1. Introduction Galactic cosmic rays (GCR) are high-energy particles (mainly protons) originating from outside the heliosphere but generally from within our galaxy. Most GCR are thought to be accelerated by shocks associated with supernova explosions that occur approximately every 50 years in our galaxy [Ackermann et al., 2013]. As the GCR travel into the heliosphere, they become deflected by the interplanetary spiral magnetic field and, for lower energies, are prevented from reaching the inner solar system. Gradient and curvature drift is an important mechanism assisting GCR entry and exit from the heliosphere [Jokipii and Kopriva, 1979]. From the global point of view, positively charged particles drift into the inner heliosphere along the heliospheric current sheet (HCS) at low latitude and leave through the two polar regions during qA < 0 (negative) solar minima. The cycle is reversed for positive qA. This drift pattern alternation plays a crucial role in the large-scale 11 year sunspot cycle and the 22 year magnetic cycle variations of GCR intensities [Kota and Jokipii, 1983]. Cross-field diffusion becomes more prominent when GCR approach the inner heliosphere within a few AU from the Sun [Kota and Jokipii, 1995]. At the same time, GCR locally gain or lose energies traveling through compression and rarefaction regions within the heliosphere, including the termination shock [Langner and Potgieter, 2004], the transient features in the heliosheath [Luo et al., 2011], and corotating interaction regions or CIRs [Kota and Jokipii, 1991]. Theoretical investigations of GCR modulation using global MHD steady state models of the heliosphere have been conducted by several groups [Ball et al., 2005; Florinski and Pogorelov, 2009; Strauss et al., 2013]. These models were based on a flat HCS concept, which is easier to implement in a time-independent simulation. The actual HCS is wavy, and the changes in its tilt angle control the long-term cosmic ray modulation to a large degree. Its most prominent effect is the alternation between peaked and plateaued cosmic ray intensity-time profiles that occurs every 11 years [Jokipii and Thomas, 1981; Potgieter and Moraal, 1985]. The wavy HCS is responsible for the sectored magnetic field structure in inner heliosheath, with profound implications for GCR transport in that region [Florinski et al., 2012]. In addition to these global trends, it was discovered that short time scale solar wind events, such as merged interaction regions (MIRs) could cause transient GCR depressions and even accumulate to produce the observed 11 year and 22 year cycle

GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

2411

Journal of Geophysical Research: Space Physics

10.1002/2013JA019546

variations [Burlaga et al., 1984; Potgieter and le Roux, 1994]. MIRs affect cosmic ray intensities mainly through reduced diffusion [Perko and Fisk, 1983], owing to the enhanced magnetic field strength and turbulence. CIRs are formed by the interaction between high- and slow-speed solar wind streams [Sarabhai, 1963], where the plasma is compressed and heated. Their structure is well developed by approximately 1 AU from the Sun [Siscoe, 1972; Gosling et al., 1972]. When the streams emanating from the Sun are approximately time stationary, these compression regions form spirals in the solar equatorial plane that corotate with the Sun. The leading edge of a CIR is a forward wave propagating into the slower stream ahead, while the tailing edge is a reverse wave propagating back into the trailing high-speed stream. By 2 AU, the forward and reverse waves typically steepen and become shocks [Hundhausen and Gosling, 1976]. These forward and reverse shock pairs propagate into undisturbed solar wind and eventually interact with the neighboring CIRs. At heliocentric distance greater than 8 AU, CIRs merge and combine to form the more complex structures called corotating merged interaction regions (CMIRs) [Burlaga et al., 1993]. CMIRs are reported to be dominant structures at low latitudes between 8 and 12 AU [Burlaga et al., 1984]. As the associated shocks weaken with the increasing heliocentric distance beyond 10–15 AU, CMIRs are replaced by corotating pressure enhancement regions [Gazis, 2000]. They are expected to evolve and interact with each other and form the “wave interaction zone” beyond some 25 AU [Burlaga et al., 1983]. However, recurrent forward shock-like structures were still observed at distances of 42–52 AU [Lazarus et al., 1999]. Observations also showed large amplitude fluctuations of solar wind parameters and magnetic field in the heliosheath region [Burlaga et al., 2010; Richardson and Wang, 2011] that could be considered CMIR remnants. Other kinds of merged interaction regions include local MIRs (LMIRs) and global MIRs (GMIRs) [Burlaga et al., 1993]. The former are produced from local transient solar wind streams on limited scales and are capable of causing short-term decreases in cosmic ray intensities. GMIRs have larger sizes and longer durations compared with LMIRs. They are known to produce large step decreases in cosmic ray intensities [Burlaga et al., 1993, 2003] and dominate GCR modulation during times of high solar activity, whereas the HCS tilt angle effects control modulation during low solar activity periods [Potgieter and le Roux, 1994]. CMIRs are formed from merging of neighboring CIRs during propagation and are usually accompanied by rarefaction regions, which cause successive decreases and increases in cosmic ray intensities. The net modulation effect of CMIRs over long timescales is expected to be small. Although transient CIRs and CMIRs affect short-term and small-scale variations of cosmic ray intensities, based on multiple observations [Richardson, 2004], it was also shown, using computer simulations [Kota and Jokipii, 1991; McKibben et al., 1999] that they do not change the global cosmic ray distribution in the heliosphere significantly. Despite this, questions remain about the modulation efficiency of CIRs on short timescales [e.g., Heber et al., 1999]. We also need a better understanding of the mechanisms of short-term periodic GCR variations present in observations. The 27 day variations in GCR intensities have been observed for many decades since the first announcement by Forbush [1938]. It was soon realized that the conditions in the interplanetary medium, caused by solar activity, were responsible for the periodic variation (see, e.g., the review by Simpson [1998]). Different mechanisms were proposed to cause the onset of GCR depression, including solar wind speed increases at stream edges, magnetic sector boundaries, magnetic field enhancements, and stream interfaces, or SIs (see the review by Richardson [2004]). It is currently thought that interfaces between fast and slow solar wind streams and the leading edges of CIRs are responsible for the depression onset [Richardson et al., 1996]. The magnetic field does not appear to order the GCR structures, and in some cases the multiple crossings of HCS are not accompanied by GCR depression. Clearly, the relationship between GCR intensity depressions and the CIR phenomena discussed above needs further investigation. The recent solar minimum provided an unprecedented opportunity to study these phenomena under relatively steady conditions. The recurrent variation of GCR intensities during this period have been extensively measured using ground-based neutron monitors [Modzelewska and Alania, 2013] and space probes [Leske et al., 2013]. For this work we have conducted a series of numerical experiments to elucidate the relationship between CIRs and GCR fluxes. CIR formation and evolution was previously studied using MHD models [Pizzo, 1994; Riley et al., 2001; Usmanov and Goldstein, 2006; Borovikov et al., 2012; Provornikova et al., 2012]. Here we use similar techniques to obtain the global heliospheric CIR structure for the region upstream of the termination shock, based on an idealized concept of the solar wind whose flow is time stationary in the corotating frame. The resulting plasma background is time dependent in the stationary frame (see below). Based on this variable plasma background, we trace the stochastic trajectories of large numbers of GCR pseudo-particles GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

2412

Journal of Geophysical Research: Space Physics

10.1002/2013JA019546

Figure 1. (a) Tilted-dipole plasma flow geometry used in the CIR simulation. (b) Color contour of the plasma velocity at the inner boundary r = 0.5 AU.

and compute the modulated GCR spectra in the heliosphere. We also perform a quantitative comparison between the results of our numerical simulations and the observational data.

2. Numerical Models 2.1. MHD Modeling of CIR Evolution Following on our previous CIR research [Florinski et al., 2013], we used the tilted-dipole flow geometry of Pizzo [1991] at the inner boundary to produce CIRs in a global MHD simulation. The problem geometry is illustrated in Figure 1. In this figure, 0xyz is the stationary coordinate frame, and 0x ′ y′ z′ is the corotating frame. Assumed to be perpendicular to the ecliptic plane xy, z is the solar rotation axis, while z′ points along the magnetic axis in the Northern Hemisphere; 𝛾 is the angle between the two axes (i.e., the dipole tilt angle). Fixed at 30◦ in our simulations, 𝛽 is the latitude of the fast-slow transition centerline in the corotating coordinate system. The point S(𝜃, 𝜑) lies on the line through the center of the fast-slow transition. Its coordinates in the 0xyz frame may be obtained from x 2 + (x tan 𝜑 cos 𝛾 + z sin 𝛾)2 = cos2 𝛽,

z cos 𝛾 − x tan 𝜑 sin 𝛾 = sin 𝛽,

(1)

where tan 𝜑 = y∕x and z = sin 𝜃 . One can then derive an expression for 𝜃(𝜑) and compute the polar angle as 𝜑 varies with solar rotation. The latitude of the HCS in the 0xyz frame can be obtained by taking 𝛽 = 0 in equation (1). Similar to the method used in Borovikov et al. [2012], a cubic polynomial f (𝜃 ′ ) is introduced to smoothen the transition between the fast and slow streams, 3

2

f (𝜃 ′ ) = a1 𝜃 ′ + a2 𝜃 ′ + a3 𝜃 ′ + a4 ,

(2)

where f is some physical quantity, such as density or pressure, 𝜃 ′ ∈ [𝛽 − Δ𝜃 ′ ∕2, 𝛽 + Δ𝜃 ′ ∕2] is the angle defined in the corotating coordinate system 0x ′ y′ z′ , and Δ𝜃 ′ is the latitude of the transition, set to 40◦ in our model in order to fit the observations. The four coefficients a1 , a2 , a3 , and a4 are defined according to the following boundary conditions: f (𝛽 − Δ𝜃 ′ ∕2) = fslow , f ′ (𝛽 − Δ𝜃 ′ ∕2) = 0, f (𝛽 + Δ𝜃 ′ ∕2) = ffast , and f ′ (𝛽 + Δ𝜃 ′ ∕2) = 0, where ffast and fslow correspond to the values in the fast and slow solar wind, respectively. Note that the above expressions are written for the Northern Hemisphere; in the Southern Hemisphere, ffast and fslow should be interchanged. The standard Parker solar wind solution was used at the inner simulation boundary. The magnetic field is given by Br = B0 (r0 ∕r)2 , B𝜃 = 0, B𝜑 = −(Ωr∕ur )Br sin 𝜃 , where the suffix “0” denotes 1 AU values. To model the cycle 23/24 solar minimum, we used solar wind parameters derived from observations, as the boundary conditions for the code. The properties of the fast solar wind were extracted from Ulysses data [Ebert et al., 2009], while the slow wind properties were derived by averaging the Wind and ACE observations over the duration of the solar minimum [Jian et al., 2011]. Because recurrent transient events are often observed at GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

2413

Journal of Geophysical Research: Space Physics

10.1002/2013JA019546

Table 1. Solar Wind Condition at 1 AU Adopted for the 23/24 Solar Cycle Minimum and the Derived Boundary Condition at 0.5 AU Parameter

SW(Fast, 1 AU)

SW(Slow, 1 AU)

SW(Fast, 0.5 AU)

SW(Slow, 0.5 AU)

n, cm−3 |𝐮|, km s−1 T, K |𝐁r |, nT

1.61 745 219,000 2.4

5.50 334 64,100 2.4

5.99 728 528,400 9.6

20.84 323 154,700 9.6

1 AU as a result of solar rotation, the average measured solar wind speed is actually higher than the slow solar wind speed. To correct for this, we took daily averages of ACE plasma velocity values and removed the points lying above the solar minimum average of 388 km/s reported in Jian et al. [2011]. The data were then averaged over the period 2007.0–2010.0 to obtain the velocity of the slow solar wind. These values are then rescaled to the inner boundary using the standard Parker solution for the solar wind [Parker, 1958]. The resulting solar wind parameters at 1 AU and the corresponding boundary conditions at 0.5 AU are listed in Table 1. The solar dipole tilt axis 𝛾 was obtained from observations (http://wso.stanford.edu). We used a constant 25 day (sidereal) solar rotation period. The computational domain in the model was restricted to the interplanetary space and supersonic flow. We chose the inner boundary at 0.5 AU, well beyond the critical point (usually under 5 solar radii). The outer boundary was located at 80 AU, which is inside the measured distance to the termination shock, namely 94 AU for Voyager 1 [Stone et al., 2005] and 84 AU for Voyager 2 [Stone et al., 2009]. The MHD simulations were performed on a geodesic grid (see Florinski et al. [2013] for more details about this model). In the radial direction, 512 spherical layers were used; a level 6 geodesic grid with 40,962 hexagonal cells covered the surface of the sphere. The MHD code was run for about 1 year simulation time to establish a periodic CIR structure over the entire simulation domain. Because of the periodic inner boundary condition, plasma variations at each point in space have the periodicity of the solar rotation (25 days). A single snapshot from the simulation is then saved and used as background for the cosmic ray transport simulation. The time dependence is reintroduced by rotating this solution about the solar polar axis. 2.2. Stochastic Transport Model and GCR Modulation Our GCR transport model is based on Parker’s equation [Park, 1965] ( ) 𝜕f 𝜕f 𝜕f 𝜕 ∇ ⋅ 𝐮 𝜕f + (ui + vd,i ) = 0, − 𝜅ij − 𝜕t 𝜕xi 𝜕xi 𝜕xj 3 𝜕 lnp

(3)

where the pitch angle-averaged phase-space density f (𝐫, p) is a function of spatial position 𝐫 and the magnitude of momentum p. The plasma convection velocity is 𝐮, and the drift velocity 𝐯d includes the gradient and curvature drifts for an isotropic distribution of particles, ( ) pcw 𝐁 ∇× 2 𝐯d = (4) 3e B where w is the velocity of the particle and e is its charge. The diffusion tensor in equation (3) is 𝜅ij = 𝜅⟂ 𝛿ij + (𝜅∥ − 𝜅⟂ )bi bj , where 𝐛 = 𝐁∕B. The parallel and perpendicular diffusion components 𝜅∥,⟂ are functions of the particle’s momentum and the turbulent fluctuation energy. Diffusion is closely related to the magnetic fluctuations in the solar wind. We assume that the turbulence properties of the heliospheric plasma can be described with two basic parameters, the mean square fluctuating magnetic field intensity ⟨𝛿B2 ⟩, and the “bendover length” lb that marks the transition between the (flat) energy-containing range and the power law inertial region of the turbulent power spectrum [Zank et al., 1998]. The bendover length approximately corresponds to the characteristic size of large-scale eddies in the flow [Zank et al., 1996]. Similar to the previous treatments [le Roux and Fichtner, 1999; Florinski and Pogorelov, 2009], we assumed that 𝜌 ⟨𝛿B2 ⟩ = ⟨𝛿B2 ⟩0 (5) 𝜌0 at all latitudes, where ⟨𝛿B2 ⟩0 = 3.0 × 10−10 G2 is the turbulent field variance at 1 AU in the ecliptic plane and 𝜌0 is a characteristic solar wind density, taken to be 5 cm−3 . This scaling gives a reasonable agreement with the results of the turbulent transport models [Breech et al., 2008]. GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

2414

Journal of Geophysical Research: Space Physics

10.1002/2013JA019546

The parallel diffusion coefficient may be computed from the standard quasi-linear theory for magnetostatic slab turbulence [e.g., Dröge, 2000], 𝜅∥ =

wB2 rg2 4

1

∫0

𝜇(1 − 𝜇 2 )d𝜇 , Psl (1∕𝜇rg )

(6)

where rg is the gyroradius, 𝜇 is the cosine of the particle pitch angle, and Psl is the turbulent power spectral density in the slab component (the component with the wave vector parallel to the mean magnetic field). The latter is computed from the assumed power spectral density profile. For the combined energy-inertial power spectrum the expression for the parallel diffusion coefficient is [Zank et al., 1998; Giacalone and Jokipii, 1999] [ ( )5∕3 ] 72 Ωlb 3v 3 3𝜋 1+ 𝜅∥ = csc (7) , 5 7 v 20Ω2 lb A2sl where lb = 1.2618lc , lc = 0.03AU is the correlation length, Ω is the particle gyrofrequency, and A2sl = ⟨𝛿B2sl ⟩∕B2 is the amplitude of the slab turbulence. The intensity of the slab component is assumed to be 20% of the total ⟨𝛿B2 ⟩, the two-dimensional component making up the remainder [Bieber et al., 1996]. It is to be understood that the scalar quantity |B| in equations (4), (7), and the expression for Asl (but not in 𝐛) should be interpreted as (B2 + ⟨𝛿B2 ⟩)1∕2 , which ensures that the turbulent ratio never exceeds one. Cross-field diffusion depends on the presence of transverse fluctuations because in a slab turbulent field particles retrace their paths due to scattering, leading to a subdiffusive overall transport [Kota and Jokipii, 2000]. Diffusion is recovered if transverse structure is added to the turbulent field [Qin et al., 2002]. Here we use the simple model of le Roux et al. [1999], where 𝜅⟂ = aA22D 𝜅∥ ,

(8)

where a is a phenomenological constant, set to 0.02 in these simulations, and A2D = ⟨𝛿B22D ⟩∕B2 is the amplitude of 2-D turbulence. Considered without the numerical factor a, equation (8) described the situation where particles propagate diffusively along the magnetic field lines that have diffusive transverse displacement relative to the mean magnetic field 𝐁 and where the particle’s correlation length is equal to the correlation length of the dominant 2-D turbulence [le Roux and Webb, 2007]. The empirical numerical factor is needed to make the results consistent with the measured amount of heliospheric modulation [Florinski and Pogorelov, 2009]. To solve the transport equation (3), we follow the Monte Carlo method for stochastic processes [Jokipii and Levy, 1977; Zhang, 1999]. The stochastic (from the diffusion term) phase-space trajectories of the Parker equation are integrated backward in time until external boundary is reached. The equations actually being solved are the characteristics of equation (3), ( ) 𝜕𝜅ij dxi = − ui + vd,i − (9) dt + Aij dWj (t) 𝜕xj ∇⋅𝐮 dt (10) 3 where Aik Akj = 2𝜅ij and Wi (t) is a multidimensional Wiener (random walk) process. We start from a location at (𝐫, p) and integrate a large number (N) of trajectories until they reach the outer boundary. The particle distribution function can be obtained by averaging over the exit points, d ln p =

f (𝐫, p) =

N 1∑ f (p ) N n=1 ob ob,n

(11)

where fob is the cosmic ray spectrum assumed at the outer boundary and pob,n is the momentum of the nth at the exit point. Because the outer boundary in these simulations is not the heliopause, fob is not the interstellar GCR spectrum but a fit to Voyager observations near the termination shock (see Figure 7 below). GCR transport model is implemented on the same grid as the MHD simulation, which greatly facilitates the calculation of the transport coefficients in equations (9) and (10) from the plasma data. Because the CIR structure rotates with the Sun, the transport coefficients are transformed to the stationary frame along the trajectory at every step of the Monte Carlo simulation. GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

2415

Journal of Geophysical Research: Space Physics

10.1002/2013JA019546

Figure 2. Simulated (a) plasma velocity and (b) magnetic field strength in the meridional plane, showing the CIR structure in the solar wind, for period 1. The wavy heliospheric current sheet is showing with a red line.

3. Simulation Results 3.1. The CIR Structure in the Heliosphere We focus on the time interval between January 2007 and December 2009 as representative of the solar minimum between solar cycles 23 and 24. This interval corresponds to Bartels rotation periods 2367 through 2407. It is further subdivided into two periods: period 1 from fractional year 2007.0 to 2008.5 and period 2 from 2008.5 to 2010.0. The reason for this distinction is the decrease in the HCS tilt angle during the second half of the solar minimum. The tilt angles we used are based on the classic PFSS (Potential Field Source Surface Models) model, which are larger, by about 10◦ , than the angles obtained from the radial model. The PFSS model uses a line of sight boundary condition at the photosphere and includes a significant polar field correction [Altschuler and Newkirk, 1969]. The radial model uses a radial boundary condition at the photosphere and requires no polar field correction. The results from PFSS model are better at reproducing the periodic variations of plasma variables near Earth, because the very small radial model tilt angles of ∼9◦ produce no CIRs at low latitudes. Because the emphasis here is on 27 day periodic variations, we use long-term tilt angle averages of 32.5◦ for period 1 and 22.0◦ for period 2, obtained by averaging the values reported by the Wilcox Solar Observatory. Figure 2 shows color contour plots of simulated velocity (Figure 2a) and magnetic field (Figure 2b) magnitudes in the meridional plane for period 1. As seen in the schematic diagram in Figure 1a, the fast-slow transition region rotates about the z axis, when observed from a fixed point in space. The fast and slow winds alternate with the solar sidereal rotation period of 25 days (27 days if viewed from Earth). As a result of fast solar wind streams overtaking the preceding slow streams, the well-known CIR structure appears at low and middle latitudes, clearly visible in the figure. The latitudinal extension of CIRs has been thought to be a controlling factor for GCR shielding during the solar minimum[Rouillard and Lockwood, 2007]. Because during period 1 the HCS tilt angle (32.5◦ ) is larger than the latitude of the fast-slow transition region center (𝛽 = 30◦ in Figure 1), the high-speed solar wind can reach the equatorial plane at certain time intervals during the rotation period. The slow solar wind eventually disappears altogether at larger heliocentric distances. As a result, the shape of the HCS (shown with a red line) is significantly distorted at these distances. Figure 3 shows the same physical quantities as Figure 2 but for period 2. Here the tilt angle is smaller, and the slow streams dominate the low latitudes, at all distances. As expected, the wavy HCS appears much smoother, and its latitudinal extent is smaller than in period 1. Because the HCS is the principle GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

2416

Journal of Geophysical Research: Space Physics

10.1002/2013JA019546

Figure 3. Same as Figure 2 but for period 2.

mechanism of cosmic ray transport during the negative solar minima, we expect more GCR to reach the inner heliosphere during period 2 than during period 1. Magnetic field strength variations in CIRs could have a prominent effect on GCR transport. Magnetic field increases correspond to “magnetic barriers” that impede GCR propagation as a result of a local decrease in the diffusion coefficients. Conversely, a depression of the magnetic field strength will facilitate GCR transport, compared to the undisturbed regions. The magnetic field contour plots in Figures 2 and 3 show an alternating radial pattern of the magnetic field magnitude increases and decreases, implying a complicated GCR transport behavior. The red lines in these figures show the cross section of the HCS in the meridional plane. Its shape was obtained using a kinematic method: the HCS location is obtained by tracing the streamlines of plasma flow originating on the neutral line at the inner boundary, similar to the technique described in Czechowski et al. [2010].

Figure 4. (a) Logarithm of the thermal pressure in the meridional (xz) plane, and (b) the profiles of key plasma variables along the direction 25◦ north latitude in the RTN coordinate system. Four shock pairs are identified, each consisting of a forward shock (FS) and a reverse shock (RS), shown with blue dashed lines. The red dotted lines mark the SIs. CMIR1 and CMIR2 are two CMIRs located at middle and low latitudes, respectively.

GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

2417

Journal of Geophysical Research: Space Physics

10.1002/2013JA019546

Figure 5. Same as Figure 4 but for period 2.

Compared with the flat current sheet, the HCS structure shown in Figures 2 and 3 is more obstructive for GCR transport because the particles are required to follow a longer path [e.g., Jokipii and Thomas, 1981]. It was also suggested that diffusion may effectively disrupt the drift process and lead to diffusive short circuiting of HCS by the cosmic rays when the HCS tilt angle becomes very large [Strauss et al., 2012]. Because the interplanetary magnetic field was weaker during the cycle 23/24 solar minimum than during the preceding negative solar minima, and the tilt angle correspondingly larger [e.g., Jian et al., 2011], diffusive effects are expected to play a more important role in GCR transport in the most recent minimum.

V (km/s)

Lat (°)

Figures 4 and 5 show the simulated radial profiles of the solar wind parameters for period 1 and 2, respectively. Left panels show the thermal pressure distribution using a logarithmic color map, with a superimposed HCS shown in red. The well-developed shock pairs are labeled in these figures as forward shocks (FS) and reverse shocks (RS). The leading edge (FS) of a CIR propagates into the slower plasma ahead, while the trailing edge (RS) propagates back into the trailing high-speed plasma. CIR shock pairs are separated from each other by rarefaction (low pressure) regions at distances less than 20 AU. Beyond that, adjacent CIRs merge to form CMIRs. One such merging area, identified as a pressure enhancement in the figure, is labeled CMIR1 in the left panel. 50 Another CMIR (CMIR2 ) appears at 0 lower latitudes, inside of 10 AU. -50 800 700 600 500 400 10 8 6 4 2 40 20 0 -20 40 20 0 -20 -40 2007.0

2007.5

2008.0

2008.5

2009.0

Years

Figure 6. (top) The Ulysses’ latitude (blue) and heliocentric distance (red). The other panels compare Ulysses observation of solar wind velocity, number density, and the two magnetic field components (blue) during its third obit from January 2007 to June 2009, and our CIR simulation results (red) for the cycle 23/24 solar minimum. The dashed vertical line separates periods 1 and 2.

GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

Theoretically, stream interfaces are tangential discontinuities separating the plasmas originated from the slow and fast solar winds. They are distinguished as abrupt and simultaneous drops in density and rises in temperature within CIRs [Belcher and Davis, 1971; Gosling et al., 1978]. In our MHD simulation, SIs are identified as local vorticity maxima, as seen in the bottom panels of Figures 4 and 5. The vertical dotted blue lines are the CIR shock fronts, the arrows showing their directions of propagation. We can see that FS3 and RS4 merge at a radial distance around 26 AU. Other CIR characteristics can be observed from 2418

Journal of Geophysical Research: Space Physics

Spectrum Intensity (/sec-m2-Mev-sterad)

10.00

10.1002/2013JA019546

the plasma profiles. For example, VN changes in directions across the stream interface, corresponding to south-north flows.

In order to validate our CIR simulation results, we compared the time-dependent model output with 1.00 the daily average Ulysses observations from January 2007 to June 2009 (Ulysses data are not available after June 2009). As shown in Figure 6, during this time interval Ulysses traveled from high southern latitudes at about 0.10 2.5 AU to low northern latitudes of Local Intersteller Spectrum about 20◦ at 4.7 AU. The simulated Modulated spectrum at Voyager 2 Spectrum near Earth (2007.0-2008.5) plasma velocity, density, and two Spectrum near Earth (2008.5-2010.0) magnetic field components qualitatively agree with the observations at low and high latitudes except within 0.01 0.01 0.10 1.00 10.00 the fast-slow transition region, where Energy (GeV) the amplitude of the fluctuations are much larger in the simulations. Figure 7. Proton spectra during the recent solar minimum. The unmodulated LIS is shown with a green line. The black line is the modulated These results are similar to those spectrum assumed at the outer boundary of our simulation, which is conreported previously by Borovikov et sistent with Voyager 2 measurement at 85 AU (pink triangles). The red and al. [2012] and are representations of blue lines are simulation results for the two periods, respectively. The range the periodic crossings of CIRs along of PAMELA observations during the past solar minimum is shown with the Ulysses trajectory. Despite the vertical dotted lines. larger-than-observed variations, the averaged simulated data approximate the solar wind conditions quite well during the recent solar minimum. 3.2. GCR Transport in CIRs Because the external boundary of the simulation is inside the termination shock, we alter the local interstellar spectrum (LIS) [Webber and Higbie, 2010] to account for the modulation that occurs at larger distances. In order to fit Voyager measurement near the termination shock, we use the following polynomial approximation for the proton spectrum: ) ( J = 11.0 T −2.6 ∕ 1.0 + 3.80 T −1.22 + 1.3 T −2.8 + 0.0087 T −4.32 (12) in unit of (sec ⋅ m2 ⋅ MeV ⋅ sterad)−1 , where T is the proton energy in GeV. This spectrum, plotted in Figure 7 with a black line, provides a reasonable fit to the Voyager 2 spectrum at 85 AU [Webber et al., 2008] (triangles). Using J = fp2 , equation (13) yields the spectrum at the boundary fob in equation (11). The unmodulated LIS is shown with a green line in Figure 7 for comparison. The red line in Figure 7 shows the simulated proton spectrum at 1 AU in the equatorial plane during period 1 from 2007.0 to 2008.5, while the blue line shows the same for the period 2008.5 to 2010.0. These results are consistent with the GCR observations, shown by vertical dotted segments, which cover the intensity variations observed by PAMELA from October 2006 to December 2009 [Adriani et al., 2013]. As the HCS tilt angle decreased from 32.5◦ to 22.0◦ , the spectrum showed an increase during period 2. This trend is caused by the enhanced drift of GCR particles along the wavy HCS. In addition to the increased role of drifts along the HCS, the variation in diffusion coefficient may have also played a role in the intensity increase, which will be discussed below. Figure 8 compares the simulated (red) and observed (blue) periodic variations of GCR intensity and plasma properties near Earth for period 1. In this figure, the 2 h averaged GCR data were obtained by the NAIN neutron monitor (http://www.nmdb.edu), while the daily averaged plasma observations were taken from the OMNIWEB database (http://omniweb.gsfc.nasa.gov/ow.html). The red line in the top panel is simply the fixed HCS tilt angle of 32.5◦ , while the blue curve shows the measured tilt angles based on the WSO observations. GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

2419

Journal of Geophysical Research: Space Physics

10.1002/2013JA019546

Figure 8. The 27 day periodic variations of GCR (E = 400 MeV) intensity and solar wind plasma from observations (blue) and simulations (red) during period 1 with a HCS tilt angle of 32.5◦ . From top to bottom: HCS tilt angle, the OMNI heliographic latitude (yellow), GCR variation in percent, plasma number density, thermal temperature, solar wind speed, magnetic field magnitude, and the RTN longitude (blue) and the simulated magnetic polarity (red). The vertical dashed lines mark the crossings of the equatorial plane by the OMNI virtual spacecraft.

The observed tilt angle drops slightly below average in the middle of the period and increases to about 5◦ above the average for the remainder of period 1. The measured GCR intensity shows a decrease and a subsequent rise, in phase with the tilt angle changes. The heliographic latitude of Earth is plotted in the top panel as a measure of the distance of the OMNI virtual spacecraft from the solar equator shown in Figure 1. The magnetic longitude, plotted in Figure 8 (bottom), is the angle measured in the RT-plane from the radial direction, positive toward BT . Sharp changes in this angle are used to locate the HCS in the data, shown with blue circle signs. The simulated magnetic polarity (positive is outward, negative is inward) is shown in the same panel for comparisons. The observed solar wind plasma varies quasi-periodically with large amplitude, at a primary period of 27 days. GCR variations measured on the ground have a similar periodicity, which has been verified through the harmonic analysis during the recent solar minimum [Modzelewska and Alania, 2013]. The diurnal variation of GCR intensity caused by Earth’s rotation will not be discussed here. The percentage rate of GCR intensity variation is about 1–3% (from peak to valley) from the neutron monitor, with the maximum occurring at around 2008.0. For comparison, we simulated the time-dependent variation of 400 MeV protons during the study period. The simulated variation rates of GCR intensity have amplitudes of about 3%. Note that neutron monitors are sensitive to GCR penetrating the Earth’s atmosphere with energies roughly above several GeV [e.g. Alanko et al., 2003], which is much larger than the reference energy of 400 MeV used in the simulation. Thus, the simulated GCR depression sizes are not expected to be comparable with those from the neutron monitor. Typically, the average GCR intensity depression size observed by neutron monitor is about one third of the guard-rate depressions observed by spacecraft at heliocentric distances from 0.9 to 1 AU [Iucci et al., 1979]. Figure 8 does not show the spacecraft guard-rate depression data. The ACE/SIS GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

2420

Journal of Geophysical Research: Space Physics

10.1002/2013JA019546

Figure 9. Same as Figure 8 for a HCS tilt angle of 22.0◦ .

data reported previously showed about 5.3% depression for low energy (roughly 100–400 MeV) GCR in 2007–2008 [Leske et al., 2011]. This is about 2 times the 3% variation obtained in Figure 8. This discrepancy remains an outstanding issue that will be addressed in future work. The orbit of OMNI virtual spacecraft is in the ecliptic plane, which is inclined by 7.2◦ relative to the solar equatorial plane. The corresponding latitude of the spacecraft is shown with a yellow line in Figure 8 (top). The tilt between the two planes causes a gradual change of the amplitudes of solar wind plasma variations during the period. The vertical dashed lines show the time when the spacecraft was located in the equatorial plane. For example, the observations show the amplitudes are at their smallest at about the year 2007.43. The amplitudes gradually increase and attain their maxima when the spacecraft reaches its highest latitude at time 2007.69. This long-term latitude-dependent trend can be clearly identified for the plasma speed from observations and simulations and for the GCR intensity variations from the simulations. During a synodic rotation period of 27 days, there are two main depressions in the GCR intensity seen in the simulations. When the spacecraft reaches its highest heliographic latitude, the depressions have magnitudes of about 1% and 3%. In the equatorial plane, the two depressions have equal magnitudes of 2%. This intriguing variation of the GCR depressions is caused by the rotational asymmetry of CIRs at different latitudes. GCR observations do not clearly show this latitude-dependent property. Both the magnetic longitude (blue circle signs) and magnetic polarity (red square wave profile) jump at the current sheets, nearly simultaneously during most of the second half period. The simulated solar wind plasma and GCR periodic variations are consistent with observations during the first period from 2007.0 to 2008.5. Figure 9 shows a comparisons between the simulations and observations performed for period 2. The HCS tilt angle clearly shows an anticorrelation with the average GCR intensities in the top two panels. The amplitudes of the GCR intensity and solar wind plasma variations are smaller than in Figure 8 because of the smaller HCS tilt angle of 22.0◦ , as well as due to the decline in solar activity. For this case, the maximum GCR variation is about 1.5%, which is about half of what it was for period 1. The latitude-dependent variation GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

2421

Journal of Geophysical Research: Space Physics

10.1002/2013JA019546

Figure 10. A subset of data shown in Figure 8. Two adjacent stream interfaces (SI1 and SI2 ) and magnetic sector boundaries (SB1 , SB2 ) are shown with vertical dashed lines. The green curve is plasma vorticity from the simulations, used to detect stream interfaces.

of GCR intensity was also observed is these simulations. The slow solar wind dominates the low latitudes, covering the OMNI virtual trajectories. Similar to the comparison during period 1, the simulated solar wind plasma and GCR intensity variations are in general agreement with the observations. More refined results could be obtained from short-term analysis of the simulations and observations during these two periods. Figure 10 shows a subset of Figure 8, sampled from 2007.76 to 2007.90. To facilitate comparison between the observation and the simulation, 1 h resolution solar wind data from OMNI were used. The observed tilt angle was nearly constant, as shown in the top panel. The latitudes of the spacecraft varied between 6.9◦ and 2.8◦ north. The magnetic polarity from the simulation (bottom) is a good fit to the magnetic RTN longitude of the observation. The vertical dashed blue and red lines show the magnetic sector boundaries (SB), marking the HCS crossings. Stream interfaces (SIs) are identified as an abrupt drop in density with a rise in temperature [Gosling et al., 1978]. For the simulated data, extrema of vorticity (green curve) were used to locate the SIs, similar to the method used previously in Figures 4 and 5. From the observations, we can see that GCR intensity went down after the passages of the two stream interfaces (blue SI1 and SI2 ), while the two HCS passages (blue SB1 and SB2 ) apparently are not correlated with the GCR variations during this period, and no obviously related GCR intensity depression is observed. Similarly, the two simulated stream interfaces (red SI1 and SI2 ) affect the GCR intensity variation, and their passages are associated with the onset of the two deeper depressions of GCR. These results demonstrate that stream interfaces do play a more important role in the GCR depression onset, compared with the proximity to the HCS, during period 1 of the recent solar minimum. Two kinds of GCR intensity depressions could be identified in the data shown in Figure 10: the short term and the long term. Typically, the former last less than a day, while the latter has much longer durations. For example, the intensity dropped sharply within 1 day after the passages of blue SI1 and SI2 . Subsequently, the GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

2422

Journal of Geophysical Research: Space Physics

10.1002/2013JA019546

Figure 11. Same as Figure 10 for a HCS tilt angle of 22.0◦ .

decrease became moderate, lasting several days before the GCR intensity bottomed out. Note that because of computational constraints, the simulated daily GCR intensity has lower temporal resolution than the data. The short-term depressions will be the subject of future work. Moving now to period 2, Figure 11 is a zoomed in view of Figure 9, sampled from 2009.24 to 2009.38. Here the simulated GCR intensity variations are in good agreement with the observational result. The spacecraft was located at Southern Hemisphere, with its latitudes ranging from −6.8◦ to −2.7◦ . The simulated magnetic polarity does not fit the magnetic RTN longitude well during this period. In contrast to period 1 (cf. Figure 10), the roles of the HCS and stream interfaces are reversed in period 2 (from observations). The HCS is now more important for the GCR variations than stream interfaces. The observations show that the two marked passages of HCS (blue SB1 and SB2 ) are accompanied by visible depressions of GCR intensity, while the depressions associated with the stream interfaces (blue SI1 and SI2 ) are much less distinct. Our simulation results have similar properties, except the passage of the stream interface (red SI1 ) still apparently results in a two-stage depression of GCR intensity as described earlier with regards to in Figure 10. Since the HCS tilt angle varied from period 1 to period 2, we believe that the enhanced role of the HCS crossings is directly related to the decrease of the HCS tilt angle from 32.5◦ to 22.0◦ . Theoretically, GCR intensity should peak near the HCS due to the inward drift of ions during the A<0 solar cycle. The solar rotation variations in GCR intensities are expected from the resulting latitudinal intensity gradient combined with the cyclically changing latitudinal distance from the rotating tilted HCS [Kota and Jokipii, 1995]. In addition, GCR intensity is expected to decrease after the passage of a stream interface because of the decrease of the diffusion coefficients [Intriligator et al., 2001]. Both the passages of the HCS and stream interfaces cause GCR intensity depressions, seen both in the observations and simulations. We conclude that both drift and diffusion play important roles in accounting for GCR intensity depressions. The question is, which mechanism dominates for the two sample cases shown in Figures 10 and 11. GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

2423

Journal of Geophysical Research: Space Physics 4

10.1002/2013JA019546

(A)

Krr K//

GCR (%)

2

0

-2 -4 2007.78 4

2007.80

2007.82

2007.84

2007.86

2007.88

2007.90

(B)

GCR (%)

2

0

-2 -4 2009.24

2009.26

2009.28

2009.30

2009.32

2009.34

2009.36

Figure 12. GCR (E = 400 MeV) intensity variations (black), the radial diffusion component 𝜅rr (solid blue), and the parallel diffusion coefficient 𝜅∥ (dotted blue) during periods shown in (a) Figures 10 and (b) 11. The stream interfaces (SI1 and SI2 ) and sector boundaries (SB1 and SB2 ) are labeled.

To make the call, we plot the radial component of the diffusion tensor 𝜅rr for both periods in Figure 12. Figure 12a shows that for period 1 the GCR intensity and the diffusion coefficient are well correlated. The GCR intensity and 𝜅rr both reach their maxima near the stream interfaces and subsequently decrease abruptly. The highest value of 𝜅rr is about 3 times its smallest value. As we noted previously, GCR intensity variations associated with the passages of the HCS are much smaller than those associated with stream interfaces during this period. Therefore, drift-related depressions are less distinct than those caused by diffusion changes. We conclude that diffusive effects dominate the depressions of GCR intensity during period 1. For comparison, we show the parallel diffusion coefficient 𝜅∥ in same figure. One can see that diffusion drops and recovers at stream interface SI2 (recovery is not observed for SI1 due to insufficient temporal resolution), which is associated with the increased turbulence caused by compression and velocity shear [Intriligator et al., 2001]. The corresponding modest 1 day GCR depression was observed near SI2 , which is classified as a short-term depression. The following long-term depression has a duration of about 7 days. These two-stage depressions agree with the observations as described in the analysis of Figure 10 above. Figure 12b shows the variations of the GCR intensity, 𝜅rr , and 𝜅∥ for the interval shown in Figure 11. Comparing with Figure 12a, we find that the maxima of 𝜅rr and 𝜅∥ are much lower, and the ratio between the peak and the valley of 𝜅rr decreases to about 1.5. The decrease in diffusion is from larger Asl caused by denser slow solar wind dominating the low latitudes during period 2 (cf. equation (5)). As a result, stream interfaces have smaller role than in Figure 12a. For example, during the passages of SB2 and SI2 , the maxima of GCR intensity occur near the HCS. There is a large depression after the passage of the HCS. The stream interfaces play a secondary role, indicating that diffusion is not as important as drift during this period. In the end, GCR variations are caused by an interplay between diffusion and drift effects, which are represented by the passages of stream interfaces and the current sheet, respectively. In order to isolate the effects of CIRs on GCR modulation, we carried out three additional simulation, in which the CIRs were switched off completely. Here the flow pattern of fast and slow solar wind was computed with the HCS tilt angle 𝛾 = 0, but the wavy HCS was computed with the actual value of 𝛾 as in the full CIR simulations in order to retain the drift transport. The resulting spectra are shown in Figure 13. As mentioned previously, GCR intensity measured at 1 AU increases with a decreasing HCS tilt angle. The green curve gives the result obtained with a flat current sheet. We then increased the HCS tilt angle to the values reported for periods 1 and 2. The results are shown with dashed lines in the same figure. Comparing the GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

2424

Journal of Geophysical Research: Space Physics

10.1002/2013JA019546

Spectrum Intensity (/sec-m2-Mev-sterad)

results obtained with the same HCS tilt angle with CIRs (solid curves) and without CIRs (dashed curves), we see that CIRs actually facilitated GCR 1.0 transport into the inner heliosphere during both time intervals. Taking period 1 as an example, the intensity of 450 MeV protons was 34% of the unmodulated intensity in the simulation containing CIRs (red solid line) versus 25% in the simulation without CIRs (red dashed line). For period 0.1 2, the two intensity ratios with and without CIRs were 36% and 27%, respectively. There is an apparent increase of about 9% for the intensity of 450 MeV GCR after CIRs are taken into account. The difference is the 0.01 0.10 1.00 10.00 largest (∼12%) at 1 GeV. The effect is Energy (GeV) evidently due to a difference in diffuFigure 13. Energy spectra of GCR protons obtained near Earth using sion coefficients shown in Figure 12. models with CIRs (solid lines) and without CIRs (dashed lines) during the The net modulation effects of CIRs A < 0 solar minimum for two different tilt angles (plotted with differon GCR would be zero if the diffuent colors). As in Figure 7, the vertical dotted lines indicate the range of sion and drift coefficients remained PAMELA observations. unchanged. In that case energetic particles experience energy gains at compressions and losses at rarefaction regions, in turn, as they pass through CIRs. As the HCS tilt angle increases, the efficiency of GCR drift along the HCS decreases, but the diffusion coefficient increase compensates for that. As seen in Figure 12, despite local decreases in diffusion coefficients that act as “modulation barrier” inside narrow regions near CIR stream interfaces (e.g., SI2 in the top panel), the overall effects of CIRs is to facilitate diffusive GCR transport owing to an increase in diffusion coefficients in much of the modulation volume. Our simulation results indicates that the net modulation effects of CIRs are significant for GCR transport due to the associated enhancement in turbulence levels.

4. Discussion Previous statistical analyses of cosmic ray observations have shown that their intensity depressions are, in general, closely associated with the leading edges (forward waves or forward shocks) of CIRs and with their stream interfaces. The corotating structure of the HCS is another important factor that may cause GCR variations at a fixed location owing to a latitudinal gradient of GCR intensities [Kota and Jokipii, 1995]. However, the HCS appears to have a less dramatic effect, and the observed intensities do not always track the passages of the HCS [cf. Alania et al., 2011]. Earlier publications on the subject presented conflicting results (see the review paper by Richardson [2004]). In this work we reproduced the GCR recurrent variations using computer simulations that combine global MHD and stochastic transport models. Our simulation results that include variations of solar wind plasma properties and GCR intensities were compared with the observations during the recent cycle 23/24 solar minimum. Based on these results, it transpires that the tilt angle of the HCS is the deciding factor that determines whether stream interfaces or HCS crossings have a more prominent effect on GCR intensity depressions. When the tilt angle was relatively large (32.5◦ in period 1), stream interfaces were the dominant factor. Conversely, HCS crossings became more important for relatively small tilt angles (22◦ during period 2). The simulation results reveal that the passages of HCS (or stream interfaces) still cause minor GCR variations when stream interfaces (HCS) crossings dominate. This may explain the conflicting interpretations of the trigger mechanism for GCR depression onset based on observations. Despite a qualitative agreement between observations and our simulations (see Figures 10 and 11), one should keep in mind that the simulation results were obtained for an idealized scenario and will not agree with the observations in detail. The anticorrelation between GCR intensity and the solar wind speed GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

2425

Journal of Geophysical Research: Space Physics

10.1002/2013JA019546

reported from observations is not well reproduced in our simulations. The enhanced solar wind speed does not contribute to simulated GCR intensity depressions. For this model we did not make the dominant convection assumption [Richardson et al., 1996], where the scattering mean-free path is independent of azimuth and the solar wind speed is the only parameter varying with azimuth. The simplified treatment in their work resulted in good anticorrelation between the solar wind speed and cosmic ray intensity. It was empirically found that the product of the solar wind velocity 𝐮 and the IMF 𝐁 tends to be proportional to the modulation parameter 𝜁 ∼ uB [e.g., Alania et al., 2011]. From the transport theory, the modulation parameter is 𝜁 ∼ uL∕𝜅 , where L is the characteristic length and 𝜅 the diffusion coefficient. If 𝜅 ∼ 1∕B, which is often assumed in transport models, then indeed 𝜁 ∼ uB. Our model does not appear to exhibit this so-called VB effect. In our simulations 𝜅 is based on quasi-linear theory and its dependence on the magnetic field is approximately B5∕3 ∕⟨𝛿B2 ⟩ for small heliocentric distances and excluding the polar regions. We plan to experiment with different expressions for the diffusion coefficient in the future but prefer to keep the present model physics based. One idealization is the predetermined CIR geometry imposed at the inner boundary, where the fast and slow solar wind streams are separated by a wide transition in the corotating frame (see Figure 1). The Sun was mostly quiet during the solar minimum, but recurrent transient events were observed, making the evolution and propagation of the solar wind more complicated than in our simple model. The simulated wavy HCS structure is oversimplified, as can be deduced from Figures 8 and 9. There were more observed HCS crossings than in the simulations because multiple crossings may occur during one solar rotation. Changes in the tilt angle during each period were also ignored. As a consequence, the inward drift rate of GCR ions could be smaller than in the model. Note also that CIR shock pairs are not fully developed by 1 AU in our model, as shown in Figures 4 and 5, which disagrees with some observations. Consequently, in reality the role of the HCS tilt angle may be less essential than in our simulations. A second serious simplification is the assumption of constant HCS tilt angles during each half of the solar minimum. From observations, there was almost no tilt angle change during period 1 (Figure 8 (top)) but a pronounced decrease and recovery during period 2 (Figure 9 (top)). These gradual changes were probably responsible for the simultaneous GCR intensity increase. However, the subject is outside the scope of this paper, which focuses on 27 day variations. Because the present model assumes a corotating structure of the flow, plasma results from a single time snapshot is sufficient to simulate GCR transport over any length of time by simply rotating the background. Relaxing the assumption of constant tilt angle would require frequent plasma simulation data sampling from many solar rotations with varying tilt angles, dramatically increasing the storage demand. For example, 1 day sampling would incur over 1000 snapshots for the complete solar minimum, requiring of the order of 1 TB of disk storage. Subsequent work will determine whether less frequent sampling might suffice for a longer-term modulation problem. It was reported that the diffusion coefficient apparently increased during period 2, and the corresponding average value is expected to be larger than that in period 1 from observations [Mewaldt et al., 2010]. This scenario is somewhat different from our results shown in Figure 12 and may be caused by two factors. First, the simulated HCS tilt angle was fixed during the two periods for simplicity, leading to the absence of diffusion change in response to the observed gradual change of the tilt angle (blue curve in Figure 9) during period 2. When the tilt angle dropped below 30◦ , the fraction of fast solar wind reaching the equatorial plane at 1 AU decreases, as does the amount of velocity shear in the flow. Because shear is an important source of turbulence [Zank et al., 1998], its absence would lead to an increase in the parallel diffusion coefficient near the ecliptic. However, shear persists at high latitudes [Erdos and Balogh, 2005], and so the parallel diffusion coefficient there is expected to be lower. The change in diffusion is therefore nonuniform and required a better understanding of the time-dependent turbulence transport process. Second, the averaged plasma properties of slow solar wind are based on statistical results covering the period from July 2008 to June 2009 [Jian et al., 2011] and may not fully representative of the first half of the solar minimum. The fast solar wind properties based on Ulysses-averaged data [Ebert et al., 2009] may not be sufficiently accurate for period 2, because Ulysses was at low latitudes (less than 30◦ ) since the beginning of 2009. The actual solar wind is certainly much more complicated than in our simplified model. However, we believe our simple, but physically consistent, approach does lead to more insight and understanding of the physics of the 27 day variations in cosmic ray intensities. GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

2426

Journal of Geophysical Research: Space Physics

10.1002/2013JA019546

Based on previous investigations [Richardson, 2004], CIR-related “magnetic barriers” or magnetic field enhancements are not always a critical factor organizing GCR intensity depressions. This is consistent with our simulation results shown in Figures 10 and 11. Whereas the magnetic field strength peaks at stream interfaces, the enhancements in the magnetic field prior to stream interfaces do not trigger the onset of GCR depressions. The role of the magnetic field is less prominent than in the works of Burlaga et al. [1985] and Kota and Jokipii [1991], where GCR intensity decrease were closely associated with magnetic field enhancements. However, contrary to our simulations, turbulent intensity variations were not taken into account in these models. Observations show that magnetic field turbulence is far from uniformly distributed within CIRs, and not simply related to the magnetic field strength [Belcher and Davis, 1971]. In the present work, the mean square turbulent magnetic field intensity ⟨𝛿B2 ⟩ was proportional to the plasma density. The resulting turbulence ratio ⟨𝛿B2 ⟩∕B2 is not exclusively dependent on the magnetic field, which is different from the diffusion coefficients used by Kota and Jokipii [1991]. The two-stage GCR depressions were identified in both data and simulations. The short-term depression has a duration typically less than 1 day and is associated with the passage of a stream interface. Generally, it is caused by a decrease in the diffusion coefficient within the CIR and in the vicinity of a stream interface as a results of compression and shear in solar wind plasma [Intriligator et al., 2001]. The following long-term depression lasts for several days (about 7 days in our simulations). As shown in Figure 12, this depression is associated with the long-term decrease of the diffusion coefficient. Fast and slow solar winds have different diffusion coefficients near Earth, and their alternation is the primary cause of the long-term diffusion variation. In this work we focused on GCR and diffusion variations at 1 AU. Investigating GCR modulation by CIRs and CMIRs at larger heliocentric distances will be the subject of our future modeling efforts.

5. Conclusions 1. We carried out a modeling investigation to determine the dominant mechanism(s) responsible for the 27 day periodic variations of galactic cosmic ray ions observed near Earth during the recent negative polarity solar minimum. The investigation involved global MHD simulations of corotating interaction regions and stochastic trajectory simulation of GCR particles. 2. Simulation results for the variations of solar wind plasma properties and ion intensities were compared with the observations performed over the 3 years, 2007 through 2009. The simulated variations of 400 MeV galactic protons are consistent with the results from ground level neutron monitors, having a depression amplitude (peak to valley) of 1–3% during the past solar minimum. 3. Our results show that stream interfaces and the associated depressions in the diffusion mean-free path are more important than the drift effects of the current sheet during periods characterized by relatively large tilt angles of the solar magnetic axis. The effect of the HCS becomes more prominent with decreasing tilt angle. This complex interplay between diffusion and drift appears to be responsible for the previously conflicting results about the roles of stream interfaces and the HCS in the onset of GCR intensity depressions. 4. The effect of CIRs on the amount of modulation at 1 AU was investigated by performing simulations without the corotating plasma structure but with the same HCS tilt angles. It was found that CIRs did facilitate the transport of GCR significantly during the recent cycle 23/24 solar minimum (without CIRs, GCR intensities would be about 10% lower). Acknowledgments This work was supported by NASA under grants NNX13AF99G, NNX10AE46G, and NNX12AH44G, and by NSF under grant AGS-0955700. We acknowledge the use of the Neutron Monitor Database, the Wilcox Solar Observatory, and the OMNI data obtained through the GSFC/SPDF OMNIWeb interface at http://omniweb. gsfc.nasa.gov. Philippa Browning thanks Horst Fichtner and an anonymous reviewer for their assistance in evaluating this paper.

GUO AND FLORINSKI

References Ackermann, M., et al. (2013), Detection of the characteristic pion-decay signature in supernova remnants, Science, 339, 807–811, doi:10.1126/science.1231160. Adriani, O., et al. (2013), Time dependence of the proton flux measured by PAMELA during the 2006 July – 2009 December solar minimum, Astrophys. J., 765, 91, doi:10.1088/0004-637X/765/2/91. Alania, M. V., R. Modzelewska, and A. Wawrzynczak (2011), On the relationship of the 27-day variations of the solar wind velocity and galactic cosmic ray intensity in minimum epoch of solar activity, Sol. Phys., 270, 629–641. Alanko, K., I. G. Usoskin, K. Mursula, and G. A. Kovaltsov (2003), Heliospheric modulation strength: Effective neutron monitor energy, Adv. Space Res., 32, 615–620. Altschuler, M. D., and G. Newkirk (1969), Magnetic fields and structure of solar corona. I. Methods of calculating coronal fields, Sol. Phys., 9, 131–149. Ball, B., M. Zhang, H. Rassoul, and T. Linde (2005), Galactic cosmic-ray modulation using a solar minimum MHD heliosphere: A stochastic particle approach, Astrophys. J., 634, 1116, doi:10.1086/496965. Bieber, J. W., W. Wanner, and W. H. Matthaeus (1996), Dominant two-dimensional solar wind turbulence with implications for cosmic ray transport, J. Geophys. Res., 101, 2511–2522.

©2014. American Geophysical Union. All Rights Reserved.

2427

Journal of Geophysical Research: Space Physics

10.1002/2013JA019546

Belcher, J. W., and L. Davis Jr. (1971), Large-amplitude Alfvén waves in the interplanetary medium, 2, J. Geophys. Res., 76(16), 3534–3563, doi:10.1029/JA076i016p03534. Borovikov, S. N., N. V. Pogorelov, and R. W. Ebert (2012), Solar rotation effects on the heliosheath flow near solar minima, Astrophys. J., 750, 42, doi:10.1088/0004-637X/750/1/42. Breech, B., W. H. Matthaeus, J. Minnie, J. W. Bieber, S. Oughton, C. W. Smith, and P. A. Isenberg (2008), Turbulence transport throughout the heliosphere, J. Geophys. Res., 113, A08105, doi:10.1029/2007JA012711. Burlaga, L. F., R. Schwenn, and H. Rosenbauer (1983), Dynamical evolution of interplanetary magnetic fields and flows between 0.3 AU and 8.5 AU: Entrainment, Geophys. Res. Lett., 10, 413–416. Burlaga, L. F., F. B. McDonald, N. F. Ness, R. Schwenn, A. J. Lazarus, and F. Mariani (1984), Interplanetary flow systems associated with cosmic-ray modulation in 1977–1980, J. Geophys. Res., 89, 6579–6587. Burlaga, L. F., F. B. McDonald, M. L. Goldstein, and A. J. Lazarus (1985), Cosmic-ray modulation and turbulent interaction near 11 AU, J. Geophys. Res., 90, 12,027–12,039. Burlaga, L. F., F. B. McDonald, and N. F. Ness (1993), Cosmic ray modulation and the distant heliospheric magnetic field: Voyager 1 and 2 observations from 1986 to 1989, J. Geophys. Res., 98, 1–11. Burlaga, L. F., N. F. Ness, F. B. McDonald, J. D. Richardson, and C. Wang (2003), Voyager 1 and 2 observations of magnetic fields and associated cosmic-ray variations from 2000 through 2001: 60 - 87 AU, Astrophys. J., 582, 540, doi:10.1086/344571. Burlaga, L. F., N. F. Ness, Y.-M. Wang, N. R. Sheeley Jr., and J. D. Richardson (2010), Observations of the magnetic field and plasma in the heliosheath by Voyager 2 from 2007.7 to 2009.4, J. Geophys. Res., 115, A08107, doi:10.1029/2009JA015239. Czechowski, A., M. Strumik, J. Grygorczuk, S. Grzedzielski, R. Ratkiewicz, and K. Scherer (2010), Structure of the heliospheric current sheet from plasma convection in time-dependent heliospheric models, Astron. Astrophys., 516, A17, doi:10.1051/0004-6361/200913542. Ebert, R. W., D. J. McComas, H. A. Elliott, R. J. Forsyth, and J. T. Gosling (2009), Bulk properties of the slow and fast solar wind and interplanetary coronal mass ejections measured by Ulysses: Three polar orbits of observations, J. Geophys. Res., 114, A01109, doi:10.1029/2008JA013631. Erdos, G., and A. Balogh (2005), In situ observations of magnetic field fluctuations, Adv. Space Res., 35, 625–635. Dröge, W. (2000), Particle scattering by magnetic fields, Space Sci. Rev., 93, 121–151, doi:10.1023/A:1026588210726. Florinski, V., and N. V. Pogorelov (2009), Four-dimensional transport of galactic cosmic rays in the outer heliosphere and heliosheath, Astrophys. J., 701, 642–651, doi:10.1088/0004-637X/701/1/642. Florinski, V., F. Alouani-Bibi, J. Kota, and X. Guo (2012), Cosmic-ray diffusion in a sectored magnetic field in the distant heliosheath, Astrophys. J., 754, 31, doi:10.1088/0004-637X/754/1/31. Florinski, V., X. Guo, D. S. Balsara, and C. Meyer (2013), Magnetohydrodynamic modeling of solar system processes on geodesic grids, Astrophys. J. Suppl., 205, 19, doi:10.1088/0067-0049/205/2/19. Forbush, S. E. (1938), On cosmic-ray effects associated with magnetic storms, Terr. Magn. Atmos. Electr., 43, 203–218. Gazis, P. R. (2000), A large-scale survey of corotating interaction regions and their successors in the outer heliosphere, J. Geophys. Res., 105, 19–33. Giacalone, J., and J. R. Jokipii (1999), The transport of cosmic rays across a turbulent magnetic field, Astrophys. J., 520, 204–214. Gosling, J. T, A. J. Hundhausen, V. Pizzo, and J. R. Asbridge (1972), Compressions and rarefactions in the solar wind: Vela 3, J. Geophys. Res., 77, 5442–5454. Gosling, J. T, J. R. Asbridge, S. J. Bame, and W. C. Feldman (1978), Solar wind stream interfaces, J. Geophys. Res., 83, 1401–1412. Heber, B., T. R. Sanderson, and M. Zhang (1999), Corotating interaction regions, Adv. Space Res., 23, 567–579. Hundhausen, A. J., and J. T. Gosling (1976), Solar wind structure at large heliocentric distances: An interpretation of Pioneer 10 observations, J. Geophys. Res., 81, 1436–1440. Intriligator, D. S., J. R. Jokipii, T. S. Horbury, J. M. Intriligator, R. J. Forsyth, H. Kunow, G. Wibberenz, and J. T. Gosling (2001), Processes associated with particle transport in corotating interaction regions and near stream interfaces, J. Geophys. Res., 106, 10,625–10,634, doi:10.1029/2000JA000070. Iucci, N., M. Parisi, M. Storini, and G. Villoresi (1979), High speed solar wind streams and galactic cosmic ray modulation, Il Nuovo Cimento C, 2, 421–438. Jian, L. K., C. T. Russell, and J. G. Luhmann (2011), Comparing solar minimum 23/24 with historical solar wind records at 1 AU, Sol. Phys., 274, 321–344, doi:10.1007/s11207-011-9737-2. Jokipii, J. R., and E. H. Levy (1977), Effects of particle drifts on the solar modulation of galactic cosmic rays, Astrophys. J., 213, L85—L88. Jokipii, J. R., and D. A. Kopriva (1979), Effects of particle drift on the transport of cosmic rays. III. Numerical models of galactic cosmic-ray modulation, Astrophys. J., 234, 384–392. Jokipii, J. R., and B. Thomas (1981), Effects of drift on the transport of cosmic rays. IV. Modulation by a wavy interplanetary current sheet, Astrophys. J., 243, 1115–1122. Kota, J., and J. R. Jokipii (1983), Effects of particle drift on the transport of cosmic rays. VI. A three-dimensional model including diffusion, Astrophys. J., 265, 573–581. Kota, J., and J. R. Jokipii (1991), The role of corotating interaction regions in cosmic-ray modulation, Geophys. Res. Lett., 18, 1979–1800. Kota, J., and J. R. Jokipii (1995), Corotating variations of cosmic rays near the south heliospheric pole, Science, 268, 1024–1025, doi:10.1126/science.268.5213.1024. Kota, J., and J. R. Jokipii (2000), Velocity correlation and the spatial diffusion coefficients of cosmic rays: Compound diffusion, Astrophys. J., 531, 1067–1070, doi:10.1086/308492. Langner, U. W., and M. S. Potgieter (2004), Effects of the solar wind termination shock and heliosheath on the heliospheric modulation of galactic and anomalous Helium, Ann. Geophys., 22, 3063–3072, doi:10.5194/angeo-22-3063-2004. Lazarus, A. J., J. D. Richardson, R. B. Decker, and F. B. McDonald (1999), Voyager 2 observations of corotating interaction regions (CIRs) in the outer heliosphere, Space Sci. Rev., 89, 53–59. le Roux, J. A., and H. Fichtner (1999), Global merged interaction regions, the heliospheric termination shock, and time-dependent cosmic ray modulation, J. Geophys. Res., 104, 4709–4730, doi:10.1029/1998JA900089. le Roux, J. A., G. P. Zank, and V. S. Ptuskin (1999), An evaluation of perpendicular diffusion models regarding cosmic ray modulation on the basis of a hydromagnetic description for solar wind turbulence, J. Geophys. Res., 104, 24,845–24,862. le Roux, J. A., and G. M. Webb (2007), Nonlinear cosmic-ray diffusive transport in combined two-dimensional and slab magnetohydrodynamic turbulence: A BGK-Boltzmann approach, Astrophys. J., 667, 930–955. Leske, R. A., A. C. Cummings, R. A. Mewaldt, E. C. Stone, T. T. von Rosenvinge, and M. E. Wiedenbeck (2011), Observed 27-day variations in cosmic ray intensities during the cycle 23/24 solar minimum, Proc. 32 Int. Cosmic Ray Conf., 11, 194–197.

GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

2428

Journal of Geophysical Research: Space Physics

10.1002/2013JA019546

Leske, R. A., A. C. Cummings, R. A. Mewaldt, and E. C. Stone (2013), Anomalous and galactic cosmic rays at 1 AU during the cycle 23/24 solar minimum, Space Sci. Rev., 176, 253–263, doi:10.1007/s11214-011-9772-1. Luo, X., M. Zhang, H. K. Rassoul, and N. V. Pogorelov (2011), Cosmic-ray modulation by the global merged interaction region in the heliosheath, Astrophys. J., 730, 13, doi:10.1088/0004-637X/730/1/13. McKibben, R. B., et al. (1999), Modulation of cosmic rays and anomalous components by CIRs, Space Sci. Rev., 89, 307–326. Mewaldt, R. A., et al. (2010), Record-setting cosmic-ray intensities in 2009 and 2010, Astrophys. J. Lett., 723, L1, doi:10.1088/2041-8205/723/1/L1. Modzelewska, R, and M. V. Alania (2013), The 27-day cosmic ray intensity variations during solar minimum 23/24, Sol. Phys., 286, 593–607, doi:10.1007/211207-013-0261-4. Parker, E. N. (1958), Dynamics of the interplanetary gas and magnetic fields, Astrophys. J., 128, 664–676. Park, E. N. (1965), The passage of energetic charged particles through interplanetary space, Planet. Space Sci., 13, 9–49. Perko, J. S., and L. A. Fisk (1983), Solar modulation of galactic cosmic rays: 5. Time-dependent modulation, J. Geophys. Res., 88, 9033–9036. Pizzo, V. J. (1991), The evolution of corotating stream fronts near the ecliptic plane in the inner solar system, 2. Three-dimensional tilted-dipole fronts, J. Geophys. Res., 96, 5405–5420. Pizzo, V. J. (1994), Global, quasi-steady dynamics of the distant solar wind, 1. Origin of north-south flows in the outer heliosphere, J. Geophys. Res., 99, 4173–4183. Potgieter, M. S., and H. Moraal (1985), A drift model for the modulation of galactic cosmic rays, Astrophys. J., 294, 425–440. Potgieter, M. S., and J. A. le Roux (1994), The long-term heliospheric modulation of galactic cosmic rays according to a time-dependent drift model with merged interaction regions, Astrophys. J., 423, 817–827. Provornikova, E., M. Opher, V. Izmodenov, and G. Toth (2012), Do corotating interaction region associated shocks survive when they propagate into the heliosheath, Astrophys. J. Lett., 756, L37, doi:10.1088/2041-8205/756/2/L37. Qin, G, W. H. Matthaeus, and J. W. Bieber (2002), Perpendicular transport of charged particles in composite model turbulence: Recovery of diffusion, Astrophys. J., 578, L117—L120, doi:10.1086/344687. Richardson, I. G., G. Wibberenz, and H. V. Cane (1996), The relationship between recurring cosmic ray depressions and corotating solar wind streams at ≤ 1 AU: IMP 8 and Helios 1 and 2 anticoincidence guard rate observations, J. Geophys. Res., 101, 13,483–13,496. Richardson, I. G. (2004), Energertic particles and corotating interaction regions in the solar wind, Space Sci. Rev., 111, 267–376. Richardson, J. D., and C. Wang (2011), Plasma in the heliosheath: 3.5 years of observations, Astrophys. J., 734, L21, doi:10.1088/2041-8205/734/1/L21. Riley, P., J. A. Linker, and Z. Mikic (2001), An empirically-driven global MHD model of the solar corona and inner heliosphere, J. Geophys. Res., 106, 15,889–15,901. Rouillard, A. P., and M. Lockwood (2007), The latitudinal effect of corotating interaction regions on galactic cosmic rays, Sol. Phys., 245, 191–206, doi:10.1007/s11207-007-9019-1. Sarabhai, V. (1963), Some consequences of nonuniformity of solar wind velocity, J. Geophys. Res., 68, 1555–1557, doi:10.1029/JZ068i005p01555. Simpson, J. A. (1998), A brief history of recurrent solar modulation of the galactic cosmic rays (1937-1990), Space Sci. Rev., 83, 169–176. Stone, E. C., A. C. Cummings, F. B. McDonald, B. C. Heikkila, N. Lai, and W. R. Webber (2005), Voyager 1 explores the termination shock region and the heliosheath beyond, Science, 309, 2017–2020, doi:10.1126/science.1117684. Stone, E. C., A. C. Cummings, F. B. McDonald, B. C. Heikkila, N. Lai, and W. R. Webber (2009), An asymmetric solar wind termination shock, Nature, 454, 71–74, doi:10.1038/nature07022. Strauss, R. D., M. S. Potgieter, I. Busching, and A. Kopp (2012), Modelling heliospheric current sheet drift in stochastic cosmic ray transport models, Astrophys. Space Sci., 339, 223–236, doi:10.1007/s10509-012-1003-z. Strauss, R. D., M. S. Potgieter, S. E. S. Ferreira, H. Fichtner, and K. Scherer (2013), Cosmic ray modulation beyond the heliopause: A hygrid modeling approach, Astrophys. J. Lett., 765, L18, doi:10.1088/2041-8205/765/1/L18. Siscoe, G. L. (1972), Structure and orientation of solar wind interaction fronts: Pioneer 6, J. Geophys. Res., 77, 27–34. Usmanov, A. V., and M. L. Goldstein (2006), A three-dimensional MHD solar wind model with pickup protons, J. Geophys. Res., 111, A07101, doi:10.1029/2005JA011533. Webber, W. R., A. C. Cummings, F. B. McDonald, E. C. Stone, B. Heikkila, and N. Lal (2008), Galactic cosmic ray H and He nuclei energy spectra measured by Voyagers 1 and 2 near the heliospheric termination shock in positive and negative solar magnetic polarity cycles, J. Geophys. Res., 113, A10108, doi:10.1029/2008JA013395. Webber, W. R., and P. R. Higbie (2010), What Voyager cosmic ray data in the outer heliosphere tells us about 10Be production in the Earths polar atmosphere in the recent past, J. Geophys. Res., 115, A05102, doi:10.1029/2009JA014532. Zank, G. P., W. H. Matthaeus, and C. W. Smith (1996), Evolution of turbulent magnetic fluctuation power with heliospheric distance, J. Geophys. Res., 101, 17,093–17,107, doi:10.1029/96JA01275. Zank, G. P., W. H. Matthaeus, J. W. Bieber, and H. Moraal (1998), The radial and latitudinal dependence of the cosmic ray diffusion in the heliosphere, J. Geophys. Res., 103, 2085–2097. Zhang, M. (1999), A Markov stochastic process theory of cosmic-ray modulation, Astrophys. J., 513, 409–420, doi:10.1086/306857.

GUO AND FLORINSKI

©2014. American Geophysical Union. All Rights Reserved.

2429