Modelling heliospheric current sheet drift in stochastic cosmic ray transport models

Abstract

Drifts are one of the major cosmic ray modulation mechanisms in the heliosphere. Three types of drifts occur in the background heliospheric magnetic field, namely curvature, gradient and current sheet drifts. The last component occurs because of the switch in magnetic field polarity across the heliospheric current sheet and is the main topic of study. We discuss and implement a new approach to model drifts in a numerical modulation model. The model employs stochastic differential equations to solve the relevant transport equation in five (three spatial, energy and time) dimensions. What is of interest is the fact that the model can handle current sheet tilt angles up to the theoretical maximum of α=90° and still remain numerically stable. We use the additional insights gained from the numerical model to investigate the effectiveness of drifts along the current sheet by examining the relationship between the current sheet path length and the cosmic ray propagation time. It is found that diffusion can disrupt the drift process very effectively, leading to diffusive short circuiting of the current sheet by the cosmic rays.

Appendices

Appendix A: Minimizing the line segment

In order to calculate the HCS drift velocity of a CR scattering center located at an arbitrary point (r,θ,ϕ), a corresponding point on the HCS (r′,θ′,ϕ′) must be found such that the distance between these two points

(23)

with

(24)

is minimized, i.e. (r′,θ′,ϕ′) must be the closest point to (r,θ,ϕ) on the HCS. The properties of the HCS at (r′,θ′,ϕ′) will determine the directionality of $\vec{v}_{\mathit{ns}}$, while L will be used to calculate its magnitude. Because θ′ is function of r′ and ϕ′, L=L(r′,ϕ′) (being unique for each point (r,θ,ϕ)) and the minimization problem above reduces to a spatially 2D problem. Analytical solutions of (23) are impossible and therefore a numerical minimization (optimization) scheme is adopted. Panel (a) of Fig. 13 illustrates the optimization problem of finding the minimum value of L: The left asterisk is an arbitrarily chosen point (r,θ,ϕ), while the right asterisk is the corresponding closest point to the HCS, i.e. (r′,θ′,ϕ′). For this specific (left) point, L(r′,ϕ′) is shown in Fig. 14 as the surface and contour graphs. It is clear that a local minimum of L(r′,ϕ′) exists.

Fig. 13

Panel (a) illustrates the optimization problem of incorporating a wavy HCS into the model. The left point is the position of a pseudo-particle (r,θ,ϕ), while the right point is the corresponding closest point to the HCS (r′,θ′,ϕ′). Panels (b) and (c) show the same scenario, but for 20 randomly selected points with α=10^∘ and α=45^∘ respectively

Fig. 14

Illustrating the minimization of L as implemented in this work. The surface and contour graphs shows L(r′,ϕ′) for the point (r,θ,ϕ) in Panel (a) of Fig. 13. The series of triangles overlaid on the contour graph shows the Nelder-Mead optimization scheme at work in calculating the minima of L, indicated by the red asterisk

In this work, it is decided to optimize L(r′,θ′) numerically by using the Nelder and Mead (1965) search algorithm (also referred to as the amoeba or downhill simplex method). In this approach, three initial point are chosen, forming the initial simplex (triangle). Through a process of elongation, reflection and contraction, the initial simplex is updated with a new simplex, but with the vertices thereof closer to the minima of the function being minimized. The process is then iterated until the final simplex is sufficiently close to the minimum of e.g. L. It is found that 10 iterations of the algorithm is sufficient to minimize L for a variety of heliospheric test cases. A visualization of the Nelder-Mead scheme is projected onto Fig. 14: The initial simplex is shown in red, all subsequent triangles by blue and the calculated minimum by the red asterisk. This numerically calculated minimum corresponds the to right asterisk shown in Panel (a) of Fig. 13. Note that the HCS projections shown in Fig. 13 are done at an azimuthal angle of ϕ and does not precisely coincide with the actual position of ϕ′ which would vary for each test point (r,θ,ϕ). In panels (b) and (c) of this figure, L is minimized for 20 randomly selected points, with both (r,θ,ϕ) and (r′,θ′,ϕ′) shown for each pair of points. Panel (b) is for α=10^∘, while panel (c) is for α=45^∘. Just from inspection it is clear that the minimization scheme can in fact minimize L under a variety of conditions. The accuracy of this algorithm is however dependent on the choice of the initial simplex, and we choose this to be close to (r,θ,ϕ). This makes the convergence of the algorithm faster, but gives less accurate solutions if L gets very large. However, CRs do not feel the effect of the HCS for this scenario (see also (19)).

The minimization scheme discussed here is valid for α∈(0,π/2]. The case when α=0 (a flat, azimuthally independent HCS) is excluded as L does not have a local minimum in terms of ϕ′. The scheme can however still be applied to this case, but with the restriction that

$$\phi' \equiv\phi: \alpha= 0$$

(25)

reducing the Nelder-Mead search algorithm to a 1D line search in terms of r and r′.

Appendix B: Benchmarks for the modulation model

This section briefly describes benchmarks of the presently described CR modulation model with results from two other spatially 3D models, namely the models discussed by Burger (2012) and Pei et al. (2012). What is also of interest, is that these models use different numerical schemes to solve the CR transport equation; Pei et al. (2012) used SDEs (similar to this model), while Burger (2012) uses a more traditional finite difference numerical method. Moreover, the drift field is also handled differently in these models, with our approach being close to the one adopted by Pei et al. (2012).

Figure 15 shows a comparison between the results of these three models. All simulations are for galactic protons energy spectra at Earth in the A<0 drift cycle using modulation parameters identical to those used by Kóta and Jokipii (1983). Two model scenarios are shown, namely the case of a flat HCS (α=0^∘) and when α=30^∘. It is clear that the results from all three models agree well, vindicating the modelling approach (both the numerical scheme and the implementation of the drift field) as used in this paper. There is however a small discrepancy between the results for the α=30^∘ scenario. This is due to different expression for θ′ being used in the different models, with the discrepancy disappearing as α→0^∘.

Fig. 15

A comparison between proton energy spectra as computed by various spatially 3D drift models, namely: This study (squares), Pei et al. (2012) (circles) and Burger (2012) (broken lines). Two modelled scenarios are shown; the case of a flat HCS (α=0^∘) and a azimuthally varying HCS of α=30^∘