GALACTIC COSMIC-RAY MODULATION IN A REALISTIC GLOBAL MAGNETOHYDRODYNAMIC HELIOSPHERE

Equation (2) is the governing MHD equation (Pogorelov et al. 2006) for the plasma flow in the heliosphere. It contains the conservation laws for mass, momentum, total energy, and magnetic flux. As for the interstellar neutral atoms, a kinetic approach is employed. The Boltzmann equation (Equation (3)) is adopted to describe these neutral atoms (Heerikhuisen et al. 2006, 2008). These two sets of equations are coupled through the source terms H^m_{p − H} and H^e_{p − H}, which illustrate the momentum and energy transfer between the neutral atoms and charged particles:

$$\begin{eqnarray} &&\frac{\partial }{\partial t}\left[ \begin{array}{@{}c@{}}\rho \\ \rho \boldsymbol {v}\\ e\\ \displaystyle\boldsymbol {B} \end{array} \right] + \nabla \cdot \left[ \begin{array}{@{}c@{}}\rho \boldsymbol {v}\\ \displaystyle\rho \boldsymbol {v}\boldsymbol {v}+p_{0}\boldsymbol {I}-\frac{\boldsymbol {B}\boldsymbol {B}}{4\pi } \\[6pt] \displaystyle(e+p_{0})\boldsymbol {v} - \boldsymbol {B}\frac{(\boldsymbol {v}\cdot \boldsymbol {B})}{4\pi } \\ \displaystyle\boldsymbol {v}\boldsymbol {B}-\boldsymbol {B}\boldsymbol {v} \end{array} \right]\nonumber\\ &&\quad = -\left[\begin{array}{@{} c @{}} \displaystyle0\\[5pt] \displaystyle\frac{\boldsymbol {B}}{4\pi }\\[8pt] \displaystyle\frac{\boldsymbol {v}\cdot \boldsymbol {B}}{4\pi }\\[6pt] \displaystyle\boldsymbol {v} \end{array} \right] \times \nabla \cdot \boldsymbol {B} + \left[\displaystyle\begin{array}{@{}c@{}} \displaystyle 0\\ \displaystyle H^{m}_{p-H}\\[6pt] \displaystyle H^{e}_{p-H}\\ \displaystyle 0 \end{array} \right] \end{eqnarray} \tag{ 2 }$$

$$\begin{equation} \frac{\partial }{\partial t} f_H + \boldsymbol {v}\cdot \nabla f_{H} + \frac{F}{m_p}\nabla _{v}\cdot f_{H} = P-L. \end{equation} \tag{ 3 }$$

The heliospheric MHD data are obtained by solving the coupled MHD-kinetic equations using the finite volume method (Pogorelov & Matsuda 1998). Our coordinate system is defined as the X-axis being along the solar rotation axis; the Z-axis belongs to the plane formed by the X-axis and the plasma velocity vector in the unperturbed local interstellar medium (LISM) flow, and points upstream into it. The Y-axis completes this orthogonal coordinate system. After the coordinate system is constructed, we choose the mesh grid in the associated spherical coordinates according to

$$\begin{eqnarray} &&R_i = 10 + \frac{1190}{2\times (e^{4}-1)}\left[e^{(\frac{i-1}{5})} + e^{(\frac{i}{55})} - 2\right] \;\;\; i=1,2,\ldots ,220\nonumber\\ \end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray} &&\theta _j = (j-0.5)\times \frac{\pi }{140} \;\;\; j=1,2, 3,\ldots , 140\quad\qquad\qquad \,\, \, \end{eqnarray} \tag{ 4b }$$

$$\begin{eqnarray} &&\phi _k = (k-1) \times \frac{2\pi }{79} \;\;\; k=1, 2, 3, \ldots , 79.\qquad\qquad\qquad \, \, \, \end{eqnarray} \tag{ 4c }$$

The polar angle θ is measured from the Z-axis, while the azimuthal angle ϕ starts from the X-axis. These angles will be converted into latitude and longitude in a heliographic system. 0° latitude is roughly in the solar equator, while 0° longitude is the direction toward the nose of the heliosphere.

These equations define a unique mesh system, and the whole spherical simulation domain with a radius of 1200 AU was divided into 220 × 140 × 79 cells. The numerical simulation allocates a unique value of solar wind density, temperature, plasma velocity, and magnetic field, namely our MHD data, into each cell.

Figure 1 shows the magnetic field data in the X–Z plane used for this study. In comparison, MHD data from Linde et al. (1998) are shown in Figure 2 (the Z-axis in Linde's MHD data is along the solar rotation axis and the X-axis is along the LISM's movement velocity and we have already rotated the figure for easy comparison). Similarly, Figures 3 and 4 show the solar wind speed of both MHD data sets at the X–Z plane. It is clear that after considering the momentum and energy transfer caused by the interstellar neutral atoms, the profile of the TS is significantly changed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1 illustrates the location of the TS in these MHD models along four different locations. It can be seen that the TS in the ecliptic plane is largely expanded in our MHD model compared with the model of Linde et al. (1998). The upwind location was moved from 70 AU to 91 AU, being more consistent with the crossing location of Voyager 1 (Stone et al. 2005) and the recent estimation made by Luo et al. (2011). Compared with the expansion in the ecliptic plane, the size of the TS in the polar direction has been reduced, which is due to the symmetrized effect caused by interstellar neutrals (Pogorelov et al. 2007). In the subsequent section, we will show that these changes affect the Galactic cosmic-ray modulation as well.

Both Figures 3 and 4 demonstrate that the TS is a closed boundary where the solar wind speed decreases from 400 km s⁻¹ to nearly 130 km s⁻¹. Inside the TS is the supersonic solar wind region where the solar wind speed is nearly symmetric.

As for the magnetic field profile, a difference between these two sets of MHD data also exists. The dark blue color denotes a vanishing magnetic field. Near the nose of the TS, a blue belt turns upward toward the northern area in our MHD data. This is the ISM-induced bending of the heliospheric current sheet, which is assumed to be flat in the supersonic SW ahead of the TS. The magnetic field will be mainly used to calculate cosmic-ray drift speed. The drift speed is very large near the current sheet because of the sudden reverse of the magnetic field direction. This is consistent with the well-known concept that cosmic rays experience very large drift speed along current sheets (Burger & Potgieter 1989).

Figure 5 provides a more straightforward demonstration of how the magnetic field varies with radial distance in the heliosphere. It shows the linear distribution of the magnetic field magnitude along three directions: the upwind direction, the Voyager 1 trajectory at 35°N heliolatitude, and the Voyager 2 trajectory at 35°S heliolatitude. It can be seen from this plot that the magnitude along these three directions decreases consistently inside the supersonic solar wind region. After the TS crossing, the magnitude of the magnetic field increases due to shock compression. However, for the upwind direction, after a short increase of the magnitude, the magnetic field experiences a sudden decrease. This is caused by the crossing of the current sheet. This crossing of the current sheet can also be seen in the Voyager 1 direction. Similar to the current sheet, another structure called the magnetic wall can also be seen in these curves. At around 150–200 AU, all three directions experience a peak in magnetic field, which we call the "magnetic wall." According to our simulation shown later, this strong magnetic field acts as a barrier that can reduce cosmic rays coming in from the ISM.

Zoom In Zoom Out Reset image size

While three-dimensional Galactic cosmic-ray transport simulation code has been developed by several groups (Kota & Jokipii 1983; Potgieter & Burger 1990), our code possesses a unique technique that has many advantages when dealing with an MHD heliosphere model. Instead of solving the Parker equation directly with the finite difference method, we solve for the stochastic phase space trajectory of a particle as a function of backward time s. This process begins from the location where we need to know the cosmic-ray flux (Zhang 1999) until it hits the boundary of the entire simulation domain, which is the outer boundary of the LISM. The solution to the Parker transport equation is the expectation value of all boundary values when the particle exits from the boundaries. The stochastic differential equation for the random particle trajectory corresponding to the Parker transport equation can be written as (Zhang 1999)

$$\begin{eqnarray} &&d\boldsymbol {x} = \sqrt{\vec{\vec{2\kappa }}} \cdot d\boldsymbol {W}(s) + (\nabla \cdot \vec{\vec{\kappa }} - \boldsymbol {V}_{{\rm sw}} - \boldsymbol {V}_d)ds \end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray} &&dp = \frac{1}{3} \nabla \cdot \boldsymbol {V}_{{\rm sw}} p ds.\qquad\qquad\qquad\qquad \end{eqnarray} \tag{ 5b }$$

During each iterating step, the location and momentum of the pseudo-particle described by the stochastic differential equation will be updated. Here, dW is a three-dimensional random Wiener differential noise term, which can be generated using a Gaussian distribution random number with a mean of zero and a standard deviation of ds. $\vec{\vec{\kappa }}$ is the diffusion tensor, and it has the following form if expressed in the local magnetic field coordinates:

$$\begin{equation} \vec{\vec{\kappa }} = \left(\begin{array}{@{}ccc@{}}\kappa _{\bot } & 0 & 0 \\ 0 & \kappa _{\bot } & 0 \\ 0 & 0 & \kappa _{\parallel } \end{array} \right) \; \;. \end{equation} \tag{ 6 }$$

In our simulation, the diffusion tensor can be transferred to our spherical coordinates using the following form:

$${\begin{equation} \vec{\vec{\kappa }} = \left(\begin{array}{@{}ccc@{}}\kappa _{\bot }\cos ^{2}\psi +\sin ^{2}\psi \kappa _{\parallel }\! &\! 0 \!&\! -\cos \psi \sin \psi \kappa _{\bot } + \sin \psi \cos \psi \kappa _{\parallel }\\ 0\! &\! 1 \!&\! 0 \\ -\sin \psi \cos \psi \kappa _{\bot } + \cos \psi \sin \psi \kappa _{\parallel } \!&\! 0 \!&\! \sin ^{2}\psi \kappa _{\bot } + \cos ^{2}\psi \kappa _{\parallel } \end{array} \right), \end{equation}} \tag{ 7 }$$

where ψ is the angle between the local magnetic field direction B and the radial direction r. Transformations to other coordinate systems can be done in a similar manner.

For the diffusion coefficient, we choose the following general form:

$$\begin{eqnarray} &&\kappa _{\parallel } = \kappa _{\parallel 0}\beta \left(\frac{p}{p_0}\right)^{b_{\parallel }}\left(\frac{B_{{\rm eq}}}{B}\right)^{a_{\parallel }}, \end{eqnarray} \tag{ 8a }$$

$$\begin{eqnarray} &&\kappa _{\bot } = \kappa _{\bot 0}\beta \left(\frac{p}{p_0}\right)^{b_{\bot }}\left(\frac{B_{{\rm eq}}}{B}\right)^{a_{\bot }}, \end{eqnarray} \tag{ 8b }$$

where p and B are the particle momentum and magnetic field strength, respectively, and β is the ratio of particle speed to the speed of light. The p₀ parameter is a reference momentum (in our case 1 GeV c⁻¹), B_eq is the IMF strength at the heliospheric equator at 1 AU. The constants κ_∥0 and κ_⊥0 determine the magnitudes of parallel and perpendicular diffusion coefficients. They are chosen to be 50 × 10²⁰ cm² s⁻¹ and 5 × 10²⁰ cm² s⁻¹, respectively. The exponents are chosen as a_{∥, ⊥} = 0.5 and b_{∥, ⊥} = 1. The choice of the above diffusion coefficients is a little arbitrary, but it is approximately consistent with the overall modulation level inferred by various observations and a few other simulations of cosmic-ray modulation.

As for the drift velocity, we use the classical form (Jokipii et al. 1981)

$$\begin{equation*} \boldsymbol {V}_d = \frac{pv}{3q}\nabla \times \frac{\boldsymbol {B}}{B^2} \;. \end{equation*}$$

Both the diffusion coefficient and drift velocity are functions of magnetic field and plasma speed. Thus, following Equation (5), it is essential to know the magnetic field and solar wind speed in order to calculate the location and momentum increments. This is done through direct numerical differentiation of the magnetic field vector on the MHD grids. When cosmic rays transport to any random location inside the heliosphere, a function for calculating the local diffusion coefficient and drift velocity is called into the stochastic differential equation (5).

In order to obtain the magnetic field and solar wind speed at any arbitrary point (r, θ, ϕ) inside the simulation domain, the discrete MHD data need to be interpolated. The original MHD data are based on a three-dimensional spherical coordinate system. According to the mesh grid system (Equation (4)), the simulation globe is divided into 220 spherical planes and on each plane, the angle difference between the two neighboring points is 2π/79 along azimuthal direction and π/140 along the latitudinal direction.

Thus, for any arbitrary point (r, θ, ϕ), there are eight adjoining grid points in the MHD data set: (r_i, θ_i, ϕ_i), (r_i, θ_{i + 1}, ϕ_{i + 1}), (r_i, θ_{i + 1}, ϕ_i), (r_i, θ_i, θ_{i + 1}) on the spherical plane with a radius of r_i and (r_{i + 1}, θ_i, ϕ_i), (r_{i + 1}, θ_{i + 1}, ϕ_i), (r_{i + 1}, θ_i, ϕ_{i + 1}), (r_{i + 1}, θ_{i + 1}, ϕ_{i + 1}) on the spherical plane with a radius of r_{i + 1}. In addition, the following relationship exists for these quantities, r_i ⩽ r ⩽ r_{i + 1}, θ_i ⩽ θ ⩽ θ_{i + 1}, ϕ_i ⩽ ϕ ⩽ ϕ_{i + 1}.

Using the value of these eight points, our interpolation method is achieved in two steps: (1) use the bilinear interpolation method (Press et al. 1992) to obtain the values at points (r_{i − 1}, θ, ϕ) and (r_i, θ, ϕ) separately and (2) after obtaining the values at these two projecting points, employ a linear interpolation method (Press et al. 1992) to obtain the local value at point (r, θ, ϕ). In this way the value at any arbitrary point inside the heliosphere can be obtained.

In order to obtain the drift speed, we need a numerical method to calculate the curl of a vector field. In this study, following definition of the curl in the spherical coordinate system, the finite difference method is used to calculate the gradient of each vector component (Press et al. 1992). Figure 6 shows a comparison made between our numerical result and the analytical result of the drift speed using Parker's magnetic field (Jokipii et al. 1977). Our numerical method is consistent with the analytical result. As for the current sheet drift, cosmic-ray particles will experience a very large drift velocity because of the sudden change in the direction of the magnetic field. An analytical calculation of the drift speed using a sharply reversed magnetic field across the current sheet yields a δ-function. A more accurate simulation of particle drift using realistic particle speed and trajectory still yields a peak drift speed of one-sixth of the particle speed at the center of the current sheet (Burger & Potgieter 1989) and an enhanced drift speed appears within a four particle gyroradius from the current sheet. In our numerical calculation, the drift speed is not that large but the profile of the enhanced current sheet drift is much thicker because of a limited resolution of our MHD grid. It turns out that as long as the magnitude of the current sheet drift times the thickness is conserved, the effect of current sheet drift on cosmic-ray modulation remains the same. An advantage of this thickened current sheet drift is that we can still use a relatively large integration step size of ∼0.1 AU and do not miss the chance to see a possible large current sheet drift as cosmic-ray particles approach the current sheet. With the help of our 512 node high-performance computer, the numerical integration is still stable even with these large current sheet drift velocities. Thus, the classical form of drift velocity with a straightforward numerical calculation of the current sheet drift is still utilized.

Zoom In Zoom Out Reset image size

However, we encounter difficulty with a finite MHD grid resolution in the calculation of DSA using the TS. The grid size of our MHD model is roughly 1 AU at the radial distance of the TS. With numerical dissipation, the TS looks unrealistically wide at about 4 AU. In contrast, the observed thickness of the TS is only a few hundred thousand km (Richardson et al. 2008). As pointed out by Zhang (2007), the ratio of the diffusion coefficient to the product of the upstream plasma speed and the shock thickness affect the resulting accelerated particle spectrum. Only when this ratio is significantly above 1 can the DSA results be approached. With a shock ramp 4 AU thick, the ratio for most cosmic-ray particles falls below 1 and the effect of DSA is underestimated. In order to avoid any possible effect from a smoothed shock, in this study we use a hyperbolic function to sharpen the TS. Following Equation (9), the solar wind speed profile can be sharpened,

$$\begin{equation} V_{{\rm sw}}= \frac{-\tanh ((r-r_{{\rm ts}})/w) + 1 }{2}\times (V_{{\rm up}}-V_{{\rm down}}) + V_{{\rm down}}, \end{equation} \tag{ 9 }$$

where r_ts is the location of the TS and V_up and V_down are the solar wind speed upstream and downstream of the TS, respectively. The coefficient w determines the width of the TS. If we choose w = 0.1, it produces a shock profile with a ramp width of ∼0.2 AU thick. Figure 7 shows the sharpened TS and original MHD shock profile. Similarly, the magnetic field profile near the shock is also sharpened.

Zoom In Zoom Out Reset image size

Figure 7 shows some radial profiles of cosmic-ray flux for a number of different energy levels (1 GeV, 800 MeV, 600 MeV, 400 MeV, 300 MeV, 200 MeV) in the direction pointing to 55°N heliolatitude and 0° heliolongitude. Figure 8 shows the flux along the upwind direction, or exactly to the nose of the heliosphere, for a number of energy levels (1 GeV, 800 MeV, 600 MeV, 300 MeV, 200 MeV). The TS has been sharpened to less than 0.2 AU so that the DSA is not underestimated. In general, our simulations illustrate the basic trend of the cosmic-ray radial flux. It increases from the inner heliospheric region to the ISM until it reaches its corresponding LIS flux value. Both the heliosheath and the TS contribute to the modulation. (There will be a detailed discussion about the effect of DSA caused by the TS in the next section.) Both Figures 7 and 8 show that the heliosheath region contributes to the modulation. Thus, the LIS can only be obtained after the Voyager spacecraft (or future spacecraft) fully enter into the LISM.

Zoom In Zoom Out Reset image size

Based on these two simulations, we can calculate the radial gradient of cosmic-ray flux given by the following formula:

$$\begin{equation} G_{\theta }(r) = \frac{d \ln J}{dr} \times 100\%. \end{equation} \tag{ 11 }$$

In the supersonic solar wind region, the radial gradients for 200 MeV, 300 MeV, and 600 MeV along the upwind direction are 2.03% (AU⁻¹), 2.01% (AU⁻¹), and 1.83% (AU⁻¹), respectively. Using Ulysses data, Heber et al. (2002) estimated that for protons with energy greater than 250 MeV, the radial gradient was 2.2% ± 0.6% (AU⁻¹) before early 1998 and 3.5% (AU⁻¹) thereafter. Our simulation result agrees with the value before early 1998, a solar minimum year. As previously stated, this is because the MHD data we have used have a flat current sheet, representing the solar minimum condition in a positive polarity period.

The radial gradients change after the TS in the local direction. These gradients are much larger than the values in the supersonic solar wind region before becoming smaller in the heliosheath (∼50 AU in Figure 7 and ∼20 AU in Figure 8 downstream of the TS crossing). These large radial gradients occur in transition regions where magnetic field lines can intersect the TS at different longitudes. In some manner, the effect of a non-spherical shock appears to contribute quite a bit to the cosmic-ray radial gradient in this region around the TS. Once deep in the heliosheath, the radial gradients settle down to their own levels according to particle energy. At high energies, the radial gradients can become smaller than in the supersonic solar wind.

In comparison, Figure 9 shows the radial flux of 600 MeV protons along the upwind direction using two sets of MHD data (one from Linde and the other from Pogorelov). Although we use the same form of LIS and diffusion coefficients, these two sets of MHD data give different radial profiles of cosmic-ray flux. Specifically, the Galactic cosmic-ray flux inside the TS decreases significantly using the new MHD data from Pogorelov et al. (2008) with a more correct treatment of interstellar neutral atoms. Our new results here are more in line with the simulation by Florinski & Pogorelov (2009), but there are still some differences.

Zoom In Zoom Out Reset image size

After considering the momentum and energy transfer between these interstellar neutral atoms and solar wind protons, the region inside the TS expands significantly (see Table 1). The increased cosmic-ray modulation level probably comes from the geometric extension of the TS. Cosmic-ray particles encounter the supersonic solar wind in a much larger region. Due to the convective effect, solar wind expels the incoming stream of cosmic-ray particles. In accordance with this, the modulation level increases and the cosmic-ray intensity decreases. Another factor that can affect cosmic-ray modulation is the shape of the TS. The TS became less asymmetric in our new MHD model which treats interstellar neutrals. A previous study by Caballero-Lopez et al. (2004) also found that the shape of the TS affects cosmic-ray transport.

In comparison with the work done by Florinski & Pogorelov (2009), our simulated radial profiles of cosmic-ray flux begin to show an increase in radial gradient roughly at the location of the TS. This is consistent with the previous scenario. Cosmic-ray particles gain energy by crossing the TS because of the DSA. Since the cosmic-ray intensity in the lower energy region of the LIS is higher, the cosmic-ray intensity value around the TS region will show a jump in radial gradient (Kota & Jokipii 1983).

To illustrate the effect of the TS and to explain our results more clearly, we simulate cosmic-ray transport in a heliosphere with a spherical TS so that the effect is concentrated at one radial distance. Figure 10 shows how the TS affects the radial flux using our previous analytical cosmic-ray modulation code, which adapts the standard Parker's IMF model (Zhang 1999). The shock is treated as a step function, where the solar wind speed suddenly decreases. The skew Brownian motion method proposed by Zhang (2000) was also used to calculate DSA. Clearly, after the acceleration effect of the TS is removed, the radial flux decreased. The largest flux decrease happens at the TS location, and the apparent radial gradient change is seen at the TS. In addition, since a spherical shock with a radius of 92 AU is used, the change in the cosmic-ray flux gradient is exactly located at 92 AU. These results strongly support the fact that DSA affects cosmic-ray transport inside the heliosphere and the flux value around the shock increases as a result of acceleration.

Zoom In Zoom Out Reset image size

In our current MHD heliospheric model, the shock is not an abrupt boundary. The hyperbolic function tanh is used to sharpen the TS (Equation (9)). The shock width can also be changed to see how it affects the radial flux. Figure 11 shows how the radial flux for 300 MeV cosmic-ray particles varies as shocks with different widths are used. (In order to highlight the difference of the radial flux caused by different treatment of the TS, the log scale is used in this plot. In addition, since the interstellar intensity for 300 MeV cosmic-ray particles is relatively higher, Figure 11 looks different from Figure 10.) It has been shown that the skew Brownian motion method gives the largest flux value since the shock is treated as a step function and the width of the TS is infinitely thin. When the width of the TS increases from 0.2 AU to 2 AU, the flux value also decreases. The cosmic-ray intensity level is the lowest for the case without acceleration by the TS.

Zoom In Zoom Out Reset image size

Previously, Jokipii et al. (1993) and Caballero-Lopez et al. (2004) also investigated the TS effect on cosmic-ray modulation. They saw change in radial gradient across the TS. Their methods of calculation are different from ours, but we reach the same conclusion.

While the change in the radial flux gradient (Figure 10) appears at just one location using the analytical spherical shock model, the flux curve (Figures 7 and 8) obtained from the cosmic-ray transport code using the MHD plasma background shows a wide transition region at the TS. This is probably due to the asymmetric shape of the TS in the MHD data. Since the location of the TS varies on different field lines, the effect of shock acceleration can propagate to some distance away.

As pointed out by Jokipii et al. (1979), the drift effect causes the latitude difference in the radial flux. During the solar minimum period, when the polarity of the Sun's magnetic field is positive, most cosmic-ray particles enter inward from the heliospheric polar region and outward near the ecliptic plane region. The intensity of the cosmic-ray particles in the ecliptic plane is lower than in the polar region. The overall latitudinal gradient is positive. Our simulation result has confirmed this scenario. Figure 12 shows the simulated cosmic-ray radial flux for different latitudes toward the nose direction. In the inner heliosphere, the flux increases as the heliolatitude increases. However, the situation reverses in the heliosheath. The flux in the lower heliolatitude is higher. This is probably attributed to the fact that diffusion is more important than drift inside the heliosheath. Another factor that can cause this result is the northward flip of the current sheet (see Figure 1). As the current sheet flips northward, most of our simulated time-backward particles exit the outer boundary near the polar region.

Zoom In Zoom Out Reset image size

The radial flux in the southern heliosphere also changes with latitude. Figure 13 shows how the radial flux varies with latitude in the southern heliosphere. The simulation shows nearly the same trend as in the northern heliosphere. In the inner heliosphere, the radial flux increases as the latitude increases, and the latitudinal gradient is positive. The situation reverses in the outer heliosphere, where the flux near the ecliptic plane region has a higher value.

Zoom In Zoom Out Reset image size

Figure 14 shows a comparison between our calculated cosmic-ray energy spectra at several locations. At the same radial distance (either 105 AU or 65 AU), the spectra at 55°S and 55°N are almost identical, but the spectra at the flank (90° from the nose in the longitudinal direction) are significantly different from those at the nose. This tells us that the variation of cosmic-ray flux strongly depends on longitude. This conclusion was not obtained by previous studies with a two-dimensional symmetric model because the longitude effect was neglected (Haasbroek & Potgieter 1998).

Zoom In Zoom Out Reset image size

The actual longitude gradient varies with a particle's energy level. Table 2 shows the longitude gradient (the percentage of cosmic-ray flux difference per degree from the nose to flank) for a number of energies. All these gradient values are calculated based on the simulation result shown in Figure 14.

Table 2. Flux Values for Different Energy Levels at 105 AU N and (105 AU, 90°, 35°) and Associated Longitude Gradient

Energy	(105 AU, 55°, 0°) Flux	(105 AU, 55°, 90°) Flux	Longitude Gradient $G_{\theta }=\frac{d \ln J}{d\theta } \times 100 \%$
(GeV)
0.2	0.08837	0.13557	0.47%
0.25	0.08709	0.13804	0.50%
0.3	0.08196	0.13757	0.56%
0.35	0.07910	0.13803	0.60%
0.4	0.07714	0.13748	0.62%
0.45	0.07192	0.13258	0.66%

Note. 105 AU N = (105 AU, 55°, 0° nose); (105 AU, 90°, 35°) = (105 AU, 55°, 90° flank).

Download table as:  ASCIITypeset image

Since a global asymmetry exists in the north–south dimensions of the heliosphere, cosmic-ray transport behaves differently in these two hemispheres. Figure 15 illustrates how the north–south asymmetry affects the radial flux. Although the asymmetric effect exists both inside and outside the TS, it varies with latitude. Near the ecliptic plane (10°), the radial flux in these two hemispheres is nearly the same for r < 90 AU or inside the TS. The flux varies as the latitude increases. Up to 60°, the flux difference along the same latitude in the northern and southern heliosphere already becomes significant, even inside the TS.

Zoom In Zoom Out Reset image size

This north–south asymmetrical effect has also been explored by Ngobeni & Potgieter (2011). They compared the radial flux along the Voyager latitudes of 35°N and 35°S and found that for E > 200 MeV cosmic rays, the effect of north–south asymmetry is insignificant.

Figure 16 shows the simulated energy spectra at (105 AU, 55°N, 0° nose longitude) and (65 AU, 55°N, 0° nose longitude) in the northern heliosphere. For comparison, the energy spectra in the southern heliosphere, i.e., the energy spectra at locations (105 AU, 55°S, 0° nose longitude) and (65 AU, 55°S, 0° nose longitude) are also shown. This demonstrates that in the region inside the TS, there is little difference between the spectra at the northern and southern heliosphere. As the simulation location is moved outward of the TS, these spectra show greater difference, particularly in the lower energy end. The north–south asymmetry of the heliosphere seems to have a more significant effect on lower energy cosmic-ray particle transport.

Zoom In Zoom Out Reset image size

Since the shape of the LIS specifies the boundary condition in the numerical simulation, the cosmic-ray intensity value inside the heliosphere depends on the choice of input for the LIS. As of now, there is still no direct observational approach for obtaining the LIS. Only by using a diffusion model for cosmic-ray propagation in the Galaxy can the LIS form be indirectly calculated (Webber & Higbie 2009). In the following paragraphs, simulated energy spectra using different LIS forms will be shown to illustrate the LIS effects on cosmic-ray modulation.

Figure 17 demonstrates the energy spectra at 105 AU in the northern heliosphere 55° N and southern heliosphere 55° S using different LIS forms. One is from Equation (10), and the other is from the recent result of LIS obtained by Webber & Higbie (2009), which can be fitted by the following function:

$$\begin{equation} f(p) \propto \frac{p^{7.941}}{(0.1121^{0.9564}+p^{0.9564})^{13.4095}}. \end{equation} \tag{ 12 }$$

The simulation result shows that using a new LIS form does not change the previous results, once we normalized Zhang's LIS to the LIS of Webber & Higbie (2009) at the high-energy end. North–south asymmetry still does not show up at the higher energy end for the simulated spectrum plot; the flux in the southern heliosphere is slightly larger than the flux in the northern heliosphere only for the E < 300 MeV cosmic rays.

Zoom In Zoom Out Reset image size

Similarly, Figure 18 shows the energy spectra at 65 AU in the northern and southern heliosphere using different LIS forms. Regardless of whether Zhang's or Webber's LIS is used, the north–south asymmetry effect is small.

Zoom In Zoom Out Reset image size

Figure 19 shows the radial flux along the upwind direction using these two sets of LIS. Once the fluxes are normalized, the two curves are almost the same. For example, considering the case using the LIS from Webber, the flux level at 90 AU is 0.0576. Comparing this value to the LIS value 0.1526, the whole heliosheath and TS region contribute (0.1526 − 0.0576)/0.1526 ≈ 62.3% to the global modulation. On the other hand, using Zhang's LIS (Equation (10)), the LIS value is 0.2575, and the flux at 90 AU is 0.0982. The modulation outside the TS is (0.2575–0.0982)/0.2575 ≈ 61.9%, which is nearly the same as the previous value. Therefore, regardless of LIS form, the relative contribution from different regions is the same. In addition, the radial gradient inside the TS using the Webber–Higbie LIS is 1.87% AU⁻¹, which is similar to the radial gradient value 1.83% AU⁻¹ using Zhang's LIS.

Zoom In Zoom Out Reset image size

Since our cosmic-ray transport model is three-dimensional, a realistic simulation of the radial flux along the Voyager spacecraft's direction can be made. In our simulation, the direction of Voyager 1 is assumed along (35°N, 0° (nose)), and Voyager 2 along (30°S, 40° (nose)). Although the latitude varies along Voyager 1 and 2's real trajectories, the simulation shows that the latitude variation does not affect the simulated radial flux. Figure 20 illustrates the radial flux comparison between the real Voyager 2 trajectory and the simplified trajectory without the latitude change. Within a 5% difference, these two flux curves are nearly the same. This difference does not affect our qualitative study.

Zoom In Zoom Out Reset image size

Figure 21 shows the radial profiles of cosmic-ray flux along the Voyager 1 and Voyager 2 directions. It demonstrates obviously different radial fluxes. Because the latitude variation does not significantly affect the result, this implies that the longitudinal difference causes the radial gradient difference in the Voyager 1 and Voyager 2 directions. However, with previous two-dimensional simulations for Voyager 1 and 2, the longitudinal difference in cosmic-ray radial flux may be easily neglected (Manuel et al. 2011).

Zoom In Zoom Out Reset image size

If we assume the speeds of the Voyager spacecraft are constant, the radius and time relationship for Voyager 1 and 2 can be derived as

$$\begin{eqnarray} &&{\rm Time\;(year)} = 0.2768\times r({\rm AU}) + 1978.9 \:({\rm for}\;\; {\rm {\textit {Voyager}}}\, \textit {1})\nonumber\\ \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} &&{\rm Time \;(year)} = 0.3334\times r({\rm AU}) + 1979.9 \:({\rm for}\;\; {\rm {\textit {Voyager}}}\, \textit {2}).\nonumber\\ \end{eqnarray} \tag{ 14 }$$

Then we can predict the cosmic-ray fluxes along Voyager 1 and 2's directions as a function of time (Figure 22). The prediction shows that Voyager 1's flux is higher than Voyager 2's before 2012, which agrees with observations, since the cosmic-ray subsystem (CRS) measurement value for Voyager 2 never exceeds that of Voyager 1. However, the apparent radial gradient derived from a direct comparison of the cosmic-ray flux at Voyager 1 and Voyager 2 tends to underestimate the true radial gradient. For example, in 2012 January, Voyager 1 was about 120 AU far from the Sun, while Voyager 2 was 20 AU behind, about 98 AU from the Sun. The simulated 180 MeV cosmic-ray flux at 120 AU (Voyager 1's location) is 0.237 and 0.234 at 98 AU (Voyager 2's location). Thus, the apparent radial gradient nearly vanishes (0.01% per AU) if one uses the two Voyagers' cosmic-ray flux difference. However, based on the simulation result, the radial gradient along Voyager 1 is about 0.2% AU⁻¹ at 120 AU and about 0.3% AU⁻¹ along Voyager 2 at 98 AU (nearly 30 times larger than the apparent radial gradient). This demonstrates how longitude affects the cosmic-ray radial profile. After 2012, the flux values for the two spacecraft will become comparable, even though there is still a radial gradient in Figure 21. At some point, Voyager 2 may even see a higher cosmic-ray flux than Voyager 1, but this does not mean that the true radial gradient of the cosmic-ray flux is negative.

Zoom In Zoom Out Reset image size

Figure 23 shows the simulated proton spectrum at 2008 Voyager 1's location. The cosmic-ray measurement from Voyager 1 (Webber et al. 2008) is also shown in the same plot. In this simulation, the LIS used is from Webber & Higbie (2009). It can be seen from the plot that our simulated spectrum is roughly within the measured spectrum of the proton.

Zoom In Zoom Out Reset image size