Magnetic Trapping of Galactic Cosmic Rays in the Outer Heliosheath and Their Preferential Entry into the Heliosphere

This paper examines the geometry of interstellar magnetic field lines close to the boundary of the heliosphere in the direction of the unperturbed local interstellar magnetic field, where the field lines are spread apart by the heliopause (HP). Such field parting establishes a region of weaker magnetic field of about 300 au in size in the northern hemisphere that acts as a giant magnetic trap affecting the propagation of galactic cosmic rays (GCRs). The choice of an analytic model of the magnetic field in the very local interstellar medium allows us to qualitatively study the resulting magnetic field draping pattern while avoiding unphysical dissipation across the HP-impeding numerical magnetohydrodynamic (MHD) models. We investigate GCR transport in the region exterior to the heliosphere, including the magnetic trap, subject to guiding center drifts, pitch angle scattering, and perpendicular diffusion. The transport coefficients were derived from Voyager 1 observations of magnetic turbulence in the VLISM. Our results predict a ring current of energetic ions drifting around the interior of the magnetic trap. It is also demonstrated that GCRs cross the HP for the first time preferentially through a crescent-shaped region between the magnetic trap and the upwind direction. The paper includes results of MHD modeling of the heliosphere that provide the coordinates of the center of the magnetic trap in ecliptic coordinates. In addition to the heliosphere, we examine several extreme field draping configurations that could describe the astrospheres of other stars.


Introduction
The heliosphere is the best-known example of astrospheres that are regions carved out of interstellar space by a stellar wind.The Sun is known to be traveling relative to the surrounding interstellar cloud with a speed comparable to the fast magnetosonic speed in the cloud (Witte 2004;Frisch et al. 2009;Möbius et al. 2012;McComas et al. 2015;Zirnstein et al. 2016).The heliosphere is believed to have a comet-like shape, blunt in the direction of the Sun's velocity vector and extruded into a long (and possibly inhomogeneous) heliotail in the opposite direction (Parker 1961;Dessler 1967;Baranov et al. 1976;Wallis & Dryer 1976;Ripken & Fahr 1983;Suess & Nerney 1990;Baranov & Malama 1993;Zank et al. 1996Zank et al. , 2013;;Fahr et al. 2000;Izmodenov et al. 2005;Pogorelov et al. 2009).It was argued that the tail portion of the heliosphere could be bifurcated and turbulent (Yu 1974;Opher et al. 2015;Pogorelov et al. 2015;Opher et al. 2017;Kornbleuth et al. 2021) and some proposed a more round shape and a short tail based on observations of energetic neutral atoms (ENAs; Dialynas et al. 2017).In this paper, we adopt the conventional picture and focus our attention on the forwardfacing part of the heliospheric boundary, known as the heliopause (HP) and its flanks, but largely ignore the extended tail region.
The magnetic field far upstream of the heliosphere, B LISM (LISM = local interstellar medium, the region sufficiently far away that it is not affected by the presence of the heliosphere), is thought to be aligned approximately with the direction toward the center of the IBEX ribbon, according to the predominant theory of the ribbon based on secondary ENAs (McComas et al. 2009;Chalov et al. 2010;Heerikhuisen et al. 2010).For straight magnetic field lines in the distant LISM, this defines two directions: parallel and antiparallel to B LISM .While the model presented below is not affected by the reversal of the field direction, we define it in a way that is consistent with the Voyager observations, from south to north (Izmodenov & Alexashov 2020).
The moving heliosphere pushes the magnetic field lines apart thus creating a diamagnetic cavity surrounded by a region of compressed magnetic field known as draped field in the very local interstellar medium (VLISM) that is defined as the region affected by the heliosphere through the disturbances emitted by the HP and ENAs of heliospheric origin on outward trajectories.Another term is outer heliosheath (OHS), which is defined as the region between the bow shock, or bow wave, and the HP.In what follows we will use the two terms interchangeably.Draping tends to result in compression leading to a stronger field compared to its unperturbed value, but the opposite is true in the direction parallel to the LISM field.In that direction, the magnetic field could become considerably weaker than B LISM reaching a value of zero at the so-called null point in the limit of ideal MHD.This general trend could already be inferred from the early conceptual models of the heliosphere-LISM interaction (Holzer 1989).
The shape of the heliosphere is affected by the compression and bending of the field lines, and the entire structure is inherently three-dimensional.Three-dimensional MHD models are the preferred tool to study the global structure of the heliosphere and surrounding LISM.However, MHD models tend to experience difficulties near stationary tangential discontinuities, where grid dissipation leads to field diffusion that erases the difference between the LISM field and its heliospheric counterpart.The field direction appears to be most susceptible to this numerical effect, which results in a gradual rotation of the field well ahead of the HP, to eventually align with the heliospheric field (Opher & Drake 2013;Opher et al. 2017;Guo et al. 2019).This effect is easy to identify in simulation results because it disappears on switching off the heliospheric magnetic field.While it is possible that an interchange instability is responsible for the alignment of the fields at the tangential discontinuity itself (Krimigis et al. 2013;Florinski 2015), this process operates on much shorter, sub-astronomical unit scales, which are typically unresolved in a global MHD model.
The work of Chalov et al. (2010) was not subject to magnetic field diffusion and/or alignment by virtue of using a flowadapted mesh fitted to the global-scale discontinuities, including the HP, which is treated as a tangential discontinuity.Their results, showing a well-defined magnetic trap in the northern hemisphere, were an inspiration for this work.Another paper worth mentioning in this context is Strauss et al. (2013), who identified a VLISM region in the northern hemisphere where the interstellar magnetic field lines approach the surface of the HP at a nearly 90°angle, and argued that the location could be the preferred GCR access point to the interior of the heliosphere.The present paper represents a systematic analysis of the influence of the magnetic field draping in general, and of the magnetic trap structure in particular, on the transport of galactic cosmic rays (GCRs) throughout the VLISM, and on their incidence on the surface of the HP.
In this paper, we use an analytic model of magnetic field topology in the VLISM.The first analytic model of magnetic field draping was constructed by Parker (1961) as a superposition of uniform and dipole fields.That model was later improved by Whang (2010) by replacing a single dipole with a string of dipoles arranged along a semi-infinite line.In contrast to these constructed vacuum models, Isenberg et al. (2015) and Röken et al. (2015) developed models where the magnetic field was obtained by solving Faraday's induction equation for the plasma in a region of space exterior to the Rankine half-body.The latter is the simplest approximation to the shape for the HP, produced by a superposition of a point source and a uniform flow.The model of Isenberg et al. (2015) requires an evaluation of an integral along streamlines of the flow, which is timeconsuming for the type of application we consider below.Röken et al. (2015) were able to reduce their result to a closed algebraic expression that only requires an evaluation of a pair of elliptic integrals, which is an ideal situation for modeling charged particle transport without a computational grid, thus achieving arbitrary high resolution.The abovementioned models suffered from an infinite pileup of magnetic fields at the flow stagnation point, making them unsuitable for studies of particle propagation near the nose of the HP.The improved model, presented here, removes this obstacle by renormalizing the flow-normal component of the magnetic field based on a mathematical algorithm described in Section 2.
The paper also reports on the progress in developing a multiphysics framework for energetic particle transport in turbulent plasmas.GCR transport in the global heliospheric context is commonly investigated using nearly isotropic (also known as Parker) formalism (Ball et al. 2005;Florinski & Pogorelov 2009;Luo et al. 2013Luo et al. , 2016;;Strauss et al. 2013;Guo & Florinski 2014).While this approximation is mostly valid in the interior of the heliosphere, it can break down in the VLISM, where turbulence is weak (Burlaga et al. 2015(Burlaga et al. , 2018;;Fraternale et al. 2019) and pitch-angle scattering is infrequent (Florinski et al. 2023).Observations of the OHS by two different instruments on board Voyager 1 revealed that GCR protons experienced episodic depletions at pitch angles near 90°, lasting several months (Krimigis et al. 2013;Rankin et al. 2019).These depletions were apparently connected with weak shock-like compressions injected into the OHS as a result of large-scale solar wind structures, such as global merged interaction regions (GMIRs), impinging on the surface of the HP (Gurnett et al. 2015;Burlaga & Ness 2016).According to models, cosmic-ray particles could become trapped downstream of convex shocks, where they lose energy due to flow expansion (Kóta & Jokipii 2017;Zhang & Pogorelov 2020).The latter two papers used the so-called focused transport physics that accounted for magnetic focusing and reflection.Our new modeling framework generalizes this approach to permit simulations of the same physical phenomenon using a variety of transport methods, including, but not limited to, the two mentioned above.This opens up possibilities to numerically investigate a wide variety of charged particle transport problems with the same code base.
While the existence of the magnetic trap is predicted by models, observational confirmation is currently lacking.Voyager 1 is traveling in the northern direction and therefore in the general direction of the trap, while Voyager 2ʼs trajectory is directed away from it (see Section 3 for a diagram).If any signatures in GCR fluxes and/or anisotropies were present in Voyager 1 observations, they would have occurred shortly after the HP crossing.As we show below, magnetic field lines approaching the HP at nearly normal incidence spread out to cover much of the forward-facing portion of its surface.For this reason, particles observed near the HP would have likely passed through the magnetic trap interior before being detected.The resulting anisotropies are expected to be small, and outside the Voyager sensitivity range.Nevertheless, a more sensitive instrument mounted on a spacecraft on a trajectory in the hydrogen deflection plane (HDP, see below) could in principle observe a ring current produced by energetic particles drifting about the perimeter of the magnetic trap (Section 4), as well as directly measure the decrease in the magnetic field strength.In this paper, we show that the dominant fraction of cosmic rays entering the heliosphere is coming from one particular direction, and enter the heliosphere preferentially in the region between the nose of the HP and the null point (Section 6).These features make the direction toward the magnetic trap a priority for in situ exploration and the planners of the next generation of interstellar missions should give it serious consideration.

An Improved Analytic Magnetic Field Draping Model
Throughout this paper, we use a coordinate system with the Sun at the origin and the z-axis directed antiparallel to the LISM flow vector.All calculations use cylindrical (s, j, z) coordinates and are performed in a plane j = const.Röken et al. (2015) obtained an exact analytic solution to the ideal MHD Faraday's law for a flow around a Rankine half-body, which is a result of an interaction between two incompressible and irrotational flows, a point source representing the stellar wind (beyond the termination shock, since the incompressibility condition is invalid for the supersonic wind), and the subsonic interstellar flow that has passed through the bow shock or wave.This is the simplest possible description of the heliosphere-LISM interaction originally used by Parker (1961) to construct the first comet-shaped model of the global heliosphere.The flow potential is expandable into Legendre polynomials, of which only the first two terms are nonzero; the resulting flow field is where z 0 is the distance from the origin (i.e., the Sun) to the stagnation point on the HP, = + r s z 2 2 , and is the flow velocity at z = ∞ .Solution (1) is axially symmetric with respect to the LISM flow direction.The surface of the HP is a Rankine half-body, which is a quartic surface of revolution given by No forces (pressure gradients, magnetic stresses, etc.) are assumed to be acting on the flow.A more elaborate model now exists featuring a subsonic compressible flow velocity field (Kleimann et al. 2017); however, for the present qualitative analysis, the simplest incompressible description of the LISM flow is sufficient.Following Isenberg et al. (2015), one can decompose the magnetic field in the unperturbed flow B LISM into a flowparallel (longitudinal) component B l LISM and a flow-normal (transverse) component B t LISM and follow their advection around the obstacle separately.The initially parallel component will always remain parallel to the flow and therefore can be written as

LISM
This magnetic field reaches zero at the stagnation point.The initially perpendicular magnetic field component develops a longitudinal component of its own after having been advected around the HP.This component is ordered by isochrones, which are surfaces of constant flow travel time from some reference surface.Far upstream in the unperturbed LISM, these surfaces are planes normal to the flow vector.On approach to the HP, the surfaces become deformed into a bell shape.The HP itself is the limit isochrone corresponding to infinite travel time.This geometry is illustrated in Figure 1.We label the isochrone surfaces by the asymptotic value ζ = z(s → ∞).For a given spatial position (s, j, z), pointed to by an arrow in Figure 1, the label can be calculated as where a is the impact parameter (distance from the streamline to the z-axis far upstream) and F and E are incomplete elliptic integrals of the first and second kind, respectively.For the HP itself ζ = −∞.The transverse field component is given by Röken et al. (2015) as where The parallel and perpendicular magnetic field draping patterns are demonstrated in Figure 2. The left panel essentially shows the streamlines of the flow parted by the HP, while the right panel shows five isochrone surfaces that start as planes far upstream and become progressively more deformed with a decreasing value of ζ.
In the vicinity of the HP the ratio s/a → ∞ and the s and j field components in Equation (7) grow without bounds.This is a direct consequence of using the passive field transport limit.Fortunately, this infinite field pileup can be eliminated by the following mathematical procedure.Observe that far upstream of the heliosphere we are free to specify a perpendicular magnetic field with an arbitrary z dependence because all such solutions have zero divergence and satisfy Faraday's law.This also means that one can apply an arbitrary scaling factor to the transverse field based on the isochrone label ζ.In particular, one can choose the scaling factor to be a taper function constructed in such a way as to make the field finite at the stagnation point.In order to do this we require the expression for the transverse field (Equation (7)) in the limit of s → 0. This is readily shown to be We also require the z-axis intersection point of an isochrone with a label ζ (point z a in Figure 1).The intersection point is obtained by taking the limit s → 0 in Equation (4), which results in This nonlinear equation must be solved iteratively to obtain z a .
We also define a transition isochrone with the z-intersection value of z t ; the field is modified only for the points that have z a < z t .The unmodified transverse field is then multiplied by the ratio using Equation (9).The parameter α is given by where K is the desired ratio of the transverse component of the field at the nose of the HP B t (0, z 0 ) and the value at infinity B t LISM .Figure 3 demonstrates the use of such a taper using K = 8/3 and a purely transverse LISM magnetic field with a magnitude = B 3 t LISM μG.Note that the field may experience a kink at the transition isochrone.The effect is quite small for our choice z t = 1.3 z 0 .
Because the taper function is applied on a per-isochrone basis, the parts of the magnetic field thus modified will no longer tend to the undisturbed B LISM at large crosswind and tailwind distances.This fact, however, is of no relevance to the present investigation, which focuses on the upwind region of the VLISM.
The actual undisturbed magnetic field B LISM is neither parallel nor perpendicular to the direction of the LISM plasma flow u LISM relative to the Sun.The latter direction is thought to coincide with the inflow direction of interstellar neutral helium because helium atoms, unlike hydrogen, are only weakly affected by charge exchange when passing through the heliospheric interface.The helium inflow defines the −z direction in the present model.The flow of interstellar hydrogen becomes slightly deflected relative to the helium as a result of charge-exchanging with interstellar protons traveling around the HP.H and He velocity vectors define the hydrogen deflection plane (HDP).The direction toward the IBEX ribbon center belongs to the HDP and makes an angle of 132°with the helium inflow vector (Funsten et al. 2013). 15Because, according to the secondary ENA model, ribbon emission arrives from directions where B is normal to the line of sight, B LISM must lie close to the ribbon center.Using the global MHD model, Zirnstein et al. (2016) found that a simulation with a B LISM −u LISM angle (hereafter, the BV angle, θ BV ) of 140°and B LISM = 3 μG best reproduced the geometry of the ribbon.Similar values were used by Kim et al. (2017)  Figure 4 compares the field line draping patterns in the axisymmetric model for BV angles of 120°(left) and 140°( right).The figure shows the BV plane, which is the plane of symmetry in this model.The field lines pile up more densely in the southern hemisphere, but are more sparse in the north, where the region of weak field is clearly visible next to the surface of the HP.The field drops to zero at one point on the HP where the field line strikes it at 90°(the null point); this effect is similar to the plasma velocity becoming zero at the stagnation point.The blue-colored region in the north is the magnetic trap.It is the largest magnetic structure in space around the solar system, with a typical size of 100-200 au.In reality, the HP is compressed in the direction normal to the field, and the resulting configuration becomes asymmetric.Nonetheless, several recent global computer simulations that include magnetic forces on the plasma prominently feature the magnetic trap lying in the northern hemisphere (Chalov et al. 2010;Izmodenov & Alexashov 2020).The magnetic trap is expected to contain a population of GCRs undergoing bounce and azimuthal drift motion conserving the respective adiabatic invariants.In Section 4, we examine these periodic motions and identify their respective timescales.
Of the two configurations shown in Figure 4, we prefer the one on the right with θ BV = 140°because it yields fairly sensible values of the magnetic field magnitude just beyond the HP in the Voyager 1 and 2 directions, 4 and 6 μG, respectively, which is not too far off from the measured values (Burlaga et al. 2018(Burlaga et al. , 2019)).The BV angle of 120°, on the other hand, yields a model B value of 6 μG at Voyager 1, which is considerably larger than measured.Given the limitations of the axisymmetric model, however, one should not make inferences about the actual direction of B in the LISM from this argument.For the 140°BV angle, the distance between the Sun and the null point is 180 au.

Numerical Validation of the Analytic Field Draping Model
The real heliosphere is not axially symmetric, mainly because of magnetic stresses on the HP (but also due to north-south asymmetries in the solar wind and the tilt between the solar rotation axis and the VLISM flow direction).The purpose of this section is to verify whether the analytic model is a valid qualitative representation of the actual heliosphere, according to our best understanding of the physics of the interaction between the solar wind and the LISM.To this end, we performed a three-dimensional computer simulation of the global heliosphere using the Moscow model, which combines an MHD solution for the plasma flow with a Monte Carlo simulation of the neutral atoms' intensities (Izmodenov & Alexashov 2015, 2020).In that simulation, we took B LISM = 2.95 μG, θ BV = 140°, and the BV plane (coincident with the HDP) was rotated by 45°relative to the meridional (also known as polar) plane.Additionally, the interstellar temperature was T LISM = 7500 K, which is slightly higher than the 6530 K value used in Izmodenov & Alexashov (2020), and the updated charge exchange cross section by Bzowski & Heerikhuisen (2020) was employed.In the simulation, the 22 yr (1995-2017) averaged solar cycle conditions were utilized at the inner boundary (1 au) as described in detail in Appendix A of Izmodenov & Alexashov (2020).
Figure 5 shows the magnetic field magnitude and magnetic lines of force draped around the surface of the HP based on a simulation using the Moscow MHD neutral model.Clearly, there is a region of magnetic field enhancement in the southern hemisphere, and the magnetic trap is in the north.Both these features are reproduced by the analytic magnetic field draping model.However, as one would expect, the numerical model predicts a strong asymmetry of the HP with respect to the zaxis, a feature that is not captured by the currently used analytic model (but could potentially be addressed through its asymmetric generalization; Kleimann et al. 2016).The interstellar magnetic field pressure pushes the HP and termination shock toward the Sun, and an asymmetric structure of the HP is produced by the difference in magnetic pressure (Izmodenov et al. 2005).The null point in the MHD model is ∼460 au from the Sun, which is 2.5 times farther than in the analytic model (180 au).Note that the field lines appear to intersect the HP in the northern hemisphere; this is a plotting artifact caused by projection onto a plane (the actual 3D field lines are not confined to the plane).
We also show the cardinal directions and planes (that determine the boundary conditions for the model) as they would appear to an observer located at the origin (the Sun) on the all-sky map in the ecliptic (J2000) coordinates shown in µ Figure 3. Magnetic field strength on the axis as a function of distance from the HP.The dashed line corresponds to the unmodified field (Equation ( 7)), while the solid line shows the field with the taper applied.The transition isochrone has z t = 1.3r 0 (143 au) and yields the largest amplification factor of 8/3 relative to the unperturbed field, whose magnitude is 3 μG.The ecliptic latitude of the null point is similar to that of B LISM , but their longitudes differ by 90°.This is a result of the bending of the magnetic field lines in the direction of the nose.It implies that, whereas Voyager 1 is moving in the general direction of B LISM (black star in Figure 6), it is quite a distance away from the center of the magnetic trap shown with the "×" sign.
It could be argued that the HP is dynamic and turbulent because of variations in the solar activity (e.g., Pogorelov et al. 2021), and that no well-defined null point exists.However, the region of a weaker magnetic field should be a resilient feature in the models, whether directly data driven or based on longerterm averaged conditions in the solar wind.The shape of the magnetic trap and the geometry of the VLISM magnetic field could change with time, possibly enabling particles to enter and escape the trap region.

Adiabatic Motion of Charged Particles in the Magnetic Trap
Cosmic rays arriving from the north can readily access the interior of the trap by traveling along the magnetic field lines, regardless of their initial pitch angle in the LISM.On entering the trap their pitch angle will decrease to conserve the magnetic moment.In the vicinity of the null point, the particle is traveling nearly parallel to the field.Passing through the trap, the particle will eventually reach a region of a stronger field and be reflected back and trace its inbound trajectory in reverse.The region will also contain trapped particles with large pitch angles.In equilibrium and under adiabatic conditions these two populations add up exactly to the isotropic distribution in the distant LISM.
Several effects can change the situation.Pitch-angle scattering would couple the passing through and trapped populations.Disturbances from the HP could disrupt the geometry of the trap allowing some particles to escape, and those entering from the distant LISM take their place.Another important phenomenon is the presence of the heliosphere.Cosmic rays stored in the trap and undertaking periodic motion would spend significant amounts of time near the HP and have a greater chance of crossing its surface into the inner heliosheath (IHS).The time to escape the trap can be estimated as τ = λ ∥ /v, where λ ∥ ∼ 3 × 10 3 -10 4 au at 1 GeV (Florinski et al. 2023) is the scattering mean free path and v ; c is the speed of the particle.This yields typical trapping times in the range between 2 weeks and about 2 months at that energy.It is likely that the trapped particles comprise the bulk of cosmic rays entering the heliosphere.Some of these cosmic rays would be deflected by the magnetic barrier and return to the LISM.In the process, those particles would have undergone scattering in the turbulent IHS and left with a different pitch angle.Other particles might enter the supersonic solar wind region where they could lose energy in an expanding flow (or accelerate at  the terminations shock) and return to the LISM with different energy (Florinski & Pogorelov 2009;Strauss et al. 2011;Luo et al. 2016).
The motion of a single cosmic-ray particle in the limit of weakly nonuniform and slowly varying fields (r L = L, Ω −1 = T), where r L is the Larmor radius, Ω is the cyclotron frequency, and L and T are the characteristic length and timescales of the magnetic field variations (∼100 au and ∼1 yr, respectively), can be described with the help of relativistic equations of motion for the guiding center.This approach is more computationally efficient than full orbit integration and is less error-prone because there is no need to resolve gyromotion.For a particle with mass m, charge q, and momentum components p ∥ and p ⊥ with respect to the magnetic field direction ˆ= b B B, these equations read for a time-independent magnetic field as (Grebogi & Littlejohn 1984;Tao et al. 2007) where the modified electromagnetic fields are defined as 2 1 2 is the relativistic factor, and . Equation (13) contains the electric, curvature, gradient, and polarization drift motions, while Equation (14) includes the effects of parallel electric field and focus.Nonadiabatic effects, such as pitch-angle scattering, will be introduced in Section 5 and the results discussed in Section 6.
Figure 7 shows a guiding center path of a 1 GeV proton injected into the trap nearly perpendicular to the magnetic field line.The injection point was 10 au away from the centerline of the trap, defined here as the field line striking the null point, hereafter referred to as the null line (blue in Figure 7).The trajectory completed one full drift period lasting 8 yr; it is not closed because of the electric drift (this was verified by tracing the same particle with the electric field set to zero, which produced a fully periodic trajectory).Trajectories started farther away from the centerline take progressively longer to complete a full cycle; at the 1 GeV energy, the shortest time to orbit the centerline is 4-5 yr, while farther away the period increases to several decades.Even the shortest period is significantly larger than the typical time between scatterings, so one would not expect the particles at this energy to be contained long enough to complete a single orbit.Generally, the drift velocity is proportional to γ (or the rigidity R), while the quasi-linear scattering mean free path is proportional to R 1/3 (Jokipii 1971).For example, a 100 GeV proton would complete the orbit shown in Figure 7 in under 2 months, which is longer than the scattering time.The average mirroring cycle for this trajectory was about 9 days, over which a relativistic proton would travel a distance of 1500 au.
A useful measure of the strength of the magnetic trap is the mirror ratio (B B max min ) or the loss cone half-angle ( -B B sin 1 max ).We computed the values for the loss cone half-angle for all positions on the surface of a Rankine half-body that is slightly larger than the HP in the model (by 1 au at the nose and somewhat more in the tail).They are presented in Figure 8 as a color map superimposed on the Rankine surface.
Here red signifies very little confinement, while blue denotes only the particles with very small pitch angles that can escape from the trap.As one would expect, the loss cone is the narrowest near the null point, where the field strength is effectively zero.At its widest part (near the HP) the magnetic trap in the analytic model is some 300 au across and is elongated in the meridional direction.
Figure 6.Cardinal directions relevant to the VLISM-heliospheric interface as seen from the Sun.Solid color lines show the solar equatorial plane (cyan), the meridional plane (gray), the HDP (magenta), and the IBEX ribbon small circle (green).Also shown are the directions of interstellar hydrogen and helium flows, the hypothetical direction of B LISM , the IBEX ribbon center, the directions toward Voyagers 1 and 2, and the position of the magnetic null predicted by the MHD model.Lines show the solar equatorial and polar (meridional) planes and the HDP.The line starting at the "×" sign (the magnetic null) connecting it to the "+" sign (the magnetic field in the unperturbed LISM) plots the coordinates of the null field line; the color along the line indicates heliocentric distance, and the color map inset shows the magnetic field strength around the null point with a clear minimum at the center.

Scattering and Perpendicular Transport
The basic approach of adding pitch-angle scattering and perpendicular diffusion to a deterministic transport equation is similar to that in Florinski et al. (2013) and Zhang & Pogorelov (2020).Essentially, we add two stochastic processes in pitch angle and displacement normal to the magnetic field.For this work, we have generalized the earlier model by replacing its isotropic scattering model with a pitch-angle-dependent scattering coefficient.In addition to Equation (13)-( 15), one solves a pair of Langevin equations, and where μ is the pitch-angle cosine, D μμ is the pitch-angle Fokker-Planck coefficient, D ⊥ is the perpendicular Fokker-Planck coefficient, W μ and W ⊥ are the Wiener (random walk) processes in pitch angle and the normal direction, and ∇ ⊥ is the gradient operator acting in the direction normal to B. The advective terms in Equations ( 18) and (19) appear because of the dependence of the Fokker-Planck coefficients on phase-space coordinates.In practice, at every step during the trajectory integration, the numerical code first constructs an orthogonal coordinate system ¢ ¢ ¢ x y z with ˆ ¢ e B z and the remaining axes chosen at random in the normal plane.In this field-aligned frame, two 1D Wiener processes act independently on the ¢ x and ¢ y axes.The deterministic term in Equation ( 19) is combined with Equation (13) into a single advective transport term.We emphasize that the resulting model is not equivalent to solving the focused transport equation.One major difference is that in the former the particle momentum is measured in the inertial frame, while in the latter it is measured in the plasma frame.
The diffusion model used in this work is based on Florinski et al. (2023), who proposed a model of magnetic fluctuations in the VLISM consisting of three components, including Alfvénic (index "A"), transverse 2D ("T") and longitudinal 2D ("L").The fluctuations measured by Voyager 1 were very weak compared with the mean field, but had a power-law spectral distribution indicating the high wavenumber end of the inertial range (Burlaga et al. 2015(Burlaga et al. , 2018)).Florinski et al. (2023) showed that for prevailing VLISM conditions pitch-angle scattering is reasonably well described by quasi-linear theory (Jokipii 1966), while perpendicular transport is in agreement with the weakly nonlinear theory, originally developed by Shalchi & Schlickeiser (2004).Assuming the turbulent spectrum is an uninterrupted power law between k = k 0 and infinity and ignoring wave propagation effects, the pitch-angle scattering coefficient is where γ is the power law of the inertial range of the turbulence, and ( ) d á ñ B A 2 is the turbulent variance of the Alfvénic component in the parallel wavenumber spectral range between k 0 and infinity.With the addition of the longitudinal component (not considered in previous work), the expression for the perpendicular Fokker-Planck coefficient is  where v is the particle's velocity and ( the variances of the transverse and longitudinal turbulence components, respectively, in the perpendicular wavenumber spectral range between k 0 and infinity (note that we are using the same k 0 for the slab and 2D components here).At gigaelectronvolt energies, the latter has only a small effect on perpendicular diffusion, so the second term in Equation (21) can be often omitted.
The wavenumber k 0 in principle corresponds to the outer scale of the turbulence.This scale could be larger than 1000 au based on observations (Burlaga et al. 2018), which is certainly larger than the scale of the system we are considering.In that case, particles would experience faster-than-diffusive transport over shorter time periods (Florinski et al. 2023).However, recent theoretical considerations based on ion-neutral collision damping and neutral viscous damping place an upper bound to the VLISM turbulent injection scale at ∼200 au (Xu & Li 2022).Therefore, a conventional diffusive model was used with a range of plausible values for k 0 corresponding to the size of the large turbulent eddies in the system.
The real HP is not a smooth and steady surface, but a dynamically evolving structure with some regions developing vortex or plume-like features.According to dynamic MHD models of the heliosphere, the typical size of these features is between 10 and 100 au (Florinski et al. 2005;Borovikov et al. 2008;Pogorelov et al. 2017Pogorelov et al. , 2021)).However, using this as an estimate of the turbulence outer scale parameter k 0 results in a perpendicular diffusion coefficient that is too large to explain the steep gradient of GCR intensity measured by the Voyagers at the HP (Stone et al. 2013;Guo & Florinski 2014;Luo et al. 2016;Stone et al. 2019) because ~-D k 0 1 according to Equation (21).For this reason, we also included the case with variable k 0 in which the largest turbulent length scale drops to 1 au very close to the HP.The background magnetic field B changes significantly throughout most of the simulation domain.Thus, the transport coefficients are not constant owing to their dependence on B. Since we do not have a model of turbulence transport, we chose to keep the ratio 〈δB 2 〉/B 2 constant everywhere for a given k 0 .
Perpendicular diffusion enables particles to enter the heliosphere.The analytic magnetic field model we use does not extend into the interior of the heliosphere where the flow pattern is considerably more complicated (Yu 1974;Suess & Nerney 1990;Ratkiewicz 1992).In MHD-based simulations, where the heliosphere is just another region in the flow (Ball et al. 2005;Florinski & Pogorelov 2009;Strauss et al. 2011;Luo et al. 2016), particles crossing the HP eventually exit back into the VLISM, but their energy at the exit point is different from their initial energy because of cooling or acceleration experienced in expanding or convergent plasma flows, respectively.We cannot replicate this effect in our model, so the HP must be necessarily treated as an absorbing boundary.We are therefore mostly interested in the distribution of the first crossing points on the HP, which should provide the locations of the regions of preferred access.

Simulation Results
We performed Monte Carlo simulations by tracing particle trajectories in the VLISM using the stochastic equations arising from the combination of Equations ( 13)-( 15), with Equations ( 18) and (19).The background was static and given by the description in Section 2 with z 0 = 117 au since this value gave the correct HP crossing distance for Voyager 1.We imposed two absorbing boundaries, one at the surface of the HP and another at a spherical surface centered on the Sun with radius R outer = 500 au.From the calculations in Florinski et al. (2023) for the second interval in Burlaga et al. (2018), we estimate 〈δB 2 〉 ≈ 3.45 × 10 −14 G 2 for k 0 ≈ 2π/7 au −1 with 7% slab turbulence, 68% 2D transverse, and 25% 2D longitudinal.Therefore, we used several combinations of diffusion parameters, all consistent with these numbers, which are summarized in Table 1.Although yielding the same slope for the spectral density of magnetic turbulence, and thus the same parallel diffusion coefficient, each of these cases represents different relative values of parallel and perpendicular diffusion.
In the variable k 0 case, the outer scale of turbulence was variable in the space between the HP and another Rankine halfbody with z 0 = 118 au.The outer scale of turbulence was linearly interpolated between l 0 = 1 au and l 0 = 10 au within this region based on the z 0 value of the Rankine half-body that would contain the particle's current position, and the turbulence scale parameter was then computed via k 0 = 2π/l 0 .Beyond the second Rankine half-body, the outer scale of turbulence was held constant at 10 au.
We performed both forward-in-time and backward-in-time simulations with several million trajectories per run.All trajectories were binned based on their initial and final positions.For the forward-in-time runs, the particles were initialized on a spherical sector near the outer spherical boundary (except for the special cases described in the Appendix).Trajectories were followed until either of the absorbing boundaries was encountered, and only the particles crossing the HP were recorded and binned.For the backwardin-time simulations, particles were injected near the surface of the HP, and the HP crossing trajectories were rejected (i.e., giving a weight of zero) because in backward-in-time simulations we are only interested in the particles that arrived from the VLISM.
The forward-in-time and backward-in-time simulations generate different statistics.The forward-in-time runs calculate the probability that an interstellar particle would cross into the heliosphere for the first time, as a function of their starting position within a uniform and isotropic initial distribution at the outer boundary, as well as the distribution of first crossing points on the HP.Note that the former is a probability density on the outer boundary, while the latter is defined on the HP surface.In absolute terms, these probabilities depend on the value of R outer , with larger R outer values yielding a smaller solid angle of incoming directions for the cosmic rays, but the pattern should remain the same regardless of the boundary location, namely, a bulk of the probability centered around the magnetic field line but skewed toward the B LISM direction (see panels (g)-(i) in Figure 9).The backward-in-time simulations generate the number density of interstellar particles at a given   and Energetic Charged particle TRansport on Unstructured Meshes; Florinski & Alonso Guzman 2023). 16SPECTRUM was designed prioritizing modularity and scalability, so that the same infrastructure of code can readily model distinct species of particles, each perhaps obeying different transport physics, in a variety of electromagnetic environments, using a range of numerical methods from which the user can choose to balance accuracy and efficiency.For example, the same underlying code that we use to run parallel simulations of cosmic rays undergoing scattering and diffusion through a guiding center approximation in the turbulent VLISM with an adaptive integration method and absorbing boundary conditions could, by merely changing a few lines of code in the main file, utilize a single processor to model electrons traveling in the Earth's magnetosphere with custom initial and boundary conditions using the full Newton-Lorentz equation and a less accurate integration method for computational speed.It is also worth highlighting that multiple distributions can be binned simultaneously throughout one simulation, maximizing the useful output of each run.
The deterministic parts of the trajectories were integrated using a fifth-order accurate, adaptive step method of Dormand & Prince (1980), while the Maruyama scheme (Kloeden & Platen 1992) was used for the stochastic portion, which has first-order accuracy in the weak sense.In all cases, particles were 1 GeV protons with an isotropic initial distribution in velocity space.The software was configured using mostly the stock SPECTRUM components consisting of the Trajec-toryGuidingDiffScatt integrator (guiding center physics with scattering and perpendicular diffusion), the BackgroundVLISMBochum magnetic field background, two boundaries, BoundarySphereAbsorb at R outer and BoundaryRankineAbsorb on the HP, and the Initi-alMomentumShell (monoenergetic, isotropic source) injected distribution.The collected statistics and the initial space conditions were custom coded to obtain a density per unit area of the HP and to limit the injection area for the purpose of efficiency, as described below.
It should be noted that from a computational point of view, the percentage of useful trajectories was small (between ∼0.5% and ∼6%), but this can be misleading because the code is actually very lean and efficient.For this particular set of simulations, a majority of trajectories that were not binned only took a fraction of the time to terminate compared to the average time it took to complete a binned trajectory, so even though the percentage of binned trajectories could be 1%, a much larger fraction of the total computation time was spent integrating them.For example, in the backward-in-time simulation cases, most trajectories exited the simulation domain through the HP, which occurs extremely fast since the particles' starting positions are so close to it, and thus, these trajectories only account for a small fraction of the total computation time.Although these percentages should not be used to judge the performance of our algorithms, they still contain some valuable information.In the backward-in-time case, they can be interpreted as an estimate for what fraction of the 1 GeV population in the vicinity of the HP is from direct interstellar origin (specifically from the nose side), meaning they have not yet crossed the HP once.In the forward-in-time simulations, they serve to obtain a rough estimate of the fraction of 1 GeV cosmic rays coming from a direction near the interstellar magnetic field that penetrates the HP.

Forward-in-time Simulations
Forward-in-time simulations were performed with the initial positions of the particles distributed uniformly over a window in the shape of a spherical sector with a radius of 499.99 au and velocities distributed isotropically.That way only the trajectories with a finite, non-negligible probability of a close approach to the HP were considered.A non-recording time boundary of 1 yr was enforced, meaning that particles that have not crossed a spatial boundary within 1 yr of simulated time were terminated without binning.This was done to prevent any stray particles from consuming resources indefinitely, but very few particles ever reached the 1 yr threshold.This indicates that GCRs, even those trapped in the lower B region, do not remain in the system for long before exiting back into the LISM or crossing the HP at least once.Because the background GCR distribution is assumed to be uniform, in the forward-in-time case the resulting distribution is representative of the density of GCR first crossing points on the HP.Panels (a)-(c) in Figure 9 show the spatial distribution of heliospheric crossing locations for cosmic rays, along with their incoming directions in panels (g)-(i), for the three values of k 0 .Each binned value on the HP surface was normalized by the surface area of that bin in order to compensate for the nonuniform bin size along z.
The extent of the spherical sector was chosen based on the results from the backward-in-time simulation.In each case, 5 million particles were injected, out of which ∼0.95%, ∼2.0%, and ∼6.3% reached the HP for the turbulent length scales of 1-10, 10, and 100 au, respectively.These percentages multiplied by the normalized values in panels (g)-(i) in Figure 9 yield the probability that a pristine GCR intersects the HP at least once as a function of position.It seems reasonable that the lower k 0 values would lead to a larger number of particles ending in the HP since these correspond to higher perpendicular diffusion.The distributions of HP crossing points are shown in the top row; as discussed below, they are consistent with the number density maps from the backward-in-time runs, but not identical.
The most important result, clearly visible in the top panels in Figure 9, is that the probability distribution of first crossing points does not peak at the null point, but has a maximum closer to the nose.In fact, a particle inside the magnetic trap would have a low probability of coming directly from the LISM.These results are naturally explained by the magnetic field topology around the magnetic trap, which is discussed in more detail in Section 7. The distributions for smaller values of k 0 appear more noisy than those for large k 0 because there are fewer particles to bin when the perpendicular diffusion is large.The situation is reversed for the backward-in-time simulations discussed below.At this point, we would like to caution the reader again that the model only provides information about the first crossing points of GCRs approaching the solar system from galactic space.Almost all cosmic rays that enter the heliosphere eventually leave back into the LISM, and these particles have a finite probability of re-entering through the HP thus increasing the probability values shown.
The middle panels in Figure 9 compare the probabilities integrated along the azimuthal direction (top) and the z direction (bottom).The probability density is largest in the forward-facing part of the heliosphere (between z = 0 and z = z 0 ), but is still significant in the heliotail.One can also see that integral density is not the largest in the plane containing the magnetic trap (j = 0), even though it also contains the point of maximum density.This is a consequence of the bent shape of the high-density region that curves around the magnetic trap, so that more of it is located near j ∼ ± 40°.The large dip on the side opposite the trap corresponds to the magnetic wall that inhibits particle approach to the surface of the HP.
Notice how the GCRs of interstellar origin come from a very localized region (the heliospheric access window), as can be seen in the bottom panels in Figure 9. Since perpendicular diffusion is relatively small even for the smallest value of k 0 , the particles essentially follow the field lines when approaching the heliosphere from the LISM.The null point has a larger polar angle than the unperturbed B field (∼55°versus 40°), i.e., it lies farther down the tail.The field line passing through the magnetic null point envelops the entire surface of the HP, so it is not surprising that most of the particles entering the heliosphere originated on the field lines adjacent to the null line (see Strauss et al. 2013).The effect is obviously more prominent for large k 0 (i.e., small D ⊥ ), while smaller values of k 0 produce a wider distribution of HP intersecting trajectories on the outer boundary, meaning a larger range of latitudes and longitudes include particles that hit the HP.Another important trend is the shift of the window away from the null magnetic field line and toward the direction of B LISM .

Backward-in-time Simulations
Backward-in-time simulations were performed by initializing particles near the HP boundary, with a uniform pitch-angle distribution, and integrating their trajectories until an absorbing spatial boundary was crossed.The initial spatial distribution of particles was along a slightly larger Rankine half-body than the HP with the stagnation point distance z 0,IC > z 0 .Note that the distance between the HP and the surface of seed particles becomes larger away from the nose, but it is never larger than twice the distance between the nose points (z 0,IC −z 0 ).We did not insert any particles on this Rankine half-body farther than 499 au from the Sun in order to avoid an artificial result near the tail where it intersects with the spherical boundary.A nonrecording 1 yr time boundary was also implemented in this case, with very few trajectory realizations reaching it.
Figure 10 shows the distributions of initial (top row) and final (bottom row) positions for particles starting at a Rankine half-body with z 0,IC = 117.01 and integrated backward until they reached the outer spherical boundary on the side of the nose, indicating they are of certain interstellar origin.The results are arranged in the same way as those for the forwardin-time runs.Only those particles that exited the simulation through the nose hemisphere (z > 0) of the VLISM are considered in these results since those that exited near the flank region did so very close to the heliospheric boundary and many of them would have likely diffused into the HP farther down the tail.Each binned value on the HP surface was normalized by the number of particles initialized on that bin in order to compensate for the nonuniform spatial initialization along z.
The value of z 0,IC was chosen such that the initial position of the particle was within a Larmor radius away from the surface of the HP (r g ∼ 0.05-0.1 au for a 1 GeV proton).Other simulations using different values of z 0,IC yielded qualitatively the same results.In all backward-in-time cases, 10 milllion particles were simulated, out of which ∼3.0%, ∼0.83%, and ∼0.41% reached the VLISM boundary on the side of the nose for the turbulent length scales of 1-10, 10, and 100 au, respectively.These percentages multiplied by the normalized values in panels (a)-(c) in Figure 10 and the VLISM GCR density yield the number density of pristine GCRs around the HP.The relative percentages make sense since a lower k 0 corresponds to stronger perpendicular diffusion, which would cause a larger number of particles starting near the HP to cross this boundary before reaching the LISM.
The top row results from the backward-in-time runs agree nicely with those from the forward-in-time runs, with the density of interstellar particles that have not yet crossed the HP being concentrated between the magnetic trap and the nose.Notice how the distributions of exit points (bottom row) are also localized near the magnetic null line but shifted toward the B LISM direction, and so they do not match the probability distributions of reaching the HP shown in panels (g)-(i) in Figure 10.This contrast will be discussed in the next section.

Implications for Cosmic-Ray Entry into the Heliosphere
The forward-in-time simulations confirm the belief that GCRs incident on the HP from the nose side come from a relatively narrow region aligning with the interstellar magnetic field (Strauss et al. 2013), or rather the magnetic null field line, a sort of exit from the GCR highway into the heliosphere.Figures 9 and 10 reveal an area of low density of the first heliospheric crossing points and low number density, respectively, for pristine interstellar GCRs around the center of the magnetic trap.This effect is due to the extremely fast spreading of magnetic field lines around the null point.There seems to be a small spike of GCR crossings and density right at or near the null point (partially obstructed by the lime green null point markers), which could be exaggerated due to numerical effects, but may also be partially physical since this area contains the most perpendicular field with respect to the HP surface.
The region between the magnetic trap and the nose is indeed the most direct path into the heliosphere for cosmic rays incoming from the nose side, resulting in a higher crossing probability.Additionally, the HP crossing probability is relatively low on the side opposite the trap.This is because the region of high magnetic field acts as a mirroring point which GCRs have difficulty traversing.The results from the forward-in-time simulations display a slight but noticeable asymmetry in the distribution of heliospheric crossing points, especially for the higher k 0 cases (lower D ⊥ ).We believe this effect is physical and it is due to the drifting patterns of cosmic rays around the magnetic trap.
As we have pointed out, the forward-in-time and backwardin-time runs qualitatively share many features.However, the resulting distributions cannot be judged as identical, even allowing for statistical deviations.This is to be expected since these simulations compute different quantities.With respect to the HP maps (top panels in Figures 9 and 10), the backward-intime simulations calculate the number density of pristine GCRs around this surface, while the forward-in-time runs give the distribution of first entry points on the HP for these particles, both assuming a uniform VLISM population.In view of this, it would be fitting to say that these results are actually complementary to each other.In addition, the same model employed to find the crossing point density on the HP in the forward-in-time runs simultaneously solves the boundary value problem corresponding to the probability for a particle to exit through the HP as a function of its initial position in the VLISM, which is represented as a heatmap on the outer spherical boundary (bottom row in Figure 9).In contrast, the distributions of exit points in the backward-in-time simulations do not have a direct physical meaning.They are merely the exit points of the pseudo-trajectories used to find the number density near the HP.Nonetheless, we report them in this paper mainly to exemplify the differences in both approaches and highlight the correct interpretation for each.

Discussion and Conclusions
We have examined the interaction between the HP, shaped as a Rankine half-body, and the interstellar magnetic field draping around it through a simplified, axisymmetric model with a uniform plasma flow.As part of this effort, we showcased a mathematical technique to eliminate the unphysical singularity at the surface of the HP from the models of Isenberg et al. (2015) and Röken et al. (2015), and derived straightforward analytic expressions for the magnetic field in terms of flow-parallel and flow-perpendicular components.This particular geometry creates regions of low and high magnetic field magnitude on opposite sides of the HP, which act as a magnetic trap or bottle, and a mirroring obstacle for GCRs, respectively.The center of the magnetic trap is the null point, which is the location where a magnetic field line is normal to the surface of the HP.In this work we postulated a 3 μG magnetic field making an angle of 140°with the inflow direction in the unperturbed LISM (which corresponds to the often-used complementary BV angle of 40°); however, the results are expected to hold qualitatively for all stronger magnetic fields and smaller complimentary BV angles.
For a numerical verification of the analytic model, we performed a full MHD plus neutral atom simulation of the heliospheric interface using the same magnetic field properties at infinity as in the analytic model.The results had similar magnetic field magnitudes in the trap and at the magnetic mirror point.As expected, the full MHD heliosphere was asymmetric with a significantly larger distance to the HP in the north than in the south.In this particular simulation, the HP was some 460 au from the Sun in the direction of the null point, which is probably too far to be a target of in situ exploration by a single spacecraft using conventional propulsion, such as the Interstellar Probe (McNutt et al. 2019;Brandt et al. 2022Brandt et al. , 2023)).However, the distance would be reduced if the complimentary BV angle is less than 40°; the distance would be further reduced decades in the future if the solar cycle averaged solar wind momentum flux is on the decline, which appears to be the current trend (e.g., Verscharen et al. 2021).
Test particle simulations of 1 GeV GCR protons were performed using this analytic background model in order to explore multiple regimes of cosmic-ray transport in the VLISM.Initially, all nonadiabatic effects were switched off to investigate the shape of the unperturbed trajectories.It was found that particles placed in the trap undergo the bounce and drift motions corresponding to the conservation of the second and third adiabatic invariant, respectively.Drift results in a ring current where ions and electrons circle the magnetic trap in opposite directions.Because trapping does not depend on energy, it could be inferred that lower energy ions, such as pickup ions (PUIs) produced from heliospheric ENAs, could actually be born and confined inside the trap.Because the scattering rate in the nonrelativistic limit is proportional to R 2/3 , such particles would experience virtually no scattering on ambient fluctuations and will remain trapped for a very long time.For example, based on the model in Florinski et al. (2023), the time between scattering events for a 1 keV proton could be as large as 10 yr.If these PUIs are indeed the parent population of the IBEX ribbon ENAs, then the geometry of the trapping region is likely to have an imprint on the shape of the ribbon.The maximum ENA flux would be expected from the regions where the PUIs have a small field-aligned velocity, i.e., near the mirror points.These points form a circle around the null field line, but this circle is displaced in the tailward direction in the MHD model (see Figures 5 and 6) compared with the IBEX ribbon, which is centered approximately on B LISM .The displacement would be less if the BV angle is closer to 180°.
We note that trapping by the global magnetic field geometry has been already captured in the models that neglected turbulence (Chalov et al. 2010) or relied on weak ambient magnetic fluctuations for PUI confinement (Zirnstein et al. 2020).These authors have argued that no scattering on selfgenerated waves is required to confine the particles in the ribbon-emitting region (see also Giacalone & Jokipii 2015).If, however, waves are abundant as suggested by kinetic simulations of the ring instability (Florinski et al. 2010;Liu et al. 2012;Mousavi et al. 2022), then turbulent confinement could be the dominant mechanism for particle retention (Schwadron & McComas 2013;Isenberg 2014).In that scenario, PUIs with pitch angles close to 90°are scattered frequently; adiabatic trapping is ineffective, which means the subject of this work would be of lesser consequence for the ribbon studies.An insight could be gained by performing a comparison between the spatial extents of the ribbon and the trap regions using more detailed, data-driven MHD models of the global heliosphere.If the overlap regions have distinctive ENA signatures, it would provide support for the weak scattering scenario, and vice versa.
It was argued previously that the heliospheric bow shock could be of the slow MHD type, where the magnetic field strength decreases downstream (Florinski et al. 2004;Zieger et al. 2013).A slow bow shock would act as an upstream magnetic mirror and be a locus of mirroring points for particles with a wide range of pitch angles because of a steep magnetic field gradient.Because trapped PUIs spend more time near the mirror points, the chance of re-neutralization through charge exchange with interstellar neutrals is much higher at the mirror points, which would imply higher ENA fluxes from those regions.This could be a potential source of the IBEX ribbon.
For the next set of simulations, in conjunction with the adiabatic terms, the simulated cosmic rays were subject to stochastic weak pitch-angle scattering and perpendicular diffusion, which were implemented using Wiener processes.Both forward-in-time and backward-in-time simulations were performed for a variety of magnetic turbulence parameters for the purpose of determining the locations on the HP where GCRs can readily approach the boundary and subsequently enter the heliosphere from the VLISM.Our results indicate that GCRs incoming from the side of the nose arrive from a relatively narrow region in the VLISM generally aligned with the direction of B LISM .Cases in which perpendicular diffusion was higher naturally led to a wider spread of incoming directions.The distribution of the first HP crossing points of pristine GCRs has a clear dip near the magnetic trap and is highest in the region between the trap and nose.Moreover, the crossing probability of cosmic rays decreases along the flanks/tail with distance from the nose, and the side of the HP opposite the trap has a very low crossing density due to the presence of the magnetic mirroring region.We emphasize that the distributions presented here do not include cosmic rays that have previously crossed into the heliosphere and later escaped back into the VLISM.If they were, the total number density would be almost uniform in space.
Confinement by magnetic fields was previously called upon to explain the GCR depletions near the 90°pitch angle measured by Voyager 1.That initial idea due to E. C. Roelof (2013, unpublished) was further developed by other teams (Kóta & Jokipii 2017;Rankin et al. 2019;Hill et al. 2020;Zhang & Pogorelov 2020).It is natural to ask whether the magnetic trap discussed here could be a contributing factor.Note that it is difficult to produce an anisotropy by simply capturing some particles in a newly created magnetic trap if their prior distribution was isotropic and uniform in space (as expected of GCRs) because of Liouville's theorem.The situation is different if particles are injected into the trap at a specific location, as often happens in the Earth's magnetosphere.Particle distributions with a notch at 90°are commonly observed in the magnetosphere (Sibeck et al. 1987;Selesnick & Blake 2002), where they are believed to be a product of drift shell splitting.The latter is a consequence of the dependence of the drift trajectories on the pitch angle, so that an initially isotropic distribution becomes spatially separated and acquires a secondorder anisotropy away from the injection region.It would be difficult, however, to see how this mechanism could apply to the physical system considered in this paper.It would take several years for a particle to complete a loop around the center of the trap, during which time it would experience many scattering events.In addition, drift shell splitting could produce an enhancement at 90°, depending on where the particle injection occurred relative to the observer, but such GCR distributions were never observed.The possibility of a localized injection into the trap is also questionable in view of the high GCR mobility.
It is possible, in principle, that particles with μ ; 0 are preferentially transported across the HP owing to their large Larmor radii.Gradient drift preferentially affects particles with large pitch angles, so these particles can efficiently travel in directions normal to B, conceivably bringing them closer to the HP.This effect was used to explain the complementary pancake distributions of heliospheric particles leaking through the HP into the VLISM (Florinski et al. 2013(Florinski et al. , 2015)).However, a notch distribution requires an oppositely directed magnetic field gradient.This does not rule out the mechanism, however, because such a gradient could exist at a different location (the anisotropies do not need to be produced at the observation point), and it deserves further study.However, we think that the model of Kóta & Jokipii (2017) and its subsequent variations where cosmic rays with a large pitch angle lose energy in a post-shock expansion region is the most likely explanation of the measured anisotropies.
The region of the magnetic trap holds promise of new discoveries, and could potentially become a target of future multi-spacecraft missions, especially those using nuclear propulsion.As mentioned above, the rate of pitch-angle scattering increases weakly with rigidity, whereas the drift velocity of trapped particles increases linearly with it (see Equations ( 13)-( 16) and (20)).As a result, higher energy particles would exhibit stronger anisotropy (in the direction of the drift, counterclockwise for positive ions if viewed from the VLISM) than lower energy particles.If this effect could be measured, it would provide observational support for our model.reflected about azimuth = 0 in that case.Again only one of the two symmetrical windows is shown here.
We also investigated a simple case of a diamagnetic cylindrical surface embedded in a magnetic field perpendicular to its axis, with no plasma flow.It is easy to show that in this case the magnetic field is given by a potential solution, namely,  null point, this model has two parallel null lines at (x, y) = (± R, 0) on the opposite sides of the cylinder.The distribution of crossing points reaches a maximum in the direction of the unperturbed field (azimuthal angles of 0°and 180°; the second maximum is not seen because the particles were injected on one side only, see the discussion above).Panel (i) shows that nearly 100% of the particles injected near the xaxis were able to approach within a Larmor radius of the HP (only a single lobe about 0°is shown in panel (i), the second identical maximum is located on the opposite side at the 180°a zimuthal angle).The distribution has an almost symmetrical shape, so the drifts appear to play no role in this field configuration.
The cosmic-ray preferential HP access pattern for the oblique field discussed in the main body of this paper appears to have features from both extreme cases.If one were to start from the perpendicular case and increase θ BV , the northern lobe would develop a maximum approaching the nose, while the density of access points in the tail would decrease.The southern lobe would diminish and eventually disappear.The intermediate state is a crescent-shaped maximum visible in the top panels in Figure 9.For such angles, relatively fewer particles reach the heliosphere in the near heliotail (within a few hundred astronomical units from the Sun), but farther down the tail the situation should start to resemble the cylindrical case and the crossing point density should increase.With a further increase of the BV angle toward 180°(same as 0°) the maximum would broaden in the azimuthal direction and eventually encircle the symmetry axis as shown in panel (a) in Figure 11.

Figure 2 .
Figure 2. Field lines of B l (left) and B t (right) calculated using Equations (3) and (7), respectively.The former follows the streamlines, while the latter lies on isochrone surfaces.The surface of the HP is shown in orange.
in a datadriven model to simulate the shock waves observed in the VLISM by Voyager 1. Opher et al. (2020) used a similar θ BV , but a stronger field of up to 4.4 μG in their multicomponent MHD simulations.Izmodenov & Alexashov (2020) favored θ BV = 120°and a stronger B LISM = 3.75 μG because it provided the best match to the HP crossing distances by the two Voyager spacecraft.

Figure 6 .
Figure6.The HDP, shown in purple, is a great circle defined by the velocity vectors of the interstellar neutral helium and hydrogen, and B LISM is thought to lie in this plane.The figure traces the projection of the null field line onto the celestial sphere, from the null point to the outer boundary of the simulation at over 2000 au from the Sun.The ecliptic latitude of the null point is similar to that of B LISM , but their longitudes differ by 90°.This is a result of the bending of the magnetic field lines in the direction of the nose.It implies that, whereas Voyager 1 is moving in the general direction of B LISM (black star in Figure6), it is quite a distance away from the center of the magnetic trap shown with the "×" sign.It could be argued that the HP is dynamic and turbulent because of variations in the solar activity (e.g.,Pogorelov et al. 2021), and that no well-defined null point exists.However, the region of a weaker magnetic field should be a resilient feature

Figure 4 .
Figure 4. Magnetic field magnitude (in microgauss) and field lines in the BV plane for B LISM = 3 μG and the BV angles in the unperturbed flow of 120°(left) and 140°( right).The LISM flow is in the negative direction along the abscissa.The distance to the HP in the nose direction is 118 au.The field is enhanced in the southern hemisphere and reduced in the north.The blue-colored region in the north is the magnetic trap.The interior of the heliosphere, not covered by the present model, is colored uniformly in gray.

Figure 5 .
Figure 5. Magnetic field magnitude (in microgauss) and field lines (shown in black) in the BV plane.The simulation was performed for the BV angle of 140°a nd B LISM = 2.95 μG in the unperturbed flow using the kinetic MHD model of the heliosphere by Izmodenov & Alexashov (2020).The black dot shows the null point; its coordinates are at z = −150 au and x = 432 au.

Figure 7 .
Figure 7.One complete orbit of a 1 GeV cosmic-ray proton placed in a magnetic trap with an initial pitch angle of 80°(red line).The trajectory integration time was 3950 days.A single magnetic field line ending at the null point is shown in blue.The HP surface is colored in light gray.The direction of the drift is counterclockwise.

Figure 8 .
Figure8.Heatmap of the loss cone half-angle (in degrees) at a Rankine halfbody surface with a distance of 118 au to the nose, which is 1 au more than the surface used to delimit the HP.The null point is marked with a green dot, and several magnetic field lines at the periphery of the magnetic trap are shown for reference.

Figure 9 .
Figure 9. Final (panels (a)-(c)) and initial (panels (g)-(i)) position distributions for the forward-in-time integrated trajectories of 1 GeV protons that exited the simulation through the surface of the HP (i.e., particles that crossed into the heliosphere) for three values of k 0 .Trajectories were initiated in the spherical sector 35°< inclination <75°, −20°< azimuth <20°at 499.99 au.The density of the first crossing points is superimposed on the surface of the HP.The distribution of HP intersecting trajectories is plotted on the same surface used to initialize the particles (the heliospheric access window).Latitude and longitude are given in code coordinates, with the z-axis in the nose direction and B LISM in the xz-plane.The crossing point of the magnetic null field line through each plotted boundary is shown in lime green for reference.Panels (d)-(f), top and bottom, show the distributions of the crossing points as a function of z, integrated over the azimuthal angle j, and the distribution of the crossing points as a function of the azimuthal angle, integrated over all z.

Figure 10 .
Figure 10.Initial (panels (a)-(c)) and final (panels (g)-(i)) position distributions for backward-in-time integrated particles that exited the simulation through the outer spherical boundary on the nose side (i.e., particles of certain interstellar origin) for three values of k 0 .The distribution of the initial positions, initialized on a Rankine half-body with z 0 = 117.01au, is superimposed on the surface of the HP.The final position distribution is plotted on the surface of a sphere 500 au in radius centered at the Sun (the heliospheric access window).The crossing point of the magnetic null field line through each plotted boundary is shown in lime green for reference.Panels (d)-(f), top and bottom, show the pristine proton number density, in relative units, as a function of z, integrated over the azimuthal angle j, and the number density as a function of the azimuthal angle, integrated over all z.
the radius of the cylinder.We used R = 2z 0 = 234 au, which is the asymptotic value of the radius of the Rankine half-body for z→ − ∞.The value of B LISM was the same as in all other models.Selected field lines are shown in the xy plane in panel (c) in Figure 11.The particles were dropped 500 au away from the axis of symmetry with a radial spread of 60°(from −30 to +30).Because of azimuthal symmetry, only the cumulative distributions are shown (panels (f) and (i) in Figure11).Instead of a
The geometry of the heliosphere-VLISM interaction region from the analytic model demonstrating the isochrone-based magnetic field tapering method used in the code.The black curve is the surface of the HP.The point (labeled position) lies on the red isochrone ζ with a z-intersect value of z a .Field modification (see the text) is applied only downstream of the limit isochrone shown with a blue line.

Table 1
Diffusion Parameters Used in the Simulations