IS THE ACCELERATION OF ANOMALOUS COSMIC RAYS AFFECTED BY THE GEOMETRY OF THE TERMINATION SHOCK?

Neutral atoms that have entered the heliosphere from the surrounding local interstellar medium (LISM) could become ionized by photoionization or through charge exchange with the solar-wind ions. The newly created ions, called pickup ions (PUIs), are singly charged, in contrast to the solar wind, which is in a higher state of ionization (Fisk et al. 1974; Bame et al. 1968). PUIs become assimilated into the solar wind and travel toward the termination shock (TS), where they are believed to be accelerated by the diffusive shock acceleration (DSA) process (Krymskii 1977; Bell 1978; Blandford & Ostriker 1978). In the plasma frame, PUIs have velocities of the order of the bulk speed of the solar wind, which are much larger than the typical thermal speeds of core solar-wind ions. Because of their greater speeds and larger Larmor radii, PUIs are preferred to solar-wind ions by the shock acceleration process (Fisk et al. 1974).

PUIs accelerated to high (a few MeV n⁻¹ to ∼100 MeV n⁻¹) energies are known as anomalous cosmic rays (ACRs) because of their composition, reflecting the interstellar source material abundances. Well before their discovery, Jokipii (1968) suggested that cosmic rays could be accelerated at a shock-like boundary of the solar system. The detection of the ACR component was first reported by Garcia-Munoz et al. (1973). Fisk et al. (1974) proposed that the ACRs are interstellar neutral atoms that are singly ionized and accelerated inside the "solar cavity" (i.e., the heliosphere). Subsequently, Pesses et al. (1981) developed the first reasonably complete model of continuous ACR production by first-order Fermi acceleration at the solar-wind TS.

The Voyager 1 (V1) and Voyager 2 (V2) space probes, launched in 1977, have at the time of this writing reached the distances of 125 AU and 103 AU from the Sun, respectively. The spacecraft are instrumented to detect ACRs over most of their energy range. ACR intensities were steadily increasing prior to the shock crossing, as was widely expected (e.g., Cummings et al. 2002). However, when V1 crossed the TS in December of 2004 at 94 AU, the measured spectrum was not consistent with the single power law that what was predicted based on steady, first-order Fermi acceleration (Stone et al. 2005). Lower ACR energy intensities were much lower than the predicted source intensity, indicating that there was substantial residual modulation of the ACRs. Subsequently, ACR spectra began to unfold in the heliosheath. Between 2004 and 2007, the intensity of 12–22 MeV n⁻¹ ACR He at V1 increased by a factor of ∼21. When V2 crossed the TS in September of 2007 at 84 AU, the observed spectra were qualitatively similar, although the intermediate-energy ACR intensities were higher owing to diminished solar modulation associated with declining solar activity from 2004 to 2007 (Stone et al. 2008). In 2012 August, V1 measured a rapid decrease and eventual disappearance of heliospheric ions, including the ACRs (Krimigis et al. 2013; Stone et al. 2013; Webber & McDonald 2013). In 2013 September, this event was interpreted as a possible crossing of the heliospheric boundary, where these particles escape into the LISM (Gurnett et al. 2013).

Observations at the TS led to a broad range of alternative theoretical models of ACR acceleration. McComas & Schwadron (2006) proposed that a blunt TS accelerates ACRs at the flanks and the tail of the shock. They suggested that the flanks and the tail are connected through the heliosheath region by magnetic field lines that enable particles' continuous access to the acceleration region. Schwadron et al. (2008) have developed a theoretical model based on this concept. Similarly, the ACR simulation of Kóta & Jokipii (2008) predicted higher ACR intensities in the direction of the heliotail. Subsequently, Kóta & Jokipii (2008) were able to reproduce the unfolding of the ACR spectra in the heliosheath. In other models, ACRs were not accelerated at the TS at all. Fisk & Gloeckler (2009) suggested that ACRs are accelerated by random plasma compressions in a region near the heliopause and subsequently diffuse back into the heliosheath. Giacalone et al. (2012) argued that their acceleration mechanism is diffusive compression acceleration, rather than stochastic acceleration, although Zhang (2006) and Strauss et al. (2010) previously examined the possibility of stochastic, or second-order Fermi, acceleration of ACRs, which is a much slower process compared to shock acceleration. Another mechanism, proposed by Drake et al. (2010), accelerates ACRs in the heliosheath via reflection in contracting magnetic islands formed as a result of magnetic reconnection. Finally, Guo et al. (2010) studied the effects of large-scale sinusoidal magnetic field variations on DSA, showing that shock acceleration is suppressed at locations where the field line connection points to the shock move apart. Conversely, the regions where the connection points approach each other act as "hotspots" that accelerate particles to high energies.

Recently, Jokipii (2012) reexamined the observations made by V1 nearly 20 years ago (McDonald et al. 2000) and pointed out that these observations contain crucial information about the source of ACRs. The data showed a prominent decrease in cosmic-ray intensity after a global merged interaction region (GMIR) passed V1 at 46.1 AU. While ACRs recovered within a year, it took galactic cosmic rays (GCRs) much longer to reach their pre-GMIR levels. The two-dimensional simulation of a propagating disturbance by Jokipii & Kota (2001) was also consistent with this result. By comparing the recovery of ACRs and GCRs, McDonald et al. (2000) calculated the source of ACRs to lie 88.5 ± 7 AU from the Sun, which is consistent with the location of the TS during both V1 and V2 crossings. This result, it appears, would rule out most of the alternative acceleration mechanisms operating in the heliosheath and near the heliopause.

In this paper we investigate the acceleration of ACRs at the TS using a three-dimensional transport model that includes drift motion. We employ two axisymmetric (with respect to the interstellar flow direction) plasma models (1) with a spherical TS and (2) with a non-spherical TS, based on realistic heliospheric conditions. A combination of several simple flow geometries was used to analytically compute the plasma velocity field. Particle simulations were performed for He⁺ ions, where ACR features clearly stand out from the GCR background (the singly ionized He is unique for ACRs). We have adapted the cosmic-ray transport code of Florinski & Pogorelov (2009), based on the stochastic integration technique along the phase-space trajectories of the Parker equation, for ACR acceleration and transport. The model assumes a steady background of pre-accelerated PUIs below a fixed momentum p_min. Starting from a larger momentum, the transport equation is integrated along the phase-space trajectories backward in time until p_min is reached. In a way, our model is related to the two-dimensional models of Kóta & Jokipii (2008) and Kóta (2012). However, with our three-dimensional approach we are able to investigate ACR behavior at high latitudes and study the effects of the gradient and curvature drifts on the acceleration process.

The remainder of this paper is structured as follows. Section 2 describes the plasma and magnetic field backgrounds used in our model. The energetic charged particle transport model is presented in Section 3. Section 4 reports the results of our simulations and the analysis of the ACR energy spectra and their radial intensity profiles. Finally, Section 5 summarizes the results, discusses the limitations, and suggests some future improvements to the current model.

The solar wind is both supersonic and super-Alfvénic upstream of the TS. By the time the solar wind reaches the TS, it becomes entrained with PUIs, which are accelerated and heated by the pickup process that transfers the bulk flow energy from the solar wind to the PUIs (e.g., Richardson et al. 2008). As a result, the solar wind decelerates by ∼20% before reaching the TS (Richardson & Stone 2009). ACRs also gain energy at the TS from the upstream plasma flow. The solar wind is decelerated, compressed, and heated at the shock. The region that lies between the TS and the LISM is known as the heliosheath. The heliosheath is shaped like a blunt cone in the direction opposite to the LISM flow vector (the "upwind" or "nose" direction) and has a comet-like tail in the opposite direction. The boundary of the heliosheath, separating the heliosphere from the LISM, is called the heliopause.

Two coordinate systems are commonly used to describe the heliospheric interface, one attached to the solar rotation axis and the other to the interstellar flow vector. For convenience, both systems are used in our model. The background plasma flow is derived in the coordinate frame where the z-axis points in the upwind direction. All the equations in this paper are written in this coordinate frame, where they have a relatively simple form. Subsequently, a coordinate transformation is performed into the frame where the z-axis is aligned with the solar rotation axis. This is essentially the solar equatorial coordinate system that the reader is more familiar with. This frame is used for all the figures in this paper. The plasma flow model used here is axisymmetric about the direction of the interstellar wind, permitting the use of scalar potentials or stream functions which greatly simplifies the derivations. The fully three-dimensional magnetic field is treated in a kinematic approximation, ignoring the magnetic stresses on the plasma.

In constructing analytic models of the heliosphere, assumptions of flow incompressibility and/or potentiality are commonly used. A flow is said to be incompressible if the density of a fluid parcel is not affected by changes in the pressure (e.g., Batchelor 2000) and the flow velocity has a zero divergence ($\nabla \,{\boldsymbol \cdot}\, \mathbf {u} = 0$). A flow is called irrotational (or potential) if the curl of its velocity, known as the vorticity, is zero ($\boldsymbol {\omega }=\nabla \times \mathbf {u}=0$). Parker (1963) and Suess & Nerney (1990, 1991) used the potential flow approximation to compute the plasma velocity field in the subsonic region (i.e., the heliosheath), whereas Khabibrakhmanov & Summers (1996) used the stream function approach (see below). These authors treated the heliosheath plasma as both incompressible and irrotational. Parker (1963) estimated that the plasma pressure behind the TS is only 15% smaller than the interstellar gas pressure. Under these conditions, the density variations in the heliosheath are small, and the flow can be approximated as incompressible. We argue that the incompressibilty condition is more important in the context of ACR acceleration because a compressible heliosheath flow would accelerate the particles and compete with shock acceleration, making it impossible to investigate the effect of the shock geometry alone. Since the power-law slope, τ, of the ACR momentum spectra depends on the shock compression ratio s as τ = 3s/(s − 1), satisfying the Rankine–Hugoniot (RH) conditions at the TS is of crucial importance in our problem. Note that none of the above-mentioned models fully satisfy the RH conditions at the TS.

In the supersonic solar-wind region we assume a uniform, constant radial velocity u_u = 450 km s⁻¹. Generally, during solar-minimum conditions the slow solar wind is prevalent at low latitudes, and the fast solar wind is prevalent over the polar regions. However, this is not critical in our simulation because we are primarily interested in how the shock geometry affects the ACR spectra. The heliosheath flow velocity may be expressed as a curl of a Stokes' stream function $\psi \hat{\mathbf {e}}_\phi$, where $\hat{\mathbf {e}}_\phi$ is the unit vector along the azimuthal direction. Such a flow automatically satisfies the incompressibility condition by construction. In spherical coordinates, the velocity is expressed as

$$\begin{equation} u_r=\frac{1}{r^2\sin \theta }\frac{\partial \psi }{\partial \theta },\quad u_\theta =-\frac{1}{r\sin \theta }\frac{\partial \psi }{\partial r}, \end{equation} \tag{ 1 }$$

where ψ is Stokes' stream function. By taking the curl of Equation (1), one obtains the equation for the stream function,

$$\begin{equation} \frac{1}{r\sin \theta }\frac{\partial ^2\psi }{\partial r^2} +\frac{1}{r^3}\frac{\partial }{\partial \theta }\left(\frac{1}{\sin \theta } \frac{\partial \psi }{\partial \theta }\right)=-\omega , \end{equation} \tag{ 2 }$$

where ω is in the $\hat{\phi }$ direction. The equation for vorticity can be derived by taking the curl of the Euler equation,

$$\begin{equation} \mathbf {u}\,{\boldsymbol \cdot}\, \nabla \left(\frac{\omega }{r\sin \theta }\right)=0. \end{equation} \tag{ 3 }$$

This equation is difficult, if not impossible, to satisfy unless a finite (small) vorticity is present in the solar wind. However, the more important incompressibility condition is always satisfied by the stream function approach.

The expression we use for the stream function is a combination of several simple flow patterns (see, e.g., Batchelor 2000). We chose

$$\begin{eqnarray} \psi &=& {-}u_0 R_0^2\cos \theta -\frac{u_{\infty 1} r^2}{2}\sin ^2\theta -\frac{u_{\infty 2} R_0^3}{r}\sin ^2\theta\nonumber\\ &&+u_{\infty 3} R_0 r\sin ^2\theta +\frac{u_{\infty 4} R_0^4}{2r^2}\sin ^2\theta , \end{eqnarray} \tag{ 4 }$$

where R₀ is a constant with the dimension of length. Substituting Equation (4) into Equation (2) gives the vorticity

$$\begin{equation} \omega = \frac{2\sin \theta }{r} \left[ 2 u_{\infty 3}\frac{R_0}{r} - 2 u_{\infty 4} \left(\frac{R_0}{r}\right)^4\right]. \end{equation} \tag{ 5 }$$

In Equation (4), the first term is a point source of strength $u_0 R_0^2$ placed at the origin, the second term describes a uniform flow with speed $\mathbf {u}=u_{\infty 1}\hat{\mathbf {e}}_z$, the third term is a dipole flow of strength $u_{\infty 2} R_0^3 \hat{\mathbf {e}}_z$ at the origin, the fourth term is a Stokes flow around a sphere of radius 4R₀/3 moving with velocity $\mathbf {u}=u_{\infty 3}\hat{\mathbf {e}}_z$, and the fifth term is a spherical vortex of strength $u_{\infty 4} R_0^4$ centered at the origin.

We begin with the spherical TS model. A non-spherical TS will be introduced in the subsequent section. The spherical TS model serves as a reference against which we compare the results for other shock geometries.

2.1. Spherical Termination Shock

Consider a spherical shock of radius R₀ and constant compression ratio s. In order to satisfy the gas-dynamic RH conditions at a spherical TS, the tangential flow velocity must be zero immediately downstream of the TS. Using the above generalized stream function, this condition is satisfied if u_∞1 = u_∞2 = u_∞3 = u_∞4(= u_∞). The plasma velocities are obtained by taking the curl of Equation (4),

$$\begin{eqnarray} \nonumber u_r &=& u_0 \left(\frac{R_0}{r}\right)^2 + 2 u_{\infty }\cos \theta \left[ \frac{R_0}{r} - \frac{1}{2} - \left(\frac{R_0}{r}\right)^3 + \frac{1}{2} \left(\frac{R_0}{r}\right)^4\right], \\ u_{\theta } &=& {-}u_{\infty } \sin \theta \left[ \frac{R_0}{r} - 1 + \left(\frac{R_0}{r}\right)^3 -\ \left(\frac{R_0}{r}\right)^4\right], \\ \nonumber u_{\phi } &=& 0. \end{eqnarray} \tag{ 6 }$$

We remind the reader that these velocities are written in the spherical coordinate system with the z-axis pointing in the upwind direction. To match the locations of the heliospheric boundaries, we chose u₀ = 150 km s⁻¹ (velocity downstream of the shock), u_∞ = 400 km s⁻¹, and R₀ = 90 AU. u_∞ may be interpreted as the effective velocity at infinity (i.e., in the LISM) if the densities in the heliosheath and the LISM were the same. In reality, the interstellar density is much larger, and the velocity correspondingly smaller (e.g., Pogorelov et al. 2004), but this is of no relevance here because the ACR transport model is confined to the heliosheath region (see below). The upstream solar-wind velocity u_u = 450 km s⁻¹, which yields the shock compression ratio s = 3. Figure 1 (left panel) shows the velocity magnitude (color) overplotted with the plasma flow streamlines in the ecliptic plane (or any other plane containing the z-axis). The heliosphere has a familiar comet-like shape with the stagnation point (the nose of the heliopause) at a distance of 157 AU from the Sun. The flow velocity given by Equation (6) has non-zero vorticity, which explains the convergence of the streamlines in the tail for this geometry.

Zoom In Zoom Out Reset image size

2.2. Non-spherical Shock

For a non-spherical TS, the plasma velocity immediately downstream of the shock is no longer radial. Consider a surface of revolution (about the z-axis) given by the function R = R(θ). The tangent vector to the surface is

$$\begin{equation} \hat{\boldsymbol {\tau }}=\frac{dR \hat{\mathbf {e}}_R+Rd\theta \hat{\mathbf {e}}_\theta }{\sqrt{dR^2+R^2d\theta ^2}} =\frac{\dot{R}\hat{\mathbf {e}}_R+R\hat{\mathbf {e}}_\theta }{\sqrt{\dot{R}^2+R^2}}, \end{equation} \tag{ 7 }$$

where $\dot{R}=dR/d\theta$. We can likewise express the normal vector as

$$\begin{equation} \hat{\mathbf {n}}=\frac{R\hat{\mathbf {e}}_R-\dot{R}\hat{\mathbf {e}}_\theta }{\sqrt{\dot{R}^2+R^2}}. \end{equation} \tag{ 8 }$$

The geometry of the shock and the unit vectors are schematically shown in Figure 2.

Zoom In Zoom Out Reset image size

Mass flux conservation for a stationary shock wave given by R = R(θ), with compression ratio s(θ) and upstream radial speed u_u, requires that

$$\begin{equation} \mathbf {u}\,{\boldsymbol \cdot}\, \hat{\mathbf {n}} =\frac{1}{s}u_u\hat{\mathbf {e}}_r\,{\boldsymbol \cdot}\, \hat{\mathbf {n}}, \end{equation} \tag{ 9 }$$

$$\begin{equation} \mathbf {u}\,{\boldsymbol \cdot}\, \hat{\boldsymbol {\tau }} =u_u\hat{\mathbf {e}}_r\,{\boldsymbol \cdot}\, \hat{\boldsymbol {\tau }}. \end{equation} \tag{ 10 }$$

Writing this in component form, we obtain

$$\begin{equation} u_r R-u_\theta \frac{d R}{d\theta }-\frac{u_u R}{s}=0, \end{equation} \tag{ 11 }$$

$$\begin{equation} u_r\frac{d R}{d\theta }+u_\theta R-u_u\frac{d R}{d\theta }=0. \end{equation} \tag{ 12 }$$

This yields a first-order ordinary differential equation for R,

$$\begin{equation} \frac{d R}{d\theta }=\frac{u_\theta R}{u_u-u_r}, \end{equation} \tag{ 13 }$$

and an algebraic expression to find s from

$$\begin{equation} s=\frac{u_u R}{u_r R-u_\theta (dR/d\theta)}. \end{equation} \tag{ 14 }$$

Because R and dR/dθ vary with θ, it is not possible to obtain a solution with a constant compression ratio along the shock. We seek the values for the constants, u₀, u_∞1, u_∞2, u_∞3, u_∞4, that (1) yield the correct shock location at the nose, (2) produce a stagnation point on the heliopause that is near to that of the spherical TS model discussed earlier, (3) have the required elongation in the tail of the TS, and (4) yield a minimum variation of s(θ). Our solution has the nose of the TS at 90 AU and the tail at 135 AU. The stagnation point is within ∼1% to the value for the spherical TS, and the compression ratio deviates by only 10% from its spherical TS value of 3.0. As before, the velocities are obtained by taking the curl of Equation (4),

$$\begin{eqnarray} u_r &=& u_0 \left(\frac{R_0}{r}\right)^2 + 2 \cos \theta \left[u_{\infty _3} \frac{R_0}{r}\ - \frac{u_{\infty _1}}{2}\right.\nonumber\\ &&\left.- u_{\infty _2} \left(\frac{R_0}{r}\right)^3 + \frac{u_{\infty _4}}{2} \left(\frac{R_0}{r}\right)^4\right],\nonumber \\ u_{\theta } &=& {-}\sin \theta \left[ u_{\infty _3}\frac{R_0}{r} - u_{\infty _1} + u_{\infty _2} \left(\frac{R_0}{r}\right)^3 -u_{\infty _4} \left(\frac{R_0}{r}\right)^4\right],\nonumber \\ u_{\phi } &=& 0, \end{eqnarray} \tag{ 15 }$$

where the values for the coefficients are u₀ = 150.2 km s⁻¹, u_∞1 = 201.2 km s⁻¹, u_∞2 = 169.6 km s⁻¹, u_∞3 = 170.0 km s⁻¹, and u_∞4 = 200.2 km s⁻¹. Using Equation (15), Equation (13) can be integrated to obtain the shape of the TS as

$$\begin{equation} R = R_a \left(\frac{C-2B}{C-2B\cos \theta } \right) ^ {A/2B}, \end{equation} \tag{ 16 }$$

where R_a = 90 AU, A = (u_∞3 − u_∞1 + u_∞2 − u_∞4), B = (u_∞3 − u_∞1/2 − u_∞2 + u_∞4/2), and C = (u_u − u₀). A velocity magnitude plot in the ecliptic plane for the non-spherical TS is shown in Figure 1 (right panel). The stagnation point distance is almost the same as the spherical TS model. The TS shows significant elongation toward the tail, where its distance R = 135 AU. Compared to the spherical model, the tail is wider here because the vorticity is significantly smaller.

2.3. Magnetic Field

In the solar wind the standard Parker field (Parker 1958) is a valid zeroth-order approximation. More complicated analytic magnetic field configurations in the solar wind may be found in Scherer et al. (2010). For the heliosheath we calculate the field by solving the steady-state induction equation ∇ × (u × B) = 0 along the streamlines, similar to Nerney et al. (1991, 1995). Because $\nabla \,{\boldsymbol \cdot}\, \mathbf {u}= 0$ and $\nabla \,{\boldsymbol \cdot}\, \mathbf {B}= 0$, one obtains from the induction equation

$$\begin{equation} (\mathbf {u}\,{\boldsymbol \cdot}\, \nabla) \mathbf {B} = (\mathbf {B}\,{\boldsymbol \cdot}\, \nabla) \mathbf {u}, \end{equation} \tag{ 17 }$$

or, in spherical coordinates, using dl = (dr² + r²dθ²)^1/2,

$$\begin{eqnarray} &&\left(\begin{array}{@{}c@{}}\dot{B_r}\\ \dot{B_\theta }\\ \dot{B_\phi }\\ \end{array} \right) =\nonumber\\ &&\left(\begin{array}{@{}ccc@{}}\displaystyle\frac{1}{2u} \displaystyle\frac{\partial u_r}{ \partial r} & \displaystyle\frac{1}{2u r} \displaystyle\frac{\partial u_r}{ \partial \theta } & 0 \\ \displaystyle\frac{1}{2u} \left(\displaystyle\frac{\partial u_\theta }{ \partial r} - \displaystyle\frac{u_\theta }{r} \right) & \displaystyle\frac{1}{2u r}\left(\displaystyle\frac{\partial u_\theta }{ \partial \theta } + \displaystyle\frac{u_r}{r}\right) & 0 \\ 0 & 0 & \displaystyle\frac{u_r + u_\theta \cot \theta }{2ur}\end{array} \right)\nonumber\\ &&\quad\times\left(\begin{array}{@{}c@{}}B_r\\ {B_\theta }\\ {B_\phi }\\ \end{array} \right),\quad\quad \end{eqnarray} \tag{ 18 }$$

where $\dot{B_i} = dB_i / dl$ and u = (u_r² + u_θ²)^1/2. The magnetic field at a given grid point is calculated by integrating backward along the streamlines until the latter reaches the TS. Using the Parker field in the upstream solar-wind region, the value downstream of the shock is found from the RH conditions. We take advantage of the principle of superposition of linear equations (any linear combination of solutions is itself a solution). Two initial guesses are required, one for each component of the magnetic field, to solve the coupled equations for B_r and B_θ. The ϕ component is computed separately using a single initial guess. Suppose the initial position of integration along the streamline is l₀ and the final position is l_f. The exact solution for the field at the starting point can be written as

$$\begin{eqnarray} \nonumber B_r(l_0) &=& \alpha B_{r_1}(l_0) + \beta B_{r_2}(l_0),\\ B_\theta (l_0) &=& \alpha B_{\theta _1}(l_0) + \beta B_{\theta _2}(l_0),\\ \nonumber B_\phi (l_0) &=& \gamma B_{\phi _1}(l_0), \end{eqnarray} \tag{ 19 }$$

where the numerical subscripts refer to the two trial starting values. Then, using the computed values at the endpoints (l_f) in Equation (19), we obtain

$$\begin{eqnarray} \nonumber \alpha &=& (B_r B_{\theta 2} - B_{\theta } B_{r2}) / (B_{r1} B_{\theta 2} - B_{\theta 1}B_{r2}), \\ \beta &=& (B_r B_{\theta 1} - B_{\theta } B_{r1}) / (B_{r2} B_{\theta 1} - B_{\theta 2}B_{r1}),\\ \nonumber \gamma &=& B_{\phi } / B_{\phi 1}. \end{eqnarray} \tag{ 20 }$$

The magnetic field components without a subscript are taken at l_f, where they are already known. Once the coefficients α, β, and γ are found from Equation (20), the exact solution at the initial point l₀ is computed from Equation (19).

A flat heliospheric current sheet was used in our model to account for the drift effect during a solar minimum. The negative magnetic cycle (A < 0) conditions were used with the solar magnetic field pointing out of the Sun in the northern hemisphere and inward in the southern hemisphere. Note that the magnetic field is weak in the polar regions, which could result in very large transport coefficients. To control the magnitudes of the diffusion and drift coefficients (see below), another term, B_m, was added to the total magnitude of the magnetic field similar to Jokipii & Kota (1989) and Florinski et al. (2003). The total magnetic field is calculated as

$$\begin{eqnarray} &&B^2 = B_r^2 + B_{\theta }^2 + B_{\phi }^2 + B_m ^2. \end{eqnarray} \tag{ 21 }$$

Figure 3 shows the magnetic field magnitude with a single spiral field line plotted in the ecliptic plane (left panel) and a three-dimensional view of several selected field lines (right panel) for the non-spherical TS model. From the left panel one can see that the magnetic connection region, where the energetic particles can revisit the acceleration sites at the flanks and the tail along these field lines, spans about 5 AU in radial distance. McComas & Schwadron (2006) suggested a much wider connection region (20–30 AU), which is not reproduced in our model, or in other heliospheric MHD models (e.g., Pogorelov et al. 2006). Because the LISM region is not part of the simulation, the magnetic field is not computed there.

Zoom In Zoom Out Reset image size

2.4. Simulation Parameters

ACR transport is described via Parker's equation (Parker 1965), which is a partial differential equation (PDE) of a parabolic kind. Parker's equation describes the evolution of the probability density for a particle to be found at a point ($\bf {r}$, p) in phase space. Processes affecting the particle distribution include convection, energy change, diffusion, gradient, and curvature drifts. Most of the current work on cosmic-ray transport and acceleration is based on this equation. It is a good approximation for energetic particles if the scattering is faster than the macroscopic timescales, so that the particle distribution is close to isotropic (Jokipii 2012).

A process that fluctuates randomly in time is known as a stochastic process. It is a probabilistic process in the sense that even if the initial condition is known, there are several directions in which the process could evolve. A stochastic differential equation (SDE) is a differential equation driven by random increments at each time step. SDEs and their applications are discussed in, e.g., Gardiner (1983) and Jacobs (2010). Jokipii & Owens (1975) were the early adopters of the stochastic method to solve cosmic-ray transport problems in the solar wind. A similar approach based on SDEs was used later by Chalov et al. (1995), Zhang (1999), Ball et al. (2005), and Florinski & Pogorelov (2009). In our transport model we use a set of SDEs of the Ito type. In general, an Ito SDE for a stochastic process x(t) can be written as

$$\begin{equation} dx_i = a_i(x_i,t){dt} + b_{ij}(x_i,t)dW_j(t), \end{equation} \tag{ 22 }$$

where x_i is the phase-space coordinate, a_i(x_i, t) is a known vector, b_ij(x_i, t) is a known tensor, and dW_j(t) is an increment of the Wiener process (Gaussian random walk). According to Ito calculus, dt, dW_j, and (dW_j)² are the only terms that contribute to the solution of an SDE. By using Ito's rule, one can show that (dW_j)² = dt, and therefore the quadratic term cannot be neglected. Stochastic integrals in which the integrand is evaluated at the start of each interval are known as Ito stochastic integrals. A diffusion process could be approximated by an Ito SDE (Gardiner 1983), which is a reversible process called Ito diffusion. The Parker equation can be transformed into a Fokker-Planck equation that describes the time evolution of the probability density of x_i(t),

$$\begin{eqnarray} \frac{\partial }{\partial t} P(x_i,t) &=& {-}\frac{\partial }{\partial x_i} [a_i(x_i,t)P(x_i,t)]\nonumber\\ &&+ \frac{1}{2} \frac{\partial ^2}{\partial x_i \partial x_j} [D_{ij}(x_i,t)P(x_i,t)], \end{eqnarray} \tag{ 23 }$$

where P(x_i, t) is the probability density and D_ij(x_i, t) = b_ikb_kj. The same equation can be re-written as a backward Fokker–Plank equation (Gardiner 1983),

$$\begin{eqnarray} \frac{\partial }{\partial t} R(y_i,t) &=& a_i(y_i,-t)\frac{\partial }{\partial y} R(y_i,t)\nonumber\\ &&+ \frac{1}{2}\frac{\partial ^2}{\partial y_i \partial y_j} [ D_{ij}(y_i,-t) R(y_i,t)], \end{eqnarray} \tag{ 24 }$$

where R(y_i, t) = P(x_i, 0|y_i, −t), which is the probability density of x_i at t = 0 given the probability density of y_i for t < 0.

SDEs are often transformed into SDEs of the Ito type because of the difficulty in handling the properties of the Stratonovich integration rule (Zhang 1999). In our code we integrate Equation (22) backward in time similar to Zhang (1999), Ball et al. (2005), and Florinski & Pogorelov (2009). The backward-in-time integration is much more efficient compared to the forward-in-time integration because it ensures that only particles satisfying the boundary conditions are counted (Kota 1977; Pei et al. 2010). In this work stochastic integration is performed in Cartesian coordinates, whereas Pei et al. (2010) presented a method to integrate directly in spherical coordinates. In our simulation the displacement of a trajectory in the solar rest frame is computed as in Florinski & Pogorelov (2009), via a coordinate transformation into the field-aligned frame, where the diffusion tensor is diagonal.

We assume a background of pre-accelerated helium PUIs below a fixed value of momentum p_min, corresponding to the energy E_min = 0.5 MeV n⁻¹. The rate of injection is taken to be uniform along the shock. We then integrate along phase-space trajectories backward in time until p_min is reached. The heliopause and the spherical inner boundary are treated as absorbing boundaries for the energetic particles. Achieving p_min requires some sort of pre-acceleration. A number of pre-acceleration mechanisms have been proposed in the literature. For example, Zank et al. (2001) proposed that the PUIs undergo multiple reflections, an acceleration mechanism known as shock surfing, before the onset of DSA. Giacalone & Decker (2010) suggested that standard shock drift acceleration could accelerate low-energy PUIs. Kóta (2012) showed that PUI injection is more efficient in the nose region because the particles injected there have a better chance of reaching ACR energies because they had more time for acceleration at the shock.

The Parker equation may be written as

$$\begin{eqnarray} \frac{\partial f}{ \partial t} &+& \left(u_i + v_{d,i} \right) \frac{\partial f}{\partial x_i} - \frac{\partial }{\partial x_i } \left(\kappa _{ij} \frac{\partial f}{\partial x_j}\right)\nonumber\\ &&- \frac{\bf {\nabla\,{\boldsymbol \cdot}\, u}}{3} \frac{\partial f}{\partial \ln p} = 0, \end{eqnarray} \tag{ 25 }$$

where f is the pitch-angle-averaged phase-space density, u is the plasma velocity, κ is the diffusion tensor, and v_d is the drift velocity. Equation (25) must be converted to a backward Fokker–Planck equation for the SDE solver,

$$\begin{eqnarray} \frac{\partial f}{ \partial t} &+& \left(u_i + v_{d,i} - \frac{\partial \kappa _{ij}}{\partial x_j} \right) \frac{\partial f}{\partial x_i} - \kappa _{ij} \frac{\partial ^2 f}{\partial x_i \partial x_j}\nonumber\\ &&- \frac{\bf {\nabla\,{\boldsymbol \cdot}\, u}}{3} \frac{\partial f}{\partial \ln p} = 0. \end{eqnarray} \tag{ 26 }$$

The diffusion tensor in the solar frame is given by

$$\begin{eqnarray} &&\kappa _{ij} = \kappa _\perp \delta _{ij} + (\kappa _\parallel - \kappa _\perp) \frac{B_i B_j}{B^2}, \end{eqnarray} \tag{ 27 }$$

where κ_∥ and κ_⊥ are the parallel and the perpendicular diffusion coefficients, respectively. In this work κ_∥ was taken to be proportional to p^4/3, which corresponds to the non-relativistic limit in the Kolmogorov inertial turbulent range, and inversely proportional to the magnetic field:

$$\begin{eqnarray} \kappa _\parallel &=& \kappa _0 \left(\frac{p}{p_0}\right)^{4/3} \left(\frac{B_0}{B}\right), \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} \kappa _\perp &=& b \kappa _\parallel , \end{eqnarray} \tag{ 29 }$$

where κ₀ = 6.8 × 10²² cm² s⁻¹, B₀ is the radial magnetic field at 1 AU, p₀ = mc is the rest momentum, and b is a constant varied between 0.01 and 0.05. With these parameters, the parallel mean free path of a 0.5 MeV n⁻¹ He⁺ ion is about 9.5 AU just upstream of the TS. For numerical reasons, for the lowest energy particles, the diffusion length should satisfy the condition L_low = κ_rr(p_min)/u_r ≫ Δr just upstream of the shock, where Δr is the radial grid size. The high-energy cutoff in the spectrum is determined by the balance between the rate of acceleration and the rate of adiabatic cooling in the solar wind (Florinski & Jokipii 2003), which implies that the diffusion length at near-cutoff momentum p_cut is L_high = κ_rr(p_cut)/u_r ∼ R. The drift velocity for an isotropic distribution is given by

$$\begin{eqnarray} &&v_{d,i} = \frac{pcw}{3e} \epsilon _{ijk} \frac{\partial }{\partial x_j} \left(\frac{B_k}{B^2} \right), \end{eqnarray} \tag{ 30 }$$

where e is the electric charge of the particle and w is the velocity.

Stochastic numerical models are generally simpler and faster compared with traditional grid-based energetic particle transport models (Florinski & Pogorelov 2009). Another advantage of the SDE approach is that the sampled trajectories and propagation histories of the "particles" are available to deduce transport times, entry or exit points, and acceleration sites. Figure 4 shows two representative stochastic trajectories in the non-spherical model with final energies of 2 MeV n⁻¹ and 10 MeV n⁻¹. The 2 MeV n⁻¹ trajectory (red line) started in the ecliptic plane in the "nose" direction at 100 AU, whereas the 10 MeV n⁻¹ trajectory (blue line) started at a latitude of 45° in the tailward direction at 148 AU. The starting points are marked by the "+" signs and the endpoints by the filled circles. Since this is a backward simulation, the starting point is actually the position where ACRs are observed and the endpoint is where the PUIs were first injected into the DSA (so the endpoints always lie on the TS). The figure shows that particles move primarily along the field lines owing to parallel diffusion. They can also switch between field lines due to perpendicular diffusion (field line meandering).

Zoom In Zoom Out Reset image size

We have performed rigorous testing of the SDE model that was found to be sensitive to spatial resolution and time stepping. Below, the SDE model results are compared with the output from a conventional PDE solver for a simple test problem. Next, anomalous helium spectra at different locations on the shock and in the sheath are discussed for the two reference TS geometries. Finally, radial intensity profiles are computed for the heliospheric plasma background with a non-spherical TS.

4.1. Accuracy of the SDE Model

To assess the accuracy of our SDE transport code, we picked a simple problem with a spherical shock placed at 90 AU with compression ratio s = 3. The background flow was strictly radial, with u_r = const upstream and u_r ∼ r⁻² (i.e., incompressible) downstream, between the shock and the spherical outer boundary. For simplicity the magnetic field was taken to be radial as well. The diffusion coefficient was proportional to p² and independent of B. A uniform radial grid spacing Δr = 0.4 AU was used, with the reflective inner boundary placed at 0.4 AU and the outer absorbing boundary at 120 AU. In the SDE code, particle motion was restricted to the ecliptic plane and the drifts were switched off in order to eliminate any possible latitudinal effects from a polar grid. The "perpendicular" (i.e., azimuthal) diffusion coefficient was varied within a broad range, from zero to about 100 times the parallel diffusion coefficient (for the radial magnetic field, the directions of the parallel and perpendicular diffusive motions are reversed compared with the primarily azimuthal Parker field).

Output spectra from the SDE model were compared with the results from the spherical model of Florinski & Jokipii (2003), which solves the transport equation using an alternating direction implicit (ADI) scheme in the r and p coordinates. Figure 5 compares the spectra at the shock computed using these two very different numerical techniques. In the left panel the diffusion coefficient was constant throughout the simulation domain, whereas in the right panel it dropped across the shock by the shock compression ratio. The second scenario is a more difficult test for the code (Achterberg & Schure (2011) argued that a more accurate numerical scheme may be needed if there is a decrease in diffusivity across a shock). The SDE model results match very well with the PDE model results, even for the extreme cases of perpendicular diffusion: zero perpendicular diffusion and 100 times the parallel coefficient. The spectra consist of a power-law region at low energies and a rollover at high energies caused primarily by deceleration of particles in the expanding solar wind (Florinski & Jokipii 2003), as well as by their escape through the external boundary. The power-law slope expected for a planar shock with a compression ratio of 3.0 is plotted in the figure for comparison. These results suggest that the model is accurate enough without the corrections due to Achterberg & Schure (2011).

Zoom In Zoom Out Reset image size

4.2. ACR Energy Spectra

In this section we examine the spectral features of anomalous helium ion distributions at several representative locations in the heliosphere. The ACR differential intensity J = fp² was calculated for A < 0 solar-minimum conditions by running time-backward simulations of 10 million stochastic trajectories and binning the results into 30 bins in kinetic energy per nucleon, spaced logarithmically between 0.5 MeV n⁻¹ and 100 MeV n⁻¹. In the following simulations we used κ_⊥ = 0.01κ_∥. Figure 6 was produced with the drift term set to zero, to verify the directional uniformity of the solution (that is, to show that the spectra do not depend on the azimuthal angle). The figure shows spectra at the shock in the ecliptic plane for the spherical TS heliosphere in the nose direction (red) and the tail direction (black). The power-law slope at low energies agrees well with the theoretical value. The broad exponential rollover appears to start at about 5–10 MeV n⁻¹. The efficiency of shock acceleration drops rapidly with momentum, and the very high energy spectrum is not well resolved here. Note that the two spectra are nearly identical, which is to be expected for a spherical shock in the absence of drifts.

Zoom In Zoom Out Reset image size

Next we examine the effects of gradient and curvature drifts within the framework of the same spherical-shock model heliosphere. Figure 7 compares spectra for the spherical TS in the ecliptic plane and the polar region with and without the drift term. In the absence of drifts spectra are virtually identical, similar to Figure 6. With drifts switched on, latitudinal differentiation of the spectra becomes evident. In this sense our results agree with the standard drift paradigm, which predicts the intensity maximum (for positively charged ions) to be near the ecliptic plane for A < 0 solar-minimum conditions (e.g., Jokipii 1986; Steenkamp & Moraal 1995; Florinski et al. 2004). During the negative magnetic polarity cycles, positively charged particles drift along the TS from high to low latitudes. Because the ions are further accelerated along the way, one observes harder spectra in the ecliptic plane compared to a softer spectrum at higher latitudes. As a result, ACR helium intensities at the TS at ∼10 MeV n⁻¹ may differ by an order of magnitude between high and low latitudes.

Zoom In Zoom Out Reset image size

We now consider the non-spherical TS, again comparing the simulations with and without the drift term. The left panel of Figure 8 shows the energy spectra at different locations along the shock in the ecliptic plane, with drift effects switched off. The spectrum is much softer in the nose region. One can see a moderate increase in ACR intensities toward the flanks and a significant increase in the tail direction. This implies that the tailward part of the TS is the preferred ACR acceleration site. This result can be explained using the concept proposed by Guo et al. (2010). The authors suggested that the particles are trapped and accelerated at places where the magnetic field line connection points to the shock converge toward each other. It is evident from Figure 3 (left panel) that field line connection points are moving away from each other in the upwind (nose) direction. In contrast, in the direction of the heliotail, the connection points are approaching each other. Here, the particles are trapped between the field line connection points to the shock that act as magnetic mirrors. As the field lines are convected downstream, the trapped particles are accelerated in the region where field line connection points are moving toward each other. In the tail region, one can see a "hump" at ∼10 MeV n⁻¹, and the cutoff appears to be shifted to a higher energy.

Zoom In Zoom Out Reset image size

Figure 8 (right panel) shows the same spectra with drift motions taken into account. These results appear qualitatively similar to those without drifts but contain an enhancement in intensities, and the cutoff in the tail region extends to ∼50 MeV n⁻¹. Evidently, the enhancement is due to the drifting of particles along the TS from the polar regions toward the ecliptic plane (see also Jokipii 1986; Steenkamp & Moraal 1995; Florinski et al. 2004). In general, spectral humps are produced by a combination of drift effects, spherical geometry effects, and modified shock effects (Florinski & Jokipii 2003). The above result suggests that the unfolding of spectra observed by V1 and V2 may have been affected by the drift. As the amount of solar modulation decreased with declining solar activity, drifts started to play a more prominent role, transporting more particles to lower latitudes, closer to the Voyagers' trajectories.

Figure 9 shows the intensities of 10 MeV n⁻¹ ACR He ions along the surface of the TS in the ecliptic plane. The figure clearly shows that the efficiency of acceleration by the TS increases monotonically from nose to tail. 10 MeV n⁻¹ intensities differ by a factor of up to 30 between the two extremes. This result is again consistent with the earlier work of Kóta & Jokipii (2008).

Zoom In Zoom Out Reset image size

In Figure 10 we show spectra of anomalous He in the heliosheath along the trajectory of V1 for the case without drifts (left panel) and with drifts included (right panel). Because the Voyagers are moving approximately in the upwind (nose) direction, the spectra are similar to those in the nose direction shown in Figure 8. The spectra unfold in the heliosheath with increasing distance from the TS. We saw earlier that the particles are accelerated more efficiently in the tail region of the TS. The unfolding appears to be a result of Voyager sampling field lines connected to progressively more tailward portions of the TS, given the small perpendicular diffusion coefficient used here. The unfolding effect is qualitatively similar to the observed unfolding of the spectra by the Voyagers, but the magnitude of the effect is much smaller in the simulations. As mentioned earlier, Voyager spectra were measured over the complete declining phase of the solar cycle, with drift effects becoming progressively more important on approach to the solar minimum. In effect, one has to "combine" the two panels in Figure 10 to understand the evolution of the spectra in the heliosheath with both time and distance. A more accurate treatment of the problem would include a time-varied tilt of the current sheet to account for the drift effects during different phases of the solar cycle. Since such a model is outside the main scope of this work, we defer it for a later publication.

Zoom In Zoom Out Reset image size

4.3. ACR Radial Profiles

In order to understand the spatial variations of ACRs in the heliosheath, we computed the intensities at several representative energies as a function of heliocentric distance along the trajectory of V1, which is close to the nose direction. The simulations were performed with 1–10 million (depending on energy) stochastic trajectories, and the results were binned into 20 uniformly spaced radial intervals ranging from 70 AU to 170 AU. Figure 11 (left panel) shows the radial intensity profiles for 0.5 MeV n⁻¹, 1 MeV n⁻¹, 2 MeV n⁻¹, 5 MeV n⁻¹, and 10 MeV n⁻¹ He⁺ ions computed using κ_⊥ = 0.01κ_∥ for the non-spherical TS model, with the drift effects switched off. An exponential precursor is evident upstream of the TS, shown with a vertical dashed line. It is clear that the TS is a source of particles, and their behavior in the precursor is a well-known modulation effect (a balance between inward diffusion and outward convection with the plasma flow). The precursor broadens as the diffusion length increases with energy.

Zoom In Zoom Out Reset image size

In Table 1, e-folding distances for each energy according to the simulations are compared with the theoretical values, showing good agreement. In the theoretical calculations for a non-spherical shock we had to include the shock normal dependence on longitude. The diffusion length is then the ratio of the normal components of the diffusion coefficient and the velocity. The e-folding distances from the simulations were estimated with the help of Figure 11 by calculating the distance by which the intensity drops by a factor of e from its value at the shock. More importantly, past the TS there is an increase in intensities by a factor of two at 5 MeV n⁻¹ and a factor of three at 10 MeV n⁻¹, with a subsequent decrease on approach to the heliopause, where the ions are lost to the LISM. Like the shock precursor, this decrease is a consequence of the diffusive behavior of the particles.

Table 1. Diffusion Length Comparison for κ_⊥ = 0.01κ_∥

Energy	Theoretical	Simulated
(MeV n−1)	(AU)	(AU)
0.5	1.0	0.9
1.0	1.6	1.5
2.0	2.6	2.7
5.0	4.7	4.7
10.0	7.5	7.4

Download table as:  ASCIITypeset image

The right panel of Figure 11 is similar to the left panel, but here the results were obtained with κ_⊥ = 0.05κ_∥. In this case the intensities appear nearly constant in the heliosheath almost all the way to the heliopause. This is a consequence of enhanced radial transport reducing the contrast between the field lines connected to forward and backward portions of the shock. Evidently, diffusive transport must be strongly field aligned to produce a measurable effect on the intensities in the heliosheath. The implication is that the heliosheath plasma must be weakly turbulent, which also follows from Voyager observations (Burlaga et al. 2011).

Finally, we performed a preliminary comparison between the radial profiles observed by V1 in the heliosheath and our simulation results. Figure 12 (left panel) shows the radial profile of 4–6 MeV n⁻¹ ACR helium ions measured by V1. The right panel shows the simulated intensity of 5 MeV n⁻¹ ions for comparison. The rapid initial increase in the data between the TS and about 100 AU was likely a temporal effect due to the decreasing level of solar modulation (Cummings et al. 2008). After 100 AU, the intensity increased by a factor ∼3 before starting to go down again. The final drop at a heliocentric distance of 122 AU is known as the "heliocliff" (Krimigis et al. 2013; Stone et al. 2013; Webber & McDonald 2013) that was interpreted by some as the heliopause crossing (Gurnett et al. 2013; Swisdak et al. 2013). Such an unexpectedly early encounter of the heliopause is contrary to most MHD model predictions, where its distance is typically in the range of 140–150 AU (Pogorelov et al. 2004; Opher et al. 2006; Borovikov et al. 2011). Similarly, in our plasma background model the heliopause is located significantly farther than 122 AU. For this reason, the two figures are not quantitatively comparable. However, the total change in intensity from the end of the strong modulation period to the intensity maximum is similar in the data and the simulations.

Zoom In Zoom Out Reset image size

Despite significant development in the past eight years, the Voyagers' ACR observations at the TS and in the heliosheath are still not fully understood. Spectra at the TS appeared modulated at low energies in both V1 and V2 data. Subsequently, these spectra began to unfold in the heliosheath. In this work we theoretically modeled the spectral features of anomalous helium in the heliosheath using a three-dimensional, semi-analytic model for plasma and magnetic field background and a series of stochastic simulations of ACR phase-space trajectories. In agreement with the published results (Jokipii 1986; Steenkamp & Moraal 1995; Florinski et al. 2004), our results show that during the A < 0 minimum conditions, characteristic of the solar minimum of cycles 23/24, ACRs dominate at lower latitudes. Contrary to the model with a spherical shock, our non-spherical shock results show that more efficient acceleration of ACRs occurs on the heliotail facing part of the TS. The effect appears to be a consequence of the converging motion of magnetic field line connection points to the shock in the tail region, causing the particles to be trapped and accelerated at those locations (Guo et al. 2010).

Further, our results suggest that the Voyagers may have observed an enhancement in intensities in the heliosheath due to a drift effect. As the amount of solar modulation decreased with declining solar activity, more particles would be drifting toward the ecliptic plane, closer to the Voyager trajectories. Model-derived radial intensity profiles in the direction of V1 travel at a fixed time and show a modest enhancement in mid-energy ACRs for small values of the perpendicular diffusion coefficient (1% of κ_∥). This degree of enhancement is consistent with the count-rate increase measured by V1 between ∼102 AU and ∼115 AU. A fully time-dependent calculation including the effect of the declining phase of the solar cycle (not performed here) is expected to also reproduce the enhancement between 94 AU and 102 AU. One could place an upper limit on the total increase in intensity by subtracting the no-drift spectrum at the shock from the drift-enabled result computed for the middle of the heliosheath (left and right panels in Figure 10, respectively). This gives a factor of three increase at 5 MeV n⁻¹, which is consistent with the increase measured by V1 between ∼102 AU and ∼115 AU. Spectra computed in the V1 direction of travel some 30 AU past the TS are essentially power laws without any breaks or other features below the high-energy rollover, again qualitatively consistent with the recent observations (Cummings 2011).

It is instructive to perform some statistical analysis of particle trajectories using our model. Figure 13 shows the relative probability to reach 5 MeV n⁻¹ as a function of injection point location in the ecliptic plane. We considered two detection locations, one on the nose and the other on the tail of the TS. As is evident from this figure, the probability to reach ACR energies is higher for particles injected into the DSA closer to the location of the observer. The above result appears to contradict the simulations of Kóta (2012), where the nose direction was found to be a more favorable location for injection into the DSA process.

Zoom In Zoom Out Reset image size

Next, we calculated the acceleration times to reach 5 MeV n⁻¹ as a function of heliographic latitude, along the nose, flank, and tail directions. We performed simulations with 5 million stochastic trajectories and binned the results into 20 equally spaced latitude intervals ranging from 0° to 90° along different directions. The average acceleration times in each latitude bin were then computed and compared with the standard theoretical prediction for a planar shock (Axford 1981),

$$\begin{equation} \tau _a = \frac{3}{u_{n1} - u_{n2}} \int _{p_\mathrm{min}}^{p} \left(\frac{\kappa _{n1}}{u_{n1}} + \frac{\kappa _{n2}}{u_{n2}}\right) \frac{dp}{p}, \end{equation} \tag{ 31 }$$

where p_min and p are the initial and the final momenta, respectively, and u_n, κ_n are the velocity and the diffusion coefficient in the direction normal to the shock, respectively. Here the subscripts 1 and 2 refer to the upstream and the downstream regions, respectively. Figure 14 shows the acceleration time dependence on heliolatitude. The theoretical curves and the simulated results suggest that the acceleration times at the flanks are longer compared to the nose and the tail. In the latter two directions the normal component of the field B_n is equal to the radial field component B_r. At the flanks, B_r is no longer perpendicular to the shock, and B_n is larger because of the contribution from the (large) azimuthal magnetic field component B_ϕ. Consequently, κ_n is also larger, which results in longer acceleration times. Note that in the simulations, particles are accelerated everywhere along the TS, sampling regions of different κ_n. Equation (31) also does not take drift and adiabatic cooling in the expanding solar wind into account. For these reasons, we do not expect the simulation results to closely match the theoretically predicted acceleration times, although their dependence on heliolatitude is qualitatively similar.

Zoom In Zoom Out Reset image size

The above probability and the acceleration time calculations support our explanation of the ACR acceleration mechanism as follows. McComas & Schwadron (2006) proposed that in the nose region, field line shock connection time is insufficient to allow PUIs to become accelerated to ACR energies, whereas in the flank and tail directions there is ample time to convert PUIs to ACRs. However, as we showed above, the acceleration time at the nose is much smaller and the ACRs observed at the tail were mostly injected near the tail region. This suggests that acceleration time is not the real cause for the difference in ACR intensities in the nose–tail directions; rather, it is a result of the global magnetic field line topology discussed earlier.

Recently, McComas & Schwadron (2012) suggested that V1 entered a new region that is magnetically disconnected from the most distant parts of the TS, based on the dropouts in the ACR intensities observed prior to the heliocliff. In our model, the width of the region where field lines are connected directly to the flanks and the tail of the shock is much smaller (∼5 AU) compared to ∼20–30 AU suggested in their paper. For smaller diffusion coefficients, our results show continued increase in intensity well beyond the connection region. No rapid decreases in the heliosheath are present in our results, and the eventual decrease on the approach to the heliopause is gradual. The rapid drop of intensities observed by V1 at the "magnetic highway" boundary suggests that the perpendicular diffusion is small in the heliosheath. The V1 heliocliff encounter did not have all the signatures of the heliopause crossing (Burlaga et al. 2013; Krimigis et al. 2013; Stone et al. 2013; Webber & McDonald 2013), but the recent V1 plasma instrument observations suggest that V1 is now in the interstellar space (Gurnett et al. 2013).

Our current model has a number of limitations. It is time independent and has a flat current sheet geometry, which is a highly idealized scenario. In reality, the tilt angle is not zero even at the solar minimum; it also varies with the solar cycle. Therefore, our current results only have a qualitative similarity to V1 observation. By using a magnetohydrodynamic plasma background with realistic, time-dependent boundary conditions, one may be able to match Voyager observations more closely. Despite these limitations, our non-spherical model provides important clues to the interpretation of the unexpected Voyager ACR observations in the heliosheath. Our results indicate that very small perpendicular diffusive mean free paths are required to obtain a measurable intensity enhancement with distance from the TS in the heliosheath. This agrees with the theoretical model of Zank et al. (2004), which suggests that κ_⊥/κ_∥ ∼0.001–0.01 in the distant heliosphere. Small perpendicular diffusion implies weaker turbulence, which appears to be consistent with the observations of magnetic field fluctuations in the heliosheath (Burlaga et al. 2011).

This work was supported in part by NASA grants NNX10AE46G, NNX11AO64G, NNX12AH44G, and NNX13- AF99G and by NSF grant AGS-0955700.