Cosmic ray transport near the heliopause

In this paper we summarize our modelling efforts for cosmic rays near the heliopause, and discuss whether galactic cosmic ray modulation beyond the heliopause is possible and present an explanation for the anisotropic nature of the observed cosmic ray intensities in the very local interstellar medium. We show that (i) modulation beyond the heliopause is possible, but highly dependent on the assumed parameters (most notable, the perpendicular diffusion coefficient). Treating the heliopause as a tangential discontinuity, significantly damps this modulation effect and leads to modelled results that are similar to Voyager 1 observations. (ii) By choosing an appropriate functional form of the perpendicular diffusion coefficient on the pitch-angle level, we are able to account for the anisotropic behaviour observed for both galactic and anomalous cosmic rays in the local interstellar medium.


Introduction
Anticipating the Voyager 1 (V1) crossing of the heliopause (HP), [1] questioned the use of a Dirichlet boundary condition for galactic cosmic rays (GCRs) at the HP. Previous estimates, e.g. [2], suggested that the presence of the heliosphere cannot alter the GCR distribution directly beyond the HP, so that the GCR boundary condition, the so-called local interstellar spectrum (LIS), can be prescribed directed there. The results of [1], and a follow-up paper by [3], suggested that the heliosphere can influence the GCR intensity beyond the HP significantly and that, once V1 crossed the HP, it could continue to measure a small positive radial gradient in GCR intensity.
In August 2012, V1 did indeed cross the heliocliff region, where anomalous cosmic ray (ACR) intensities decreased and GCR increased substantially [4,5]. This crossing was later confirmed as the HP [6]. Since then, V1 has been sampling the very local interstellar medium (ISM) in-situ and is by now ∼ 10 AU beyond the HP.
In this paper we summarize our GCR modelling efforts beyond the HP, after which we critically evaluate the parameters adopted in these simulations in lieu of the V1 GCR observations. We also discuss a plausible explanation for the observed anisotropic nature of both the ACR and GCR distributions directly beyond the HP. Figure 1. The left panel shows computed GCR proton energy spectra along the V1 direction for two different HMF polarity cycles, while the right panel shows energy spectra, only for the A < 0 drift cycle, along both V1 and V2 directions.

Galactic cosmic ray modulation beyond the heliopause
We make use of the hybrid model of [7], where the Parker transport equation [9] is solved, by means of stochastic differential equations (SDEs), on top of an MHD modelled asymmetric heliospheric environment. This allows us to capture the complex geometry of the heliosphere which cannot, very easily, be described analytically. Fig. 1 shows the resulting GCR energy spectra, using the model of [7], at different radial positions. The left panel is along the V1 direction and for different drift cycles, i.e. different polarities of the heliospheric magnetic field (HMF), while the right panel shows energy spectra, for only the A < 0 polarity cycle, along both V1 and Voyager 2 (V2) trajectories. The spectra at the HP (here corresponding to the 120 AU scenario) are significantly lower than the LIS values at low energies; at 100 MeV, the HP intensities are ∼ 25% lower than the LIS values. Moreover, we find that drifts have little or no effects beyond 75 AU, while the results are also very similar along the V1 and V2 directions.

Generic results
In order to examine what the main drivers of modulation in the ISM are, we performed a parameter study, varying the parallel mean free path (λ || ) and the ratio of the perpendicular (λ ⊥ ) to the parallel mean free path through the parameter η, defined by λ ⊥ = ηλ || . The results are shown in the right panel of Fig. 2, which shows the resulting modulation fraction, i.e. the ratio j HP /j LIS , when η is kept constant and λ || is changed (the blue curve) and when λ || is kept constant and η changed (red curves). We find the modulation most sensitive to η, which is explained by the draping of the interstellar magnetic field along the HP [8], illustrated in the left panel of Fig. 2: In order for GCRs to enter the heliosphere, they must diffuse perpendicular to the mean field. For larger values of η (i.e. more effective perpendicular scattering), the amount of modulation increases, and hence, the GCR gradient beyond the HP increases. We thus emphasize the importance of λ ⊥ in determining the gradient of GCRs beyond the HP. It must be noted that these results are refuted by [10] and [11] who found that: GCR modulation is significantly less pronounced than presented here and that the modulation process in the ISM is predominantly determined by the magnitude of λ || . The fact that [7] and [11] used similar models, but came to completely different conclusions, due to different results, is puzzling and is presently investigated. On the other hand, the SDE modelling of [12] is in agreement with the results of [7].

The heliopause as a tangential discontinuity
For most of the results presented in the previous section, we assumed a generic value of η = 0.02. Although this value is a good estimate near Earth, it is not necessarily so in the outer heliosheath and perhaps especially directly at the HP. The HP represents a tangential discontinuity, e.g. [13], separating heliospheric and interstellar plasmas, and may, as illustrated below, damp λ ⊥ significantly leading to an effective GCR propagation barrier.
To investigate the presence of this barrier, it is useful to consider other tangential discontinuities present in the heliosphere, most notable, stream interfaces associated with corotating interaction regions. It is well known that these stream interfaces are CR propagation barriers [14], although the exact mechanism responsible for this is still relatively unclear. Most authors postulate that the random motion of diffusing HMF lines is damped at the tangential discontinuity [15], thereby also damping the field line random walk (FLRW) process of perpendicular diffusion [16]. A schematic illustration of this is presented in Fig. 3: The left panel shows illustrative HMF lines diffusing under normal solar wind conditions. GCRs with gyro-radii smaller than the associated turbulence correlation length will then follow these field lines (e.g. the blue curve showing the trajectory of such a particle's guiding center), leading to motion perpendicular to the mean HMF. Near a tangential discontinuity (illustrated by the red line in the right panel), the diffusion of the field lines is damped (inhibited) near the discontinuity and the FLRW diffusion process cannot take place effectively. If the same process occurs at the HP, it is very conceivable that η 0.02, which, as shown below, can have large consequences for GCR transport in the ISM: Referring back to Fig. 2, and extrapolating to η → 0, it should follow that j HP → j LIS and no modulation takes place in the ISM. Recently, [12] studied GCR modulation beyond the HP and found similar results to [7]. They did, however, also investigate the possibility of small η values with their results summarized in Fig. 4: As η becomes smaller, the gradient beyond the HP diminishes, while for very small values η ∼ 10 −10 , a step-wise increase in GCR intensity is noticeable, directly at the HP, after which LIS values are quickly reached. This step-like increase is indeed reminiscent of the V1 data [4], which also showed no measurable gradient beyond the HP. Similar conclusions were reached by [11].

Cosmic ray anisotropies in the local interstellar medium
In the ISM, the observed ACR and GCR intensities show significant anisotropies [17]. Sector diagrams of these anisotropic distributions are shown in the bottom panel of Fig. 5 for ACRs (referred to as heliospheric particles) and GCRs (referred to as galactic particles), with the red arrow indicating the magnetic field direction. Note that both distributions are still isotropic inside the heliosphere. Although it is not surprising that these distributions are anisotropic in the ISM (λ || is expected to be very large there [18], which implies weak scattering, or equivalently, a small pitch-angle diffusion coefficient D µµ ), it is interesting that the anisotropies are in the opposite sense for both distributions: ACRs show an enhancement at pitch-angles of 90 • (or, in terms of the pitch-angle cosine, µ = 0), while GCRs show a depletion near µ = 0. [19] and [20] investigated the cause of the ACR anisotropies and could reproduce the observed features. Their model was, however, only applied to ACRs. We however believe that a common process  Figure 4. Results from [12], showing the modelled GCR intensity near the HP for varying values of η. For very small values of η a step-like increase in GCR intensity is found at the HP with no measurable gradient beyond it.
should be responsible for both the ACR and GCR anisotropies and such a process, as discussed by [21], is put forward in the following section.

Model results
In [21] we proposed that these anisotropies are simply due to more effective perpendicular diffusion across the HP at certain pitch-angles. By assuming that the perpendicular diffusion coefficient, on the pitch-angle level, has the functional form we showed that, at least qualitatively, the anisotropic behaviour of both cosmic ray species could be explained. Results from this model are shown in the top panel of Fig. 5, in the form of modelled sector diagrams. The modelled behaviour is in agreement with the observations: The ACR distribution, in the ISM, peaks at µ = 0, while the GCR distribution reaches a minimum here. The reasoning is that µ = 0 particles have the highest mobility across the HP (as was the assumption when Eq. 1 was postulated): The µ = 0 ACRs can thus readily escape across the HP, and hence, there should be an enhancement of these particles in the ISM. Similarly for GCRs, the µ = 0 particles enter the heliosphere more readily and a deficiency of these particles should be observed in the ISM.
Criticism towards this explanation for the observed anisotropies is mostly related to the adhoc fashion in which the functional form of D ⊥ was chosen. Such a form can, however, be derived from diffusion theory as is discussed in the next section.
As will also be shown in Section 3.2, this form of D ⊥ results from drift effects, and as such, generally increases linearly with rigidity, D ⊥ ∼ P . The model of [21] was again used to  Figure 5. The top panel shows modelled sector diagrams, from [21], for ACRs and GCRs, inside and outside of the heliosphere. The bottom panels are the corresponding observations from [17].
model the transport of ACRs across the HP, but now for two different rigidities, and the results summarized in Fig. 6. The left panel shows the intensity for ACRs, of two different rigidities, across the HP, with the higher rigidity particles' intensity decreasing first in front of the HP, and remaining higher in the ISM, in accordance with the V1 observations of [20]. We believe that these results show additional, qualitative, observational support for our explanation of the observed anisotropies. Moreover, when the anisotropy level is calculated (the bottom right panel), it is clear that: The magnitude of the ACR anisotropy increases with radial distance beyond the HP, peaking at some, parameter dependent, distance beyond the HP (the same is true for the GCRs and in accordance with the observations of [17]) and that the anisotropy is generally higher for higher rigidity particles. It is unclear whether these rigidity predictions for ACRs can be tested against observations. Interestingly, the GCR anisotropy was observed to disappear (the distribution was isotropized) and reappear from time to time in the local ISM. It is believed that times of isotropy closely coincided with the passage of transients past the spacecraft. In our modelling framework, these observations can be explained in two ways: (i) The additional scattering (enhanced turbulence levels) present in these transient structures can isotropize the distribution function. However, additional acceleration [23] will also be needed to maintain the relatively constant levels of the isotropic GCR flux. (ii) An isotropic distribution can be maintained, at the same constant level, if the effectiveness of perpendicular diffusion is diminished; without perpendicular diffusion, no GCRs will be removed from the ISM (lost to the heliosphere). Perpendicular diffusion can become less efficient if either the magnetic topology or turbulence levels change considerably with the passage of the transient (similar to a reverse Forbush decrease occurring beyond the HP -the passage of the transient can actually increase the GCR flux relative to undisturbed periods), or if the magnetic magnitude increases significantly (see Eq. 13). It is unclear which (if any or perhaps even both) of the above-mentioned isotropization processes are most viable in the interstellar medium.

Parameter motivation
Below, we present a derivation of the applicable functional form of D ⊥ near the HP. The derivation follows, to some extent, the methodology presented in [22], which is valid for particles with gyro-radii much smaller than the scale of the turbulence considered; i.e. D ⊥ is due to non-resonant effects. As usual, the mean and turbulent HMF is decomposed as with B 0 constant and δ B = 0. The fluctuating 2D component, believed to be primarily responsible for perpendicular diffusion, is generalized to include a component along the mean field with weak turbulence assumed, δB x,y,z B 0 . In the supersonic solar wind, the fluctuations are usually observed to be transversal in nature, i.e. δB z (x, y) ≈ 0 [24]. However, observations near the HP (in the so-called stagnation region) indicate that the fluctuations are mostly longitudinal (compressional) in nature [25]. As the HP can also be considered a tangential discontinuity, it is likely that the transverse components will be strongly damped directly at the HP, thereby also damping the FLRW process. We therefore assume the 2D component to be purely longitudinal in nature, i.e. δB x , δB y δB z , and proceed to calculate D ⊥ due to non-resonant interactions, i.e. for the scenario when the particle gyro-radius is smaller than the associated turbulent correlation scales. The mean-square displacement of particles across the mean HMF, due to random walking magnetic field lines, is defined as If homogeneous, magnetostatic fluctuations are considered, we have, at late times, a Kubo formula given by Assuming that the correlation function decays exponentially as it follows that Perpendicular motion is assumed to be due to random drift motion of the guiding center across the mean field Furthermore, the decorrelation time is taken to be the time it takes a guiding center to drift across a perpendicular correlation length, l ⊥ , so that D ⊥ becomes The guiding center drift velocity is, in general, given by whereb is a unit vector pointing along B and electric fields were neglected. For purely longitudinal fluctuations, we only get a contribution from the gradient drift term, leading to D ⊥ follows directly as which follows, in principle, the assumed dependence assumed ad hoc in the previous section. Also note the linear dependence on momentum (rigidity) contained in r L .

Discussion and conclusions
Our modelling illustrates that GCR modulation beyond the HP is possible, but the amount of modulation is parameter dependent and we argue that this is mostly determined by the level of perpendicular diffusion directly at the HP. If the HP is considered as a tangential discontinuity, and the perpendicular diffusion process is effectively suppressed, it seems unlikely that any measurable amount of modulation will occur beyond the HP. The original model of [7] predicted a radial gradient of ∼ 0.25 − 0.45 %/AU, assuming a ∼ 25 − 45 % decrease occurring over ∼ 100 AU. Currently, the presence of such a small gradient cannot be ruled out observationally. It should also be kept in mind that V1 cannot be considered yet to be in the pristine ISM, and that V2 data could show a different situation. By choosing the appropriate functional form for the pitch-angle dependent perpendicular diffusion coefficient, we are able to reproduce a myriad of observed features in the ISM related to the anisotropy of both ACR and GCR species.