Theory of Cosmic Ray Transport in the Heliosphere

Modelling the transport of cosmic rays (CRs) in the heliosphere represents a global challenge in the field of heliophysics, in that such a study, if it were to be performed from first principles, requires the careful modelling of both large scale heliospheric plasma quantities (such as the global structure of the heliosphere, or the heliospheric magnetic field) and small scale plasma quantities (such as various turbulence-related quantities). Here, recent advances in our understanding of the transport of galactic cosmic rays are reviewed, with an emphasis on new developments pertaining to their transport coefficients, with a special emphasis on novel theoretical and numerical simulation results, as well as the CR transport studies that employ them. Furthermore, brief reviews are given of recent progress in CR focused transport modelling, as well as the modelling of non-diffusive CR transport.

handle both parallel and perpendicular diffusion (e.g. Owens and Jokipii 1971;Fisk 1976b), to be followed in the 1980's by 3D drift modulation models (e.g. Kota and Jokipii 1983). The latter illustrated the importance of particle drifts and the associated role of the wavy heliospheric current sheet (HCS), and gave eloquent explanations for the observed 22-year modulation cycle and the charge-sign dependent effect; see also the review by Rankin et al. (2022) on GCRs in this journal. In the 1990's, time-dependent modulation models for the TE followed (e.g. le Roux and Potgieter 1991), also suitably addressing acceleration effects caused by shock waves and stochastic motion of scattering centers (first and second order Fermi type acceleration, respectively) that were necessary to explain the acceleration of anomalous CRs (ACRs), as reviewed elsewhere in this journal by Giacalone et al. (2022). The stochastic interpretation of the CR transport equation, solving an equivalent set of stochastic differential equations (SDEs) (see, e.g., Jokipii and Owens 1975;Zhang 1999a,b) paved the way for more efficient computer utilization to solve multi-dimensional problems in CR transport. This approach to solving Parker's TE has contributed valuable insight into the understanding of the solar modulation of CRs, such as CR entry points, residence times, and acceleration sites. For reviews on the contribution made by time-dependent models, see Potgieter (1998Potgieter ( , 2013a) and on charge-sign dependence, see Potgieter (2014Potgieter ( , 2017. For a history and appreciation of the chronological development of finding solutions to Parker's TE, see the reviews by Quenby (1984), Moraal (2013), Potgieter (2013b) and Kóta (2013). These reviews provide a much more detailed discussion of aspects of CR acceleration and transport not covered here. For a comprehensive discussion on the SDE approach to solar modulation, see Strauss and Effenberger (2017).
The Parker TE is given here in terms of the omnidirectional GCR distribution function f 0 , a function of position r, momentum p, and time, by where this distribution function is related to the CR differential intensity per unit energy J T via f 0 (r, p, t) = p −2 j T (see, e.g., Moraal 2013, and references therein). The term Q represents any potential sources (or sinks) of cosmic rays, an example of which being the low energy electrons originating in the Jovian magnetosphere (see, e.g., Ferreira et al. 2003b;Nndanganeni and Potgieter 2018;Vogt et al. 2020Vogt et al. , 2022, and references therein). The terms in Eq. (1) describe various physical processes which influence the heliospheric transport of GCRs: the term v sw · ∇f describes the outward convection of CRs by the solar wind travelling at speed v sw ; CR adiabatic energy changes are described by (1/3) (∇ · v sw ) ∂f/∂ ln p; and ∇ ·(K · ∇f ) describes both cosmic ray diffusion and cosmic ray drift along the HCS and due to gradients and curvature in the 3D heliospheric magnetic field, which are incorporated into the 3D diffusion tensor K. Reviews of these processes, and our growing understanding of their influence on GCR modulation, can be found in Quenby (1984), McDonald (1998), Fisk et al. (1998), , Potgieter et al. (2001), Kóta (2013). It has become clear that diffusion and drifts play the most significant role in heliospheric GCR transport (e.g. Kota and Jokipii 1983;Potgieter and Moraal 1985;Potgieter 1996;Burger et al. 2000;Potgieter 2014;Engelbrecht and Burger 2015b;Kopp et al. 2017;Moloto et al. 2018;Shen et al. 2019). The diffusion tensor can be written in the heliospheric magnetic field (HMF) aligned coordinates as K (see, e.g., Kobylinski 2001;Burger et al. 2008;Effenberger et al. 2012), where allows one to model diffusion parallel and perpendicular to the HMF, as well as drift effects, via the diffusion (κ and κ ⊥ ) and drift (κ A ) coefficients, which can be related to corresponding mean free path (MFP, λ) lengthscales by κ = vλ/3, where v is the particle speed (see, e.g., Shalchi 2009). In the past, many studies modelled GCR diffusion coefficients in a phenomenological manner, citing a lack of theory and observations for these quantities as justification, and motivated by the fact that such an approach leads to computed GCR differential intensities in excellent agreement with observations while yielding useful insights as to global GCR modulation. Over the last few decades, significant theoretical advances have been made in our understanding of how solar wind turbulence influences the diffusion of charged particles (e.g. Shalchi 2006;Qin 2007;Shalchi 2009Shalchi , 2010Shalchi 2020Shalchi , 2021, as well as their drifts along the HCS and due to gradients in, and curvature of, the heliospheric magnetic field (e.g. Bieber and Matthaeus 1997;Burger and Visser 2010;). This coincided with significant advances in our theoretical understanding of solar wind turbulence, driven by an ever-increasing number of in situ observations (see, e.g. Marsch 1991;Frisch 1995;Goldstein 1995;Horbury et al. 2005;Matthaeus and Velli 2011;Bruno and Carbone 2016;Oughton and Engelbrecht 2021;Matthaeus 2021), as well as great refinements in turbulence transport modelling, with models progressing from the single-component models of e.g. Zhou and Matthaeus (1990), Zank et al. (1996), and Breech et al. (2008) to complex, two-component models that take into account the observed (e.g. Bieber et al. 1996;Oughton et al. 2015) anisotropic nature of solar wind turbulence (see, e.g., Oughton et al. 2006Oughton et al. , 2011Zank et al. 2012. As a result of these efforts, reasonable agreement was obtained with multiple spacecraft observations of various turbulence and turbulence-related quantities in different regions of the heliosphere (see, e.g., Smith et al. 2006;Engelbrecht and Burger 2013a;Adhikari et al. 2014;Usmanov et al. 2014;Adhikari et al. 2015Adhikari et al. , 2017aEngelbrecht and Strauss 2018), including Parker Solar Probe observations taken very close to the Sun (e.g. Adhikari et al. 2020;Chhiber et al. 2021c), making it possible to develop GCR modulation models of ever-increasing complexity. Such ab initio cosmic ray modulation modelling involves combining models for both the large scale (such as the heliospheric magnetic field) and small scale (such as the various statistical properties for solar wind turbulence, like the magnetic variance) plasma quantities in the heliosphere with transport coefficients derived from first principles in a manner that is observationally-driven and as theoretically and physically self-consistent as possible, in order to achieve agreement with spacecraft observations of cosmic ray intensities in different regions of the heliosphere, and to glean new insights into the fundamental physics underlying the transport of these particles. Consequently, it is one of the objectives of this review to briefly describe these advances in our understanding of ab initio approaches to modelling GCR transport coefficients, and to provide a short review of their application to the study of GCR transport and modulation (Sect. 6). Another objective is to review recent progress in our theoretical understanding of CR transport coefficients that have made such studies possible (Sects. 4 and 5). Furthermore, the results of various estimations of CR transport coefficients, from modelling and direct simulation, are also reviewed (Sect. 3).
Pitch-angle scattering will isotropize any (possibly anisotropic) CR distribution within a length-scale comparable to the parallel mean free path λ . One can therefore generally expect an isotropic cosmic ray distribution when examining particle behaviour at scales larger than λ . However, if a process is actively driving anisotropic behaviour, on a length-scale of L characteristic of whatever that process may be, the ratio L/λ determines when a distribution will be (nearly) isotropic. For instance, when L/λ 1 we tend to find strong deviation from isotropy. In such situations, the Parker TE is insufficient to describe CR transport, and a so-called focused transport approach, which describes particle transport on the pitch-angle level, is needed. Much recent work has been done to study the focused transport of cosmic rays, and as such, a brief review of these advances beyond the Parker formulation is presented in Sect. 7. Furthermore, an increasing amount of observational evidence has come to light, as reviewed elsewhere in this journal, that particle transport in certain circumstances becomes non-diffusive. This phenomenon presents unique theoretical and modelling challenges, which are reviewed and discussed in Sect. 8 of this review, which closes with a section devoted to a summary, and discussion of the future outlook of theoretical studies in GCR transport.

Turbulence and Particle Transport
An in-depth discussion of turbulence transport models (TTMs) and their development is beyond the scope of the present review (the reader is referred to a review by Fraternale et al. 2022 included in this journal), and as such, their treatment here will be brief, and focused on their application to the studies of CR transport. Furthermore, a motivated reader is referred to Breech et al. (2008), Oughton et al. (2011), Zank et al. (2012, Usmanov et al. (2018) for details on the derivation of specific TTMs. As a point of departure, TTMs are usually derived from compressive, single fluid MHD equations (see, e.g., Zhou and Matthaeus 1990;Breech et al. 2008), assuming a decomposition of the relevant variables into the sum of a large-scale, uniform component and a fluctuating component, represented using Elsässer variables. Scale-separated MHD equations are then written that include the effects of the large-scale flow compression and shear as well as various turbulence source terms, which may model, for example instabilities due to flow shear (e.g. Roberts et al. 1992;Zank et al. 1996;Ruffolo et al. 2020) or wave generation due to the formation of pickup ions (e.g. Lee and Ip 1987;Williams and Zank 1994;Isenberg et al. 2003). Various assumptions can be made regarding the nature of the fluctuations, i.e. whether they are incompressible (e.g. Oughton et al. 2011) or nearly so (e.g. . Various further assumptions pertaining to the closure of the resulting set of highly non-linear equations (see, e.g., discussions in Breech et al. 2008;Oughton et al. 2011;Zank et al. 2012;Oughton and Engelbrecht 2021;Adhikari et al. 2021a), such as an assumption that large-scale gradients in the solar wind are mainly radially directed, lead to a set of radial differential equations that, although not often analytically solvable (unless with simplifying assumptions, e.g. Zank et al. (1996)), can be solved using relatively straightforward numerical techniques, using as inputs either prescribed background models (i.e. the Parker (1958) heliospheric magnetic field) (see, e.g., Breech et al. 2008;Engelbrecht and Burger 2013a;Zhao et al. 2017;Adhikari et al. 2021b), or the outputs of global MHD models (e.g. Florinski and Pogorelov 2009;Usmanov et al. 2012;Wiengarten et al. 2015Wiengarten et al. , 2016Usmanov et al. 2011Usmanov et al. , 2014Usmanov et al. , 2016Usmanov et al. , 2018Chhiber et al. , 2021b, thereby providing a powerful tool for studying turbulence throughout the heliosphere (e.g. Smith et al. 2001;Ng et al. 2010;Pine et al. 2020), as well as information about the spatial behaviour of the various turbulence quantities that could be used to compute the diffusion coefficients based on theory, that would not necessarily be available from in situ spacecraft observations. Initially, turbulence transport models yielded results only pertaining to quasi-2D fluctuations (see, e.g., Oughton et al. 2006;Hunana and Zank 2010). More recent models, such as those proposed by Oughton et al. (2011), Zank et al. (2012, solve for both quasi-2D and wavelike turbulence quantities. This results in a more consistent modelling of the spatial variations Fig. 1 Comparison of magnetic variances and correlation scales computed using the Breech et al. (2008), Oughton et al. (2011 and Usmanov et al. (2016Usmanov et al. ( , 2018 turbulence transport models, employing the assumptions and boundary values used respectively by Pei et al. (2010), Engelbrecht and Burger (2013a), Adhikari et al. (2017a), and Usmanov et al. (2016, 2018. Observations shown are those reported by Zank et al. (1996), Smith et al. (2001) and Weygand et al. (2011) of anisotropic turbulence quantities and allows for an inclusion of the contribution due to pickup ions into the wavelike component (see Hunana and Zank 2010). Given the complexity of turbulence transport modelling, it is encouraging that the results yielded by several such models converge on spacecraft observations, as can be seen in Fig. 1, which compares the results obtained using the models by Breech et al. (2008), Oughton et al. (2011), and Usmanov et al. (2016, 2018. The figure shows magnetic variances (top panel) and correlation scales (bottom panel), that are two principal turbulence quantities used as inputs for theoretically-derived diffusion coefficients, as a function of heliocentric radial distance in the ecliptic plane, alongside spacecraft observations of these quantities. These models were run using the assumptions and boundary values employed by studies that specifically used their turbulence outputs to study diffusion coefficients, Pei et al. (2010), Engelbrecht and Burger (2013a), Adhikari et al. (2017b), and, respectively. These models include the influence of pickup ion generated waves, which can be seen clearly in the flattening of the radial dependences of wavelike, or slab, variances yielded by the two two-component models at radial distances beyond several astronomical units. In terms of the respective radial dependence, model outputs often disagree (especially regarding the correlation scales), owing to differences in assumptions of incompressibility of the fluctuations (e.g. , differences in parameters pertaining to the pickup ion source term, and differing boundary values (for a discussion, see Oughton and Engelbrecht 2021;Adhikari et al. 2021a). These differences would be reflected in diffusion coefficients calculated using each model's output (see Sect. 5.2), that would in turn have profound implications for CR transport (see Sect. 6).

Estimates of Particle Transport Coefficients
In parallel to the developments in our understanding of turbulence and its transport, much progress was achieved on understanding the CR transport coefficients. These studies were informed and constrained by a large number of numerical simulations using test particles in a box, as well the information on transport coefficients inferred from various heliospheric particle transport studies based on the criterion that the computed results should be in agreement with the spacecraft observations of particle intensities. In what follows, both of these approaches are briefly reviewed.

Estimates from Comparisons with Observations
The diffusion coefficients derived from scattering theories are also, to a degree, constrained by the values reported by numerous studies of particle transport modeling in comparison with observational data, particularly those studying the transport of solar energetic particles (SEPs). The best known such study is the Palmer (1982) consensus range for cosmic ray MFPs parallel and perpendicular to the interplanetary magnetic field (IMF). These estimates, however, are calculated in disparate ways. Often, they are reported as values required in SEP transport models to match the observed intensity and/or anisotropy profiles of the SEPs observed at or near Earth and beyond (see, e.g., Bieber et al. 1994;Dröge 2000;Dröge 2005;Agueda et al. 2008;Malandraki et al. 2009;Artmann et al. 2011;Tautz et al. 2016;Dröge et al. 2018;Aran et al. 2018), but values from cosmic ray observational (e.g. Bieber and Pomerantz 1983;. Other sources include transport/modulation studies of galactic and anomalous/heliospheric CRs (e.g. Dwyer et al. 1997;Burger et al. 2000;Cummings 1999;Cummings and Stone 2001), Jovian electrons (e.g. Ferreira et al. 2001a,b;Vogt et al. 2018;Nndanganeni and Potgieter 2018;Vogt et al. 2020Vogt et al. , 2022 and pickup ion anisotropy observations (Gloeckler et al. 1995). The intensity and anisotropy temporal profiles of SEPs can also be inferred from neutron monitor observations during ground level events (GLEs). These profiles can then be fit using SEP transport models of varying complexity, which in turn yields the estimates of these particles' pitch angle diffusion coefficients, from which the parallel mean free paths can be inferred (see, e.g., Bieber et al. 2002Bieber et al. , 2013. There are, however, considerable disparities in the values reported for particle MFPs (see, e.g., Bieber et al. 1994), and some caution is advised in interpreting the values arising from SEP transport models because they are strongly influenced by the various assumptions inherent in such studies. One example is the assumption of how the solar energetic particle source release rate is modeled, while another is the omission of perpendicular diffusion often made in early SEP transport studies. More recent studies, however, include this mechanism as a rule (see, e.g., Dalla et al. 2003;Zhang et al. 2009;Strauss et al. 2016;Laitinen et al. 2016;Laitinen and Dalla 2019;Chhiber et al. 2021b). Furthermore, the influence of solar cycle-related effects on these results is not fully understood. Several reports show large variations of transport coefficients over a solar cycle (e.g. Cummings and Stone 2001). This is not surprising given the strong association between SEP events and solar activity (e.g. Klein and Dalla 2017) and the observed solar cycle dependence of the HMF turbulence (e.g. Engelbrecht and Wolmarans 2020;Burger et al. 2022). Estimates based on comparisons of transport model results with Jovian electron observations can give vastly different results depending on the level of magnetic connection between the Earth and Jupiter, as well as the level of solar activity during the observation  Palmer (1982), Bieber et al. (1994), Gloeckler et al. (1995), Möbius et al. (1998), Dröge (2000), Kartavykh (2009), Dröge et al. (2014), and Vogt et al. (2020). Parallel MFP values calculated from analyses of GLEs are taken from Bieber et al. (2002), , and . Perpendicular mean free path values shown are those reported by Chenette et al. (1977), Palmer (1982) (see also Giacalone 1998), Burger et al. (2000), Dröge et al. (2014), and Vogt et al. (2020) period (Vogt et al. 2018, and references therein). It should be pointed out that many of the models have gone beyond the diffusion approximation and solved the equation of focused transport (Roelof 1969;Schlueter 1985;Ruffolo 1995), which allows for more physically consistent modeling of time profiles of SEP intensities and anisotropies (see also Sect. 7). Such studies have generally reported their results in terms of a mean free path that would correspond to a spatial diffusion coefficient.
As an illustration of the considerable variation in values reported for parallel and perpendicular MFPs in the literature, Fig. 2 shows a (by no means complete) sample of such values at 1 au as function of rigidity P , taken from various studies as indicated in the legend. Note that the original references for the data taken from Bieber et al. (1994) can be found in that paper, and that information pertaining to the individual events corresponding to the Dröge (2000) values for the parallel MFP can be found in that particular study. The parallel MFP values, shown in the top panel of Fig. 2, span a range of about two orders of magnitude, and vary considerably depending on the modelled event they correspond to (see, e.g. Giacalone 1998), with values often larger than the Palmer (1982) consensus range. Comparison of results corresponding to proton and electron studies led Bieber et al. (1994) to conclude that this consensus range is applicable to electron parallel MFPs at rigidities below 25 MV, and to protons above it, with electron MFPs displaying a rigidity-independent parallel MFP (see also Dröge 1994;Potgieter 1996;Evenson 1998;Dröge 2003). At higher rigidities, proton parallel MFPs appear to display (see, e.g., the data points from Dröge 2000) a ∼ P 1/3 dependence, expected from magnetostatic quasilinear theory (QLT) under the assumption of a turbulence power spectrum with a Kolmogorov inertial range (e.g. Bieber 2003). Although this particular theoretical result can reproduce the observed parallel MFP rigidity dependence, it yields values considerably smaller than observations if pure slab turbulence is assumed. This particular issue was resolved by Bieber et al. (1994), who found that a better agreement between theory and observations is obtained using a composite slab/2D turbulence (see, e.g., Bieber et al. 1996). Using various models of dynamical turbulence in conjunction with a slab turbulence spectral form that includes a dissipation range in QLT was shown to produce rigidity independent electron MFPs (see, e.g., Bieber et al. 1994;Schlickeiser 2002, 2003) in agreement with observations. However, it was reported by Dröge and Kartavykh (2009) that observed electron pitch angle distributions did not agree with those predicted by dynamical quasilinear theory. Figure 2 also shows three estimates for the parallel MFP calculated from analyses of three separate GLEs by Bieber et al. (2002),  and . Although corresponding to roughly the same rigidity, these values vary considerably, which may be due to inter-event variability in transport conditions. Some other events are, however, believed to involve transport in a closed interplanetary magnetic loop, for which energetic particles can have a much longer mean free path (Tranquille et al. 1987;Torsti et al. 2004) because of a greatly reduced amplitude of magnetic turbulence (Burlaga et al. 1981), and in particular slab turbulence (Leamon et al. 1998). Such events include those analysed by Ruffolo et al. (2006b), who report a parallel MFP of 1.2-2 au at 1.6 GV, and , who find a parallel MFP of 0.8 au at 3.1 GV (see also Sáiz et al. 2008). These values, however, have not been included in Fig. 2 as they may not apply to the interplanetary medium in general.
Considerably fewer observational values for the perpendicular MFP, shown in the bottom panel of Fig. 2, have been reported in the literature. The Palmer consensus range implies a constant value for this quantity over a rather broad range of rigidities at 1 au. However, it should be noted that Palmer (1982) reported a large spread of observed values for this quantity, as reported in the various studies considered in that paper (see also the discussion by Giacalone 1998). This can also be seen in the Jovian (Chenette et al. 1977;Vogt et al. 2020) and solar energetic electron (Dröge et al. 2014) perpendicular MFP observations seen in this figure, as well as in the values reported by Ferrando et al. (1993), Simpson et al. (1993). However, at higher rigidities the Burger et al. (2000) values imply a fairly rigidity-independent proton perpendicular MFP, while at lower energies, relevant to the Jovian electrons studied by Zhang et al. (2007), the authors reported a perpendicular MFP that increased with energy, with a rigidity dependence of ∼ P 0.1−0.7 .

Insights from Numerical Test Particle Simulations
Numerical test particle simulations, in general, involve solving the Newton-Lorentz equation for a large ensemble of test particles moving in the presence of a simulated background magnetic field (usually assumed to be uniform, but not always (see, e.g., Cohet and Marcowith 2016)) and synthetic turbulent magnetic fluctuations, generated under various assumptions that usually reflect the observed statistical properties of solar wind turbulence, but sometimes also with specific types of turbulence, in order to test the predictions of different scattering theories for specific turbulence scenarios. Such turbulence is usually generated as a superposition of waves with random phases from an assumed form of the turbulence power spectrum (for more detail, see, e.g., Decker and Vlahos 1986;Decker 1993;Giacalone and Jokipii 1994;Tautz 2012;Tautz and Dosch 2013), although alternative methods can be employed (e.g. Juneja et al. 1994;Mertsch 2020). Running diffusion and drift coefficients are then calculated based on the statistics of particle's displacement as a function of time; these can be subsequently compared with theory. These models are not limited to studies of diffusion coefficients alone, but have also been employed in testing the basic assumptions and inputs for various scattering theories, such as the forms that velocity autocorrelation functions assume at various levels of turbulence (Fraschetti and Giacalone 2012), or the validity of the use of the Corrsin hypothesis fundamental to many of these Fig. 3 Example of a sample particle trajectory through, and diffusion parameters calculated for pure slab turbulence superimposed on a uniform background field. The upper left hand panel shows the trajectory of a proton propagating through slab turbulence, with the background field along the z-axis. The remaining panels show the evolution of the different diffusion coefficients as a function of time. Note that coordinates x and y denote directions transverse to z, and that the unit of length in this figure is light-second theories (Tautz and Shalchi 2010;Snodin et al. 2013a). These models can be used to study the behaviour of particle transport coefficients corresponding to turbulence conditions in specific regions, such as those close to the Sun, where such insights are of particular value to studies of SEP transport (see, e.g., Laitinen et al. 2016;Laitinen and Dalla 2017;Chhiber et al. 2021b). Furthermore, results from these models can be used to estimate the behaviour of diffusion coefficients outside the heliosphere (see, e.g., Sonsrettee et al. 2015;Snodin et al. 2016b;Reichherzer et al. 2020). Figure 3 provides an example of outputs of such a code, showing a single particle trajectory, along with parallel and perpendicular diffusion coefficients and drift coefficients calculated for an ensemble of test particles, in the presence of synthetic axisymmetric magnetostatic slab turbulence. The turbulence was calculated using the methods outlined in Owens (1978) and Decker and Vlahos (1986), which generate a series of fluctuations δB n from a given power spectrum P k defined at discrete parallel wavenumber values k n = 2πn/L. The spectral form for P k was defined similarly to the spectra employed by Bieber et al. (1994) and Minnie et al. (2007a), with a wavenumber-independent energy-containing range, and a Kolmogorov inertial range. The minimum and maximum magnetic fluctuation wavelengths (λ min and λ max ) were fixed at λ min = 10 −4 au and λ max = 1 au, following . One such set of fluctuations is known as a single realisation of turbulence. A background magnetic field of B 0 = 5 nT along the z-axis was assumed, with the relative strength of the magnetic fluctuations being defined through (δB/B 0 ) 2 = 0.1. The trajectories of 2000 particles were integrated as they traversed the magnetic realisation. One set of diffusion and drift coefficients was then calculated from the average of the results of simulations performed with 25 different turbulence realisations, following the approach outlined by Minnie et al. (2007b). Each simulation is run for a total of 2000 correlation crossing times, which is the time it would take a particle of velocity V 0 to travel a distance equal to one correlation length λ c = 0.01 AU. All particles had the energy of 100 MeV with initial velocity directions uniformly distributed over a sphere. After an initial period where particles stream essentially freely in the z-direction, the proton parallel diffusion coefficient (top right panel of Fig. 3) reaches a diffusive limit at around 1.4 Ls 2 /s, where Ls denotes a light-second. The two perpendicular coefficients κ xx and κ yy are essentially identical, and do not reach a diffusive limit, a phenomenon associated with reduced-dimensional turbulence reported on by several studies (see discussion below). Drift coefficients (bottom right panel) fluctuate broadly around a running average corresponding to the weak scattering drift coefficient (black dashed line) (see Sect. 4 for more detail). Several such test particle models exist, of varying complexity, and the reader is invited to consult, e.g., , Mace et al. (2000), Tautz (2010), Dalena et al. (2012), , Lange et al. (2013), Arendt and Shalchi (2018), Mertsch (2020), for in-depth discussions of their workings, which are beyond the scope of the present study, which concerns itself rather with a brief review of some of the results of such simulations. Given the breadth of this particular field, a comprehensive review is not possible here. Giacalone and Jokipii (1994) studied perpendicular diffusion coefficients in the presence of a Kolmogorov spectrum in two and three dimensions, showing that perpendicular diffusion is suppressed in the presence of lower-dimensional 2D turbulence, as expected from earlier theoretical work (see . Shortly afterwards, Michałek and Ostrowsky (1996) studied the behaviour of parallel and momentum diffusion coefficients for various levels of isotropic Alfvénic turbulence, finding good agreement with the predictions of QLT for both coefficients at relatively low levels of turbulence. In a subsequent, comprehensive study,  investigated both parallel and perpendicular diffusion coefficients in the presence of Kolmogorov isotropic and compound turbulence, reporting that parallel mean free paths in the presence of simulated isotropic turbulence are several times smaller than those computed under the assumption of composite turbulence, and that the energy dependence of the simulated parallel mean free paths follow that expected from standard QLT, a finding also confirmed by the studies of Casse et al. (2001) and Candia and Roulet (2004). For perpendicular diffusion,  reported that their simulation results cannot be fully explained by the so-called field line random walk (FLRW) approximation. From their simulations of perpendicular diffusion coefficients in the presence of weak magnetostatic slab turbulence, Mace et al. (2000) drew a similar conclusion regarding FLRW applicability. Furthermore, these authors also found that the theory proposed by Bieber and Matthaeus (1997) could not completely reproduce their simulation results, although that theory yielded results in somewhat better agreement with simulations than the FLRW theory. It should be noted, however, that the FLRW description should be applicable for particles with gyroradii significantly smaller than the turbulence correlation scales , as in the case, for example, of SEPs. Further simulations, performed by  and employing composite 2D/slab turbulence generated using turbulence spectra that also included an energy-containing range, showed that the nonlinear guiding center (NLGC) theory introduced in that study yielded results in significantly better agreement with simulations over a broad range of energies compared with FLRW and the Bieber and Matthaeus (1997) theories. Qin et al. (2002a), however, in a study of perpendicular diffusion in the presence of slab turbulence, found that parallel scattering of the test particles suppressed perpendicular diffusion, leading to subdiffusive perpendicular transport, a result confirming the theoretical analysis of particle velocity autocorrelation functions performed by . The diffusive regime could be attained once more by the action of a modicum of 2D turbulence in the simulations, as demonstrated by Qin et al. (2002b). This result led to a modification of the NLGC theory by Shalchi (2006), motivated by the fact that the theory predicted a non-zero perpendicular diffusion coefficient in the presence of pure magnetostatic slab turbulence (see also . Ruffolo et al. (2004) investigated the behaviour of FLRW diffusion coefficients in composite turbulence, finding that magnetic field lines separate non-diffusively over relatively short distances, but diffusively over long distances. Further studies of FLRW coefficients in the presence of pure slab turbulence performed by Shalchi and Weinhorst (2009) have also shown that the behaviour of the slab spectrum energy-containing range can determine whether field line meandering is super-or subdiffusive. Minnie et al. (2007a) performed a comparative study of the parallel and perpendicular MFPs predicted by various theories with the results of simulations of the same in the presence of different levels of axisymmetric composite turbulence. These authors found that, for low levels of turbulence, the QLT accurately describes particle parallel MFPs, as well as their rigidity dependence (see also Hussein et al. 2015;Cohet and Marcowith 2016;Dundovic et al. 2020;Reichherzer et al. 2020). At higher turbulence levels, they reported that the weakly nonlinear theory, or WNLT (see Shalchi et al. 2004c) yields results in better agreement with simulations. In terms of the perpendicular MFP, Minnie et al. (2007a) reported that, at lower turbulence levels, FLRW provides a reasonable description of this quantity, while at higher turbulence levels, the NLGC result is in better agreement with simulated perpendicular MFPs. Furthermore, Qin et al. (2006) considered parallel diffusion in the presence of strong 2D turbulence, finding that such turbulence would influence the parallel diffusion of test particles, calling into question the prediction of QLT for such turbulence scenarios. Qin and Shalchi (2012) investigated the influence of larger-scale, energy range fluctuations on parallel and perpendicular mean free paths, finding that the assumed spectral index of this range has little effect on computed parallel diffusion coefficients, but a significant effect on perpendicular diffusion coefficients, as is expected from various nonlinear theories (see, however, the discussion by . Arendt and Shalchi (2018) considered the time-dependent behaviour of computed running diffusion coefficients. While finding good agreement with the unified nonlinear theory (UNLT, Shalchi 2010) in the case of perpendicular diffusion in the presence of composite turbulence, these authors reported that particle distribution functions as function of perpendicular position in the presence of slab and purely 2D turbulence are not well-approximated by Gaussians, which is a key assumption in many scattering theories. In a comprehensive study, Dundovic et al. (2020) studied the energy dependence of parallel and perpendicular diffusion coefficients in a wide variety of turbulence scenarios, studying in particular the energy dependence of the ratio of these quantities. These authors found that the diffusion coefficients are energy-dependent. Furthermore, Dundovic et al. (2020) reported that the second order QLT of Shalchi (2005) adequately describes their simulated parallel MFPs at intermediate levels of composite turbulence, and that the predictions of the UNLT are somewhat closer to their simulated perpendicular mean free paths than those of the NLGC.
Several studies have also addressed the behaviour of diffusion coefficients in the presence of turbulence with various specific properties.  studied the perpendicular diffusion of particles in the presence of nonaxisymmetric composite turbulence (see also Ruffolo et al. 2006a;Pommois et al. 2007;Tautz and Lerche 2011a;Strauss et al. 2016).
They also found that the NLGC theory provides a reasonably good description of the resulting diffusion coefficients, particularly with respect to their variation depending on the assumed level of 2D fluctuation anisotropy, as opposed to the predictions of the strongly nonlinear theory of Stawicki (2005). Hussein et al. (2015) investigated the influence of various turbulence models beyond the standard composite model, such as noisy slab turbulence, and found that the rigidity dependence of their computed diffusion coefficients was relatively insensitive to the choice of turbulence model assumed (see also Fraschetti 2016). Other studies have investigated the influence of magnetic helicity on diffusion coefficients (Tautz and Lerche 2011b), the influence of dynamical turbulence (e.g. Gammon et al. 2019), as well as the influence of intermittency. Pucci et al. (2016) reported that intermittency has a greater influence on parallel transport coefficients than on perpendicular coefficients.
Numerical simulations have also been employed to study pitch-angle diffusion coefficients. Qin and Shalchi (2009) compared the results of such simulations with the predictions of various nonlinear diffusion theories, finding that the second order QLT can reproduce simulation results, albeit for relatively weak slab turbulence, at the 90 • pitch angle, as opposed to standard QLT (see also Qin and Shalchi 2014a). For composite turbulence, these authors also showed that WNLT cannot describe the pitch angle diffusion coefficient accurately. These authors also reported difficulties in estimating these quantities at higher levels of turbulence, arising from the fact that for strong turbulence the test particles' pitch angles fluctuate too much to allow for a calculation of the pitch angle diffusion coefficient. However, several alternative techniques have been suggested to resolve these issues (e.g. Ivascenko et al. 2016;Pleumpreedaporn and Snodin 2019). Building on this prior work, Qin and Shalchi (2014b) explored the pitch angle dependence of the perpendicular diffusion coefficient, finding that for low levels of turbulence the FLRW theory predictions are accurate, but that it is better described by the UNLT at higher levels of turbulence.

Turbulence and Drift
Since their initial incorporation into numerical cosmic ray modulation models (see, e.g., Jokipii and Thomas 1981;Kota and Jokipii 1983;Potgieter and Moraal 1985;Burger and Potgieter 1989;Ferreira et al. 2003a;Alanko-Huotari et al. 2007;Pei et al. 2012;Ngobeni and Potgieter 2015;Kopp et al. 2017;Tomassetti et al. 2017;Song et al. 2021;) using various numerical techniques to model drift effects (see, e.g., Burger et al. 1985;Burger 2012;Mohlolo et al. 2022), drift of cosmic rays produced by gradients and curvature of the heliospheric magnetic field, as well as drift along the HCS, have been demonstrated to play a key role in charge sign-dependent cosmic ray modulation. Drifts have been demonstrated to be the cause of the 22-year cycle in cosmic ray intensities observed by neutron monitors and space-borne detectors (see, e.g., Webber et al. 1990;le Roux and Potgieter 1990;Reinecke and Potgieter 1994;Usoskin et al. 2001;Gieseler et al. 2017;Di Felice et al. 2017;Caballero-Lopez et al. 2019;Fu et al. 2021;Krainev et al. 2021;Mohlolo et al. 2022), accounting for the observed HCS tilt angle and magnetic polarity dependence of GCR intensities (McDonald et al. 1992;Lockwood and Webber 2005;, and to influence the observed GCR latitudinal gradients (e.g. Heber et al. 1996;Zhang 1997;Heber and Potgieter 2006). Drift effects may also play a role in GCR modulation in the heliosheath (e.g. Langner and Potgieter 2004;Potgieter and Langner 2005;Webber et al. 2008;Kóta 2016), in the transport of SEPs (e.g. Dalla et al. 2013;Battarbee et al. 2018;van den Berg et al. 2021), and in the very long-term modulation of GCRs (Moloto and Engelbrecht 2020).
It has been theoretically known for along time that the drift coefficient, which enters into the Parker (1965b) transport equation via the antisymmetric components of the diffusion tensor, is reduced from its maximal, weak-scattering value of κ ws A = vR L /3 (Forman et al. 1974) in the presence of magnetic turbulence (see, e.g. , Parker 1965a;Burger 1990;Jokipii 1993;Stawicki 2005; le Roux and Webb 2007). Numerical cosmic ray modulation studies also found that a reduction in the weak-scattering drift coefficient yields computed GCR intensities in better agreement with spacecraft observations (e.g. Burger et al. 2000), employing rigidity-dependent modifications to κ ws A constructed to yield good agreement between computed and observed GCR intensities (see, e.g., Vos and Potgieter 2016;Corti et al. 2019). A typical form for this reduction factor is (e.g. Burger et al. 2008) with P 0 = 1/ √ 2 GV a free parameter. It is interesting to note that several such studies have reported a possible solar cycle dependence in this drift reduction factor (e.g. Ndiitwani et al. 2005;Manuel et al. 2011;Aslam et al. 2019, which may be the result of solarcycle variations of the turbulent conditions in the heliosphere Moloto and Engelbrecht 2020;Engelbrecht and Wolmarans 2020;Burger et al. 2022).
Numerical test particle simulations have shed some light on the details of this effect.  and Candia and Roulet (2004) found that drift coefficients simulated for composite (slab/2D) and isotropic turbulent fluctuations in the presence of a uniform background magnetic field do decrease when turbulence levels increase. Minnie et al. (2007b) reported on the results of simulations computed for composite turbulence, assuming both a uniform background field and a magnetic field with a gradient. These authors confirmed the previous findings, and reported a rigidity dependence of the drift reduction factor. Furthermore, they reported that turbulence-reduced drift coefficients calculated with and without a gradient in the background field are very similar. Minnie et al. (2007b) investigated drift velocities in the presence of turbulence, and also reported a reduction due to turbulence, demonstrating that the drift motions of charged particles in such scenarios depend not only on the off-diagonal elements of the diffusion tensor, but also on those associated with the diffusion of these particles. Tautz and Shalchi (2012) performed a very broad range of simulations of particle drift coefficients for composite, isotropic, and pure slab turbulence. They found a reduction of said coefficients from the weak-scattering values with increased turbulence levels. The degree of reduction was found to depend on particle rigidity as well as on the spectral index of the energy-containing range of the assumed form of the 2D turbulence power spectrum of the simulated turbulent magnetic field fluctuations. Intriguingly, when pure slab turbulence is used, particle drift patterns are indeed affected, but their drift coefficients remain equal to the weak-scattering value regardless of the assumed turbulence level (see also Fig. 3), with the implication that such turbulence conditions would not lead to an effective reduction in ensemble drift effects.
Various theoretical approaches have been proposed to model the influence of turbulence on particle drifts. The focus here will be on those that yielded relatively simple, tractable results that were employed in GCR transport studies. As such, the more complicated approaches to this problem, such as those presented by Stawicki (2005) and le Roux and Webb (2007), will not be discussed here. One of the first tractable theoretical results for the drift reduction factor to be proposed was that of Bieber and Matthaeus (1997). These authors calculated the particle drift coefficient using the Taylor-Green-Kubo formalism, modelling the required velocity autocorrelation functions by assuming only moderately perturbed particle trajectories, and an exponential temporal decay of the said functions. Assuming that this temporal decay proceeds with a characteristic time τ that is a function of the FLRW perpendicular diffusion coefficient D ⊥ and the full particle speed v, Bieber and Matthaeus (1997) found that the drift coefficient can be written as with ω c the particle gyrofrequency in the unperturbed background field and ω c τ = 2R L /3D ⊥ . This expression, however, was found to disagree with the simulations of Minnie et al. (2007b), which prompted Burger and Visser (2010) to propose an expression derived from fits to those simulations, so that ω c τ = 11 Tautz and Shalchi (2012) also present an expression for the drift reduction factor using a fit to the results of their simulations. It is, however, unclear whether such fits would be applicable to turbulence conditions throughout the heliosphere, and both the Bieber and Matthaeus (1997) and Burger and Visser (2010) results would yield reduced drift coefficients in the presence of pure slab turbulence.  revisited the theory of Bieber and Matthaeus (1997), but argued that the decorrelation time associated with the perpendicular particle velocity autocorrelation function should rather be a function of the particle perpendicular velocity and perpendicular mean free path (an assumption in qualitative agreement with the results of numerical test particle simulations pertaining to this decorrelation time reported by Fraschetti and Giacalone (2012)). They found that which yields results in reasonable agreement with the simulations of Minnie et al. (2007b) and Tautz and Shalchi (2012), as the reduction factor would become unity in pure slab turbulence, due to the subdiffusive nature of perpendicular transport in this turbulence geometry (see Qin et al. 2002b;Shalchi 2006). Furthermore, Eq. (5) reproduces the expected solar cycle variation in the turbulence-reduced drift coefficient. It should be noted, however, that expressions for the turbulence-reduced drift coefficient derived following the Bieber and Matthaeus (1997) approach may not be applicable when turbulence levels are high, given the underlying assumption of relatively unperturbed particle trajectories (see also the discussion by van den Berg et al. 2021). Several studies investigated the spatial dependence of the various effects discussed above, which were found to vary considerably (see, e.g., Wiengarten et al. 2016), and subsequently were demonstrated to have a significant influence on CR transport (e.g. Engelbrecht and Burger 2015a), as well as on SEP transport (van den Berg et al. 2021). The behaviour of some of the turbulence-reduced drift coefficients discussed here will be the subject of a later section.

Theories of Particle Transport Coefficients
Next, we present a brief review of the development of various scattering and diffusion theories, followed by a comparative analysis of several transport coefficients derived from these theories at 1 au, and beyond.

Theoretical Approaches
Here we review common theoretical approaches used to estimate the diffusion coefficients in directions perpendicular to the large-scale magnetic field (κ ⊥ ) and parallel to the largescale magnetic field (κ ). This section will not consider the drift coefficient (κ A ), which was discussed in Sect. 4. While most theories attribute the diffusion of energetic charged particles to scattering from turbulent fluctuations in the magnetic field, and make use of models of magnetic turbulence, one simple framework used a theory of hard-sphere scattering together with gyration about a large-scale magnetic field (Axford 1965;Gleeson 1969). In this theory, κ = (1/3)v 2 τ , where τ is the mean free time that might be determined empirically, and there is a relationship between κ ⊥ and κ given by where ω is the gyrofrequency. If ωτ 1 (or in terms of a parallel mean free path λ and the Larmor radius R L , if λ /R L 1) then scattering is so frequent (relative to the gyrofrequency) that the gyromotion is irrelevant and κ ⊥ = κ . However, considering the measured values of the parallel mean free path as shown in Fig. 2, for energetic particles in the heliosphere the opposite is almost always true, λ /R L = ωτ 1, so that κ ⊥ κ . Various other theories, including some that are not based on hard-sphere scattering, have made use of the relationship in Eq. (6).
At about the same time, Jokipii (1966) proposed that the perpendicular diffusion of charged particles could be attributed to the random walk of magnetic field lines (i.e., the variation of the perpendicular coordinates (x, y) of a turbulent field line as a function of z), which became known as the field line random walk (FLRW) theory of perpendicular diffusion. In this view, one considers a field line diffusion coefficient that governs the field line displacement in, say, x as a function of z: Then consider that the particle moves along the field line with a constant speed v s = μv, which is treated as motion along the mean field. Then for a uniform distribution of μ between −1 and 1, one obtains the perpendicular diffusion coefficient along x as where D x can be derived from a theory of field line random walk for a given type of turbulence.
Note that the hard-sphere scattering model does not take the field line random walk into account, and the FLRW model does not take parallel scattering into account. The most successful modern theories of perpendicular diffusion have combined these two lines of thought by including the effects of both parallel scattering and FLRW on the perpendicular diffusion of particles, as will be described shortly. Now we will review the theories of diffusion coefficients of energetic charged particles, as well as the field line random walk, in a large-scale (mean) magnetic field B 0 (typically in the z-direction) with superimposed turbulent magnetic fluctuations so that It should be pointed out that in some theories b also varies with time based on the magnetic turbulence model. We will now review several common turbulence geometries used in studies of charged particle transport.

Slab Turbulence
Slab turbulence is an idealized model where the fields vary only in the direction along the mean magnetic field. To satisfy ∇ · B = 0, we must have b z = 0, i.e., slab fluctuations are transverse to the large-scale field. Such turbulence could be viewed as a broadband ensemble of Alfvén waves with a wide range of wavenumbers k z . The fluctuating field can be treated as a function of wave vector by subjecting it to a Fourier transform. Then the Fourier amplitude of a component b i is a random complex function b i (k z ). However, the squared magnitude |b i (k z )| 2 is proportional to the power spectral density S ii (k z ), which is defined as the Fourier transform of the autocorrelation function. The power spectral density is not a random quantity, but physically represents the magnetic fluctuation energy per wavenumber. Suppose an ensemble of turbulent fields has a specified functional form of the power spectrum and |b i (k z )| 2 . Then the Fourier amplitude b i (k z ) at each wavenumber involves a random complex phase, and each pattern of phases leads to a different representation of turbulence in real space. Note that spacecraft instruments typically measure the power spectral density of solar wind fluctuations as a function of frequency f . The latter can be considered proportional to the wavenumber based on the Taylor frozen-in hypothesis that states that temporal variations are caused by spatial variations as the medium convects past the observer (Taylor 1935). As such, these observations also provide the wavenumber power spectrum that is more relevant in the context of theoretical models. Jokipii (1966) initially developed the quasilinear theory (QLT), which remains the most commonly used theory of parallel diffusion. This theory considers resonant interactions between the particles' gyro-orbits with Larmor radius R L and turbulent magnetic fluctuations with the parallel wavenumber k z = 1/(|μ|R L ) (the resonance condition). This interaction leads to a diffusion coefficient D μμ in the pitch angle cosine μ given by (see also Bieber et al. 1994;Zank et al. 1998). This in turn can be related to the spatial diffusion coefficient for a nearly isotropic particle distribution as κ by (Jokipii 1968; Hasselmann and Wibberenz 1968) An important modification to QLT was the addition of dynamical turbulence. Bieber et al. (1994) considered a dynamical Lagrangian (as seen by the moving particle) correlation function of the turbulent magnetic fluctuations that retain the time coordinate dependence after a Fourier transform. They derived a broader resonance condition for pitch angle scattering, and obtained parallel mean free paths that were very different for electrons and protons at the same rigidity. Other notable modifications to QLT are the second-order QLT (Shalchi 2005) and the model that includes the effects of magnetic helicity , which, along with drift motions, can contribute to the explanation of the observed effects of the ∼22-year solar magnetic polarity cycle on the solar modulation of Galactic cosmic rays (Nuntiyakul et al. 2014).
Turning now to the perpendicular transport of cosmic rays, as mentioned above, the FLRW model of Jokipii (1966) requires knowledge of the diffusion coefficient for the magnetic field lines. Jokipii and Parker (1968) calculated the field line diffusion coefficient with a quasilinear model (not to be confused with the quasilinear theory of parallel particle diffusion), obtaining where b 2 ≡ b 2 and c is the correlation length of the fluctuating field. This result is now known to apply to fluctuations that are slab-like (see Isichenko 1991a) and is exact for slab turbulence. Jokipii and Parker (1969) then used that field line diffusion coefficient in the FLRW model of perpendicular particle diffusion. Note, however, that for pure slab turbulence, the perpendicular transport of particles is actually non-diffusive (see Sect. 8.1). The reason partially stems from a theorem by  and Jones et al. (1998) that for a magnetic field with an ignorable coordinate (i.e., for turbulence of reduced dimensionality in wavenumber space), the particles are confined to within a gyroradius of a magnetic flux surface and do not diffuse throughout space. For the case of slab turbulence, which has two ignorable coordinates, the theorem implies that particles remain within a gyroradius of a single field line. Consequently, parallel scattering of particles back and forth along the same field lines leads to a subdiffusive perpendicular transport. Nevertheless, the FLRW concept has been applied to other types of magnetic turbulence that has some variation in the perpendicular directions.

2D+Slab Turbulence
It has been shown that in the presence of a mean magnetic field, turbulence tends to develop wavenumber anisotropy with enhanced power along the perpendicular wavenumbers (k x , k y ) (Shebalin et al. 1983). An extreme type of turbulence that could result from such a process is 2D turbulence, an idealized case of transverse fluctuations that depend only on the coordinates perpendicular to the mean field: Therefore, slab turbulence is dynamically unstable in the sense that interactions among the wave modes tend to cause the slab power along the k z -axis in wavenumber space to spread in the k x and k y directions and become more like 2D turbulence. Mathematically, we can express this 2D field in terms of a scalar potential a(x, y): Then the magnetic fluctuation follows contours of constant potential, and when combined with a mean field B 0ẑ , the magnetic field lines are confined to a cylindrical surface bounded by a constant a(x, y). Thus a 2D model of turbulent fluctuations gives rise to a flux-tube structure, as in the "spaghetti" model frequently used to describe solar wind plasma (e.g., Bartley et al. 1966;Bruno et al. 2001;Borovsky 2008).  identified that the distribution of turbulent power in wavenumber space (k , k ⊥ ) resembles a Maltese cross, with relatively more power near the k z axis (as in slab turbulence) and the (k x , k y ) plane (2D turbulence). This led to the idealized 2D+slab two-component model of magnetic fluctuations in the solar wind: where as before both fluctuation components are transverse to the mean field. A physical interpretation is that the slab component represents "fossil" Alfvénic fluctuations from the Sun, which can interact to form the 2D turbulence with a flux-tube structure. Later observational research on solar wind turbulence corroborated the Maltese cross picture (Weygand et al. 2009) and clarified how the relative contributions of slab-like and 2D-like components are different for slow and fast solar wind (Dasso et al. 2005) and may vary with distance from the Sun (Bandyopadhyay and McComas 2021). The idealized 2D+slab prescription is particularly convenient for numerical modeling because generating representations of such fields only requires 1D and 2D inverse fast Fourier transforms, whereas a general fluctuation model requires a 3D transform that may be computationally expensive and requires a large amount of memory. This model has been widely employed in models of turbulence transport, as noted in Sect. 2.
In the QLT framework, pitch angle scattering is attributed to a resonant value of k z , in which case the parallel diffusion of energetic particles is affected only by the slab component. Bieber et al. (1994) noted a long-standing discrepancy between the predictions of QLT (assuming slab turbulence with the observed solar wind fluctuation amplitude) and observed mean free paths (as shown in Fig. 2; see also Palmer 1982). They then considered that only 20% of the turbulent energy of the solar wind is in slab turbulence, with the remaining 80% in 2D turbulence, so the turbulent energy relevant to parallel diffusion is only 20% of the value previously used, thereby resolving the discrepancy. This slab energy fraction was later corroborated by observational tests (Bieber et al. 1996). The QLT of parallel diffusion in 2D+slab turbulence was extended by adding nonlinear effects and coupling with perpendicular diffusion in the weakly nonlinear theory (WNLT; Shalchi et al. 2004c). Qin et al. (2006) also demonstrated computationally that a strong 2D turbulence component does affect the parallel scattering in a mixture of slab and 2D turbulence.
Theories of perpendicular diffusion in 2D+slab turbulence have accounted for the magnetic field line diffusion, which in turn was first addressed by Matthaeus et al. (1995). This work employed Corrsin's independence hypothesis (Corrsin 1959) with a diffusive model of the field line spread and decorrelation, which later work has referred to as the diffusive decorrelation (DD) approach. For the 2D component of magnetic turbulence, they obtained the field line diffusion coefficient (as expressed by Ruffolo et al. 2004;), whereλ is a newly defined "ultrascale" of turbulence, related to the k −2 ⊥ moment of the 2D power spectrum. Their result has a form of Bohm diffusion, i.e., a linear dependence on b/B 0 , in contrast with the quadratic dependence of the quasilinear result for slab turbulence (Eq. (13)). In this framework, the contributions to the field line diffusion coefficient from the slab component (Eq. (13)) and 2D component (Eq. (17)) are combined into an overall field line diffusion coefficient D as follows: Matthaeus et al. (1995) noted that when combining the slab and 2D fluctuation components, the field lines are no longer confined to flux surfaces and are able to wander in three dimensions. Ruffolo et al. (2004) extended that work to examine separation between neighboring field lines, which is related to a spreading (as opposed to displacement) of a particle distribution. The effect of nonaxisymmetry, i.e., different turbulence properties along the k x and k y directions on the field line diffusion coefficients D x and D y was examined for 2D+slab turbulence by Ruffolo et al. (2006a) and for other types of turbulence by Tautz and Lerche (2011a).  and Shalchi and Weinhorst (2009) identified conditions for non-diffusive behavior.
Another theoretical development regarding field line diffusion coefficients was to use Corrsin's hypothesis with a different assumption for field line spreading and decorrelation, known as random ballistic decorrelation (RBD) (Ghilea et al. 2011). The resulting diffusion coefficient for the 2D component of turbulence is where λ 2D is the perpendicular correlation length of the 2D turbulence, and in this case the slab and 2D contributions are simply added to obtain the total field line diffusion coefficient. Equation (19) indicates another pathway to Bohm diffusion that relates to the general arguments by Kadomtsev and Pogutse (1979). The Bohm diffusion results (with a linear dependence of D on b/B 0 ) in the theories of Matthaeus et al. (1995) and Ghilea et al. (2011), based on Corrsin's hypothesis, are in contrast with other work that predicts percolative or trapping behavior in a quasi-2D limit (Gruzinov et al. 1990;Isichenko 1991b;Vlad et al. 1998;Neuer and Spatschek 2006;Negrea et al. 2007). Indeed, simulation results by Ghilea et al. (2011) do show deviations from Eq. (19), which are attributed to trapping effects, when the slab fraction is very low. Field line diffusion coefficients from the DD and RBD approaches, together with an older self-consistent ODE approach that also relies on Corrsin's hypothesis (e.g., Saffman 1963;Taylor and McNamara 1971;Shalchi and Kourakis 2007), were compared by Snodin et al. (2013a), while the pre-diffusive evolution of the field line random walk was explored by Snodin et al. (2016a). In parallel, there were theoretical developments concerning how the field line diffusion coefficient should be incorporated into a theory of perpendicular diffusion of particles. Bieber and Matthaeus (1997) used Eq. (6) as in the classic hard-sphere scattering theory, but with the correction that the scattering time τ was no longer related to parallel scattering, but rather to the field line diffusion coefficient. Then in collaboration with other authors, they developed the nonlinear guiding center theory (NLGC; , which abandoned the framework of Eq. (6) in favor of a theory based on Corrsin's hypothesis that effectively included both parallel scattering and the field line random walk, and could accommodate dynamical turbulence as well. This theory uses a DD approach, in which the diffusion coefficient is determined using Corrsin's hypothesis for a probability distribution function based on the diffusion itself; therefore, it obtains an implicit equation for κ ⊥ with an integral over wavenumber of an integrand that involves κ ⊥ : where a 2 = 1/3, S xx is the power spectrum of b x (axisymmetry in the x-y plane is assumed), T is the effective parallel scattering time, and γ is the decay rate of the dynamical correlation spectrum.
Here particles are not assumed to strictly follow one field line; rather, their perpendicular diffusion coefficient is determined based on the mean free time for the decorrelation of the guiding center motion along x and y, which could be due to either parallel scattering or changes in the field line direction. The physical relationship between the NLGC theory and the FLRW theory of perpendicular diffusion is as follows . Particles may initially follow field lines, spreading roughly at the rate expected from FLRW theory. Subsequently parallel scattering causes particles to backtrack along the same field lines so their random walk has a "negative memory," leading to a period of transient subdiffusion (see Sect. 8). Finally, the particle motion decouples from the original field line, leading to asymptotic diffusion at long times that can be modeled by the NLGC theory.
The NLGC theory was originally tested for magnetostatic and dynamical 2D+slab turbulence ), but it can accommodate a general transverse fluctuation field. Stawicki (2005) proposed an alternative, strongly nonlinear (SNL) theory, and suggested that the NLGC and SNL should be tested by comparing their very different results for nonaxisymmetric turbulence. As a result,  worked out the NLGC theory for nonaxisymmetric turbulence and performed computer simulations in nonaxisymmetric representations of 2D+slab turbulence, which were inconsistent with SNL and reasonably consistent with NLGC. Meanwhile, Shalchi (2006) proposed a modification to NLGC to remove the direct slab contribution in the case of 2D+slab turbulence. A wide variety of possible modifications to NLGC theory, as well as other nonlinear perpendicular transport theories, were described by Shalchi (2009).
To conclude this subsection, we discuss two more modifications of NLGC that can be considered state of the art in NLGC-inspired theories. The unified nonlinear theory (UNLT) of Shalchi (2010) again uses Corrsin's hypothesis with a DD approach in which the diffusion coefficient is considered to be proportional to |μ|, yielding an implicit equation for κ ⊥ , which uses Eq. (20) and This theory provides sensible results in some limiting cases where NLGC does not, e.g., for pure slab turbulence. The other notable modification to NLGC was to implement an RBD interpretation with a backtracking correction (also making use of the modification by Shalchi 2006). This RBD/NLGC theory provides an explicit formula for κ ⊥ , using Eq. (20) with This explicit formula is more convenient for numerical simulation work than the implicit formulas from most versions of NLGC, and resulting κ ⊥ values were shown to provide a better match to simulation results than the original DD-based NLGC ). This formulation was extended to nonaxisymmetric turbulence by Strauss et al. (2016), and evaluated for different models of dynamical turbulence by Dempers and Engelbrecht (2020). The numerical results obtained for κ ⊥ from the FLRW, NLGC, UNLT, and RBD/NLGC theories are compared in Sect. 5.2 and Fig. 4 below.

Modified 2D or Slab Turbulence
We must briefly mention some generalizations of the idealized 2D and slab turbulence models that make them multidimensional. Reduced magnetohydrodynamics (RMHD) is a classic dynamical model that includes transverse turbulence with a slow variation parallel to the mean magnetic field, developed to describe tokamak plasmas (Strauss 1976) and later widely employed to model fluctuations in solar coronal magnetic loops (Longcope and Sudan 1994).  introduced an idealized magnetostatic model called noisy reduced magnetohydrodynamic (NRMHD) turbulence, with a prescribed power spectrum. While 2D turbulence only has power exactly on the (k x , k y ) plane, NRMHD turbulence spreads this over a finite volume bounded by |k z | ≤ K, making it a three-dimensional model. The "noisy" aspect is that by increasing K, the model can be make more slab-like (i.e., adding noise to RMHD turbulence), whereas the limit K → 0 yields 2D turbulence. The magnetic field line diffusion coefficients indeed display quasilinear behavior for large K and Bohm behavior for small K, as verified by computer simulations for both synthetic NRMHD and dynamical RMHD turbulence (Snodin et al. 2013b).
Later the time-differenced field line random walk in NRMHD was explored by . For particle transport, this magnetic fluctuation model has the interesting property of a hard resonance gap, in which particles with a wide range of pitch angles have no resonant interactions because the fluctuation power spectrum is zero at |k z | > K. This allows the effects of nonlinear scattering to be seen more clearly, as in computer simulations by Shalchi and Hussein (2014) and . An analogous "noisy slab" model was introduced by Shalchi (2015), in which the slab power along the k z -axis is extended over the volume of a cylinder about that axis, again making the model three-dimensional and allowing control over the degree to which it approximates slab turbulence.

Isotropic Turbulence
Finally we will say a few words about magnetic field line and particle diffusion in isotropic turbulence. These have been used to model turbulent magnetic fields in the Galaxy. In particular, the case of zero mean field (B 0 = 0) is the only truly isotropic case, with no preferred directions and no distinction between parallel and perpendicular diffusion, which may apply to highly chaotic magnetic fields that lack an organizing large-scale field. For this reason, and because of its conceptual simplicity, historically particle diffusion in isotropic turbulence received extensive attention even for cosmic ray transport in the heliosphere. Quasi-linear theory for isotropic magnetosonic turbulence was developed by Schlickeiser and Miller (1998). Because isotropic turbulence is compressive, an important contribution to scattering is provided by the Landau resonance where a particle becomes trapped in a region of weaker magnetic field when its velocity is approximately equal to the phase speed of the wave, a process also known as transit time damping (Fisk 1976a). This process is most effective at small values of the pitch angle cosine. Schlickeiser and Miller (1998) computed the pitch angle scattering coefficient for a mixture of slab and isotropic modes and concluded that the resulting parallel diffusion coefficient could be significantly larger than in the case of pure slab turbulence.
Theories of the field line random walk, based on Corrsin's hypothesis, were presented by Sonsrettee et al. (2015Sonsrettee et al. ( , 2016. From numerical simulation of particle transport, it is generally agreed that there are various regimes (e.g., for varying particle rigidity) in which different theoretical formulations are appropriate (Casse et al. 2001;Candia and Roulet 2004;Snodin et al. 2016b;Reichherzer et al. 2020). For example, Casse et al. (2001) concluded that in some cases quasilinear theory (QLT) works well, but in other cases the field line random walk needs to be taken into account. For the case of zero mean field,  found good agreement between numerical results and a quasilinear approach for very high particle rigidity (with the unperturbed orbit as a straight line) and a modified version of QLT for low rigidity. Gammon et al. (2019) used a theoretical approach that included dynamical turbulence effects.

Compressive Turbulence
Relatively fewer results for transport coefficients for the case of compressive turbulence, where turbulent fluctuations occur along the background field, have been published. Given the predominance of this type of turbulence in the heliosheath (see, e.g., Burlaga et al. 2014Burlaga et al. , 2018, such results would be particularly relevant to the self-consistent modelling of particle transport in this region. Lerche and Schlickeiser (2001) and Schlickeiser (2002) investigated the influence of fast magnetosonic waves and different spectral forms on pitchangle diffusion coefficients derived using the QLT approach, finding that the mean free path due to these fast waves, assuming a single power-law form for their spectrum that commences at wavenumber k m,f , with a spectral index q f , is given by (Schlickeiser 2002) with v A the Alfvén speed, and δB 2 f the variance of the fast magnetosonic fluctuations. Yan and Lazarian (2004) extend this analysis, additionally considering Alfvénic and slow magnetosonic modes (see also Yan and Lazarian 2002) in their calculation of Fokker-Planck transport coefficients, and reporting that fast modes provide the dominant contribution towards cosmic ray scattering (see also Yan and Lazarian 2008). It should be noted that these latter studies were done in the context of the galactic transport of cosmic rays.
Motivated by the need to model the influence of compressive turbulence on the perpendicular diffusion coefficient in particular, given the importance of this quantity to particle transport in the outer heliosphere, Strauss et al. (2016) modeled the transverse displacement of a particle using the well-known TGK approach, thereby writing the instantaneous perpendicular diffusion coefficient D ⊥ as a function of the gyroperiod-averaged perpendicular particle velocity as well as a perpendicular velocity autocorrelation function decorrelation time following an approach similar to that employed by, e.g., Bieber and Matthaeus (1997). Strauss et al. (2016) then proceeded to approximate the gyroaveraged perpendicular speed using the gyroperiod-averaged expression for this quantity (see Rossi and Olbert 1970), in an approach similar to that employed by Fraschetti and Jokipii (2011), and employing a simple expansion of the magnetic field into a uniform and a 3D fluctuating component, assuming also that particles encounter decorrelated turbulence by crossing one (perpendicular) correlation length l x transverse to the background field. This yields D x (μ) ≈ l x g V 2 x 1/2 and D y (μ) ≈ l y g V 2 y 1/2 , where g V denotes the relevant component of the ensemble-averaged gyroperiod-averaged velocity, with x and y denoting directions perpendicular to the uniform background field B 0 , assumed to lie in the z-direction. Assuming predominantly compressive fluctuations (without specifying the physical mechanism which gave rise to them), as well as axisymmetry, would then result in with r L being the maximal Larmor radius, and δB 2 z the variance associated with the compressive fluctuations. It is interesting to note that the pitch angle dependence of the above expression was found by Strauss and Fichtner (2014) to qualitatively reproduce the anisotropy in energetic particle distributions observed by Voyager 1 beyond the heliopause (see Krimigis et al. 2013).

Results for Transport Coefficients
Many studies have explored the expected behaviour of various theoretical expressions for cosmic ray diffusion MFPs in the heliosphere, both spatially and temporally. Some studies used turbulence quantities as direct inputs for expressions derived from various scattering theories (see, e.g., Erdős and Balogh 2005;Engelbrecht and Wolmarans 2020), while others assumed pre-specified models for large scale heliospheric plasma quantities, such as a Parker heliospheric magnetic field, and employ a suitable turbulence transport model to calculate the relevant turbulence quantities. This is the approach taken by e.g. Zank et al. (1998), Pei et al. (2010), Engelbrecht and Burger (2013a,b), Engelbrecht and Burger (2015a,b), Zhao et al. (2017), and Adhikari et al. (2021a). In a more recent development outputs from MHD codes are used to model the behaviour of large scale heliospheric plasma quantities. These models can be coupled with a turbulence transport model to produce outputs which can be used to calculate GCR transport parameters. This approach was taken for different MHD codes and turbulence transport models, by Wiengarten et al. (2016) and Chhiber et al. ( , 2021b. Nevertheless, the above studies report considerable variation in the diffusion, and turbulence-reduced drift, coefficients obtained using various scattering theories, finding also considerable variation due to assumptions made regarding the values and behaviour of the turbulence quantities used to calculate these parameters. To demonstrate this, several expressions for the GCR proton perpendicular MFPs, derived assuming different scattering theories, and used in GCR transport studies will be compared below, employing turbulence quantities yielded by two different turbulence transport models. Note that the list of expressions described below is by no means comprehensive, and that various models, derived under different assumptions about turbulence conditions as well as using multiple scattering theories, can be found in the literature (see, e.g., Pei et al. 2010;Engelbrecht and Burger 2013a;Qin and Zhang 2014;Shalchi 2015;Gammon and Shalchi 2017;Shalchi 2018). In what is to follow, expressions derived from the field line random walk theory of Jokipii (1966) and Matthaeus et al. (1995), the nonlinear guiding centre (NLGC) theory of , the unified nonlinear theory (UNLT) of Shalchi (2010), and the random ballistic decorrelation interpretation of NLGC (RBD/NLGC) proposed by  will be compared.
The 2D turbulence power spectrum is a key input for most theories of perpendicular diffusion, and therefore its form, particularly at low wavenumbers, has a significant influence on the resulting perpendicular MFPs (see, e.g., Engelbrecht 2019b). The UNLT, RBD/NLGC, and FLRW perpendicular MFPs discussed here are derived assuming a piecewise-defined, continuous 2D spectral form with an energy-containing range displaying an observationally-motivated k −1 ⊥ wavenumber dependence (e.g. Goldstein and Roberts 1999), and an inertial range with spectral index −ν commencing at lengthscale λ 2D . In what follows, a Kolmogorov value of 5/3 is assumed for this quantity. This (modal) spectral form, given by (Engelbrecht and Burger 2015a) also has a third, 'inner' range at lengthscales below λ o . This is specified, following the arguments of , to allow for the calculation of a non-diverging 2D ultrascale (see also Matthaeus et al. 1999), with a spectral index q = 3, to ensure that the 2D magnetic fluctuations satisfy the solenoidal condition ). Although λ 2D is approximately of the same order as the magnetic correlation scale (see, e.g., Bruno and Carbone 2016), and is therefore approximated thus both here and in other studies, the lengthscale λ o at which the energy-containing range commences remains a significant source of uncertainty, exacerbated by the fact that the perpendicular MFPs derived from many scattering theories are extremely sensitive to the behaviour of this quantity (see, e.g., Engelbrecht and Burger 2015b;Engelbrecht 2019b). This is due to the fact that no direct observations of this quantity exist, which may be attributed to the interference of coherent structures with associated timescales of around a month in magnetic observations (Goldstein and Roberts 1999), and that single-point measurements are not well-suited to characterizing the spectral properties of magnetic fluctuations on such large scales ). As such, several indirect estimates for this quantity have been proposed, based on GCR modulation (Engelbrecht and Burger 2013a), the solar rotation rate (Adhikari et al. 2017a), and calculation from observed magnetic island sizes (Engelbrecht 2019b). The modulationbased estimate is employed here, with λ 0 = 12.5λ 2D , which also ensures a well-defined energy range at all radial distances. It should be noted that it is still not entirely clear how this quantity should be modelled. The normalization constant C 0 is determined using ∞ 0 2πk ⊥ S 2D (k ⊥ ) dk ⊥ = δB 2 2D (see, e.g., , where δB 2 2D is the 2D variance. For the spectral form given in Eq. (29) this yields (Engelbrecht and Burger 2015a) The 2D ultrascale can be calculated in the manner described by , which yields (e.g. Engelbrecht and Burger 2015a;) which enters into the FLRW expression for the perpendicular MFP, as well as the turbulencereduced drift coefficient proposed by Burger and Visser (2010). In brief, Engelbrecht and Burger (2015a) derive an expression for the perpendicular MFP from the UNLT, using Eq. (29): which needs to be evaluated numerically. Note that Eq. (32) is approximate, neglecting terms resulting from the integration of the inner and inertial ranges of Eq. (29) that make very small contributions to the perpendicular MFP (see also Shalchi 2013Shalchi , 2014, and implicit, in the sense that the RHS of Eq. (32) is a function of λ ⊥ . Dempers and Engelbrecht (2020) derived an expression for λ ⊥ from the RBD/NLGC theory given by (see also Strauss et al. 2016) with δB 2 T denoting the total (slab+2D) transverse variance, erfc the complementary error function, and where Equation (33) is also approximate, for the same reasons as the UNLT expression. The corresponding FLRW expression is given by (e.g. Matthaeus et al. 1995) with the 2D ultrascale as given in Eq. (31). The NLGC expression to be evaluated here is that derived by Shalchi et al. (2004a), as rewritten by Burger et al. (2008) for a general slab/2D ratio: with denoting the standard Gamma function. This expression is derived for a different form of the 2D turbulence power spectrum, which only displays a wavenumber-independent energy containing range and an inertial range with spectral index 2ν, and is included here for comparative purposes, as it has been used in many GCR transport studies. It should also be noted that this expression is similar to that proposed by Zank et al. (2004) and employed by e.g. Florinski and Pogorelov (2009), Zhao et al. (2017, Guo et al. (2021). Note that in what follows it is assumed that the parameter a, arising as an ansatz of the NLGC theory and interpreted as the probability that particles follow field lines (e.g. Shalchi 2010) or as describing the finite gyroradius effects of turbulence on the particle trajectories (Shalchi 2016), is assumed to take a value of 1/ √ 3, based on fits to the numerical test particle simulations of . Theoretical calculations of this quantity vary, ranging from a 2 ∈ [1, 2] from a derivation of the NLGC theory directly from the Newton-Lorentz equation proposed by Shalchi and Dosch (2008) to a 2 ≥ 1 (Shalchi 2016) and a 2 = 1/3 (Shalchi 2019) (see also the discussion by Arendt and Shalchi 2020;Shalchi 2021).
The perpendicular MFP expressions described above, as well as the drift reduction factor due to its dependence on λ ⊥ , require as a basic input some expression for the parallel MFP. The reason is that in some of the above theories the particle's pitch angle autocorrelation functions decorrelate after a time τ = λ /v (see, e.g. Engelbrecht 2019a). This assumption is not always valid (see Shalchi 2011) and is not directly made for the UNLT (Shalchi 2010(Shalchi , 2020. The parallel MFP expression used here was constructed by Burger et al. (2008) from the expressions derived by Teufel and Schlickeiser (2003) from QLT (Jokipii 1966), assuming a slab turbulence spectrum with a wavenumber-independent energy-containing range, and an inertial range with spectral index −s, given by where s = 5/3, R = R L k m , R L the maximal Larmor radius, and k m = 1/λ sl the wavenumber at which the slab spectrum inertial range commences (see, e.g., Zank et al. 1998;le Roux et al. 1999, for similar, parametrised QLT expressions). This expression has often been used in modulation studies (e.g. Engelbrecht and Burger 2013a; Moloto and Engelbrecht 2020), as it yields the rigidity dependence expected from what is reported for protons at higher energies from the consensus parallel MFP values (see Sect. 3.1). Furthermore, from numerical test particle simulations, the QLT yields results in reasonably good agreement with simulated parallel MFPs for relatively low levels of turbulence (see Sect. 3.2). However, it must be noted that QLT may not provide the best description of particle parallel transport, owing to limitations inherent in the theory (see, e.g., Shalchi 2009, and references therein). Several modifications to this theory, and a number of alternatives, were proposed in the literature (see, e.g., Shalchi et al. 2004c;Shalchi 2005). The results of these theories, however, are relatively intractable, and therefore Eq. (37) is used here to model λ . For the purpose of comparison, we also consider two expressions for the electron parallel MFP, constructed from expressions derived from the QLT by Teufel and Schlickeiser (2003) for a slab turbulence power spectrum with a dissipation range with spectral index p, assuming the random sweeping and damping models of dynamical turbulence (see, e.g., Bieber et al. 1994;Hussein and Shalchi 2016, for more detail). The former expression is given by where α d ∈ [0, 1] is a parameter governing the strength of dynamical effects (Teufel and Schlickeiser 2003), and v and v A the particle and Alfvén speeds, respectively. The electron parallel MFP derived under the assumption of the damping model of dynamical turbulence is given by where f 1 = 2/(p − 2) + 2/(2 − s). Here, it is assumed that p = 2.6, based on the observations at 1 au reported by Smith et al. (2006). We also use the 1 au value for k d reported by Leamon et al. (2000), although this quantity is expected to vary spatially (see, e.g., Smith et al. 2012;Bruno and Trenchi 2014;Engelbrecht and Strauss 2018;Duan et al. 2020).
The top left panel of Fig. 4 shows both proton and electron QLT parallel MFPs, along with the Palmer (1982) consensus range, to guide the eye. The proton parallel MFP displays, at lower rigidities, the expected P 1/3 dependence, which steepens to a P 2 dependence at the highest rigidities shown, this latter dependence being determined by the form assumed for the slab turbulence power spectrum at low wavenumbers. At the lowest rigidities, both electron parallel MFPs display a flatter rigidity dependence, with the damping turbulence expression yielding values almost an order of magnitude smaller than those yielded by the random sweeping expression. Comparison of this figure with Fig. 2 leads to a conclusion that the random sweeping expression probably yields unrealistically large values for λ at 1 au. The proton perpendicular mean free paths (shown in the top right panel of Fig. 4) yielded by the different theories considered here display a range of rigidity dependences, ranging from the rigidity independent FLRW result to varying rates of monotonic increase with increasing rigidity, with the UNLT and RBD/NLGC expressions becoming approximately constant with increasing rigidity beyond ∼ 1 GV, and the NLGC expression increasing, due to its λ 1/3 dependence. Values for λ ⊥ yielded by the UNLT (largest) and FLRW theory (smallest) differ by almost an order of magnitude at 1 GV, suggesting that modulation studies using these different expressions would in most probability yield very different results. The bottom left panel of Fig. 4 shows turbulence-reduced drift lengthscales calculated using the perpendicular MFPs shown in the top left panel of this figure (from Eq. (5)), along with the ad hoc drift reduction factor commonly used in modulation studies (Eq. (3)), and the maximal Larmor radius, the weak-scattering value for this quantity. At high rigidities, all reduced drift scales converge on this quantity, with the implication that turbulent HMF fluctuations would not influence the drifts of such high-energy particles. At lower rigidities, the theoretical expressions considered here for λ ⊥ yield very similar turbulence-reduced drift scales, which in turn do not greatly differ from the values yielded by the ad hoc reduced drift scale. It must be noted, however, that due to a lack of any spatial dependence, the ad hoc reduction factor in Eq.
(3) will be the same throughout the heliosphere, with the implication that it could lead to too great a reduction in drift effects in the outer heliosphere, where the turbulence levels are lower.
To demonstrate the influence of the assumed TTM on transport coefficients, the bottom right panel of Fig. 4 shows the QLT parallel (Eq. (37), blue lines) and RBD/NLGC perpendicular MFPs (Eq. (33), red lines) at a rigidity of 1 GV as a function of heliocentric radial distance, calculated using turbulence quantities computed with the two-component Oughton et al. (2011) (dashed lines),  (solid lines) and Usmanov et al. (2016Usmanov et al. ( , 2018 (dotted lines) TTMs. Note that the turbulence quantities shown in Fig. 1 are employed here, including the influence of pickup ion-generated waves. The changes wrought by this process on the large wavenumber behaviour of the slab turbulence power spectrum (see, e.g., Williams and Zank 1994;Zank 1999;Isenberg 2005;Cannon et al. 2014) have not been taken into consideration in the derivation of Eq. (37), and hence this expression may be an oversimplification (see Engelbrecht 2017, for more detail). Considerable differences can be seen in the resulting MFPs, with perpendicular MFPs differing by almost two orders of magnitude at 85 au, and parallel MFPs similarly at ∼ 4 au. Both sets of pathlengths display differing radial dependence, too, depending in part on how the influence of pickup ion-generated waves is modelled in these TTMs. As the pickup ion-generated waves only influence the slab spectrum at higher wavenumbers, it would be expected that these fluctuations would not greatly affect the scattering of GCRs, these being influenced mainly by the behaviour of the turbulence power spectra at lower wavenumbers corresponding to the energy-containing and inertial ranges. As such, the omission of pickup ion effects would lead to relatively monotonic increases/decreases in magnetic correlation scales/variances (see, e.g., Breech et al. 2008;Usmanov et al. 2012;Engelbrecht and Burger 2015a;Adhikari et al. 2017b). The behaviour of the transport coefficients in such a scenario would be somewhat similar (see, e.g., Engelbrecht and Burger 2015b) to that seen in Fig. 4 for the Usmanov et al. (2016Usmanov et al. ( , 2018 model outputs, although the outputs of the latter model in that figure were calculated including the influence of pickup ions and hence are not exactly the same as when pickup ion effects are ignored. This is reflected in the increase in the parallel MFP, due to the decrease of the slab variance at higher radial distances. The perpendicular MFP, however, remains relatively constant as a function of radial distance, because the decrease in the 2D variance is balanced out to some degree by the increase in the parallel MFP.
In concluding this section we would like to point out that despite all these developments in turbulence transport theories many contemporary GCR modulation studies that will be discussed below employed, with success, MFPs calculated using simple power law scalings for the magnetic variances and the respective correlation lengths.

Ab Initio Studies of Particle Transport in the Heliosphere
With the recent developments in our understanding of turbulence and its influence on the transport of energetic charged particles in the heliosphere, many recent studies implemented these advances so some degree, in order to address a variety of challenges to our understanding of particle transport and GCR modulation. Given the complexities implicit to solving a TTM in conjunction with a particle transport model, many studies opted to model turbulence quantities using relatively simple, and hence tractable, analytical scalings fitted to spacecraft observations of such quantities as magnetic variances that directly influence particle diffusion and drift coefficients. Burger et al. (2008) followed such an approach, using QLT and NLGC theory parallel and perpendicular MFP expressions derived by Teufel and Schlickeiser (2003) and Shalchi et al. (2004a), to demonstrate that a Fisk-type HMF model can lead to a linear relationship between GCR proton relative amplitudes and latitude gradients in qualitative agreement with observations. This model was also used by Hitge and Burger (2010) to investigate the influence of various models for Fisk-type fields on GCR 26day recurrent variations. In a similar study of the transport of GCR electrons that took into account the influence of dissipation range slab turbulence Engelbrecht and Burger (2010) reported that such small scale turbulence quantities can have a large influence on the global transport of these particles. Qin and Shen (2017) also modeled basic turbulence quantities using relatively simple scalings combined with direct spacecraft observations of large scale heliospheric plasma quantities, and used these as inputs for diffusion coefficients proposed by Qin and Zhang (2014). Their computed GCR proton intensities were in good agreement with PAMELA observations taken from 2006-2009 at Earth.
A simplified approach to modelling turbulent quantities using fits to spacecraft observations was also proven useful when studying time-dependent GCR transport. Shen and Qin (2018) and Shen et al. (2019) extended the model of Qin and Shen (2017), with the former study employing time-dependent large scale plasma quantities to reproduce the observed GCR intensities over the last several solar cycles as well as along the Voyager and Ulysses spacecraft trajectories, and the latter study demonstrating that their model can reproduce the observed intensities of multiple different species of GCR. Moloto et al. (2018) studied GCR modulation during the solar minima corresponding to the start of solar cycles 22-24, using as inputs the diffusion coefficients employed by Burger et al. (2008), as well as the turbulence-reduced drift coefficient proposed by , based on direct observations of turbulence quantities during these periods that have been spatially scaled following the results of the Oughton et al. (2011) TTM; they were successful in reproducing spacecraft observations during all three solar minima. Engelbrecht and Wolmarans (2020)  Moloto and Engelbrecht (2020) showed that using the turbulent quantities observed during solar cycle 20 as inputs to an ab initio modulation code using the same set of diffusion and drift coefficients as employed by Moloto et al. (2018) yielded model GCR proton temporal intensity profiles in qualitative agreement with the neutron monitor observations corresponding to the same time period. This resolved a long-standing conundrum presented by the seemingly anomalous relationship observed between the heliospheric magnetic field magnitude at Earth (which stayed relatively uniform as a function of time) and GCR intensities as observed using neutron monitors (which displayed a temporal profile typical of A < 0 HMF conditions). Moloto and Engelbrecht (2020) and Engelbrecht and Moloto (2021) extended this model, with the former study employing simple temporal scalings for large and small scale heliospheric plasma quantities constructed from spacecraft observations and a time-dependent HCS to compute GCR proton intensities from the 1700s to the present, achieving good agreement with observations during the space age, and demonstrating the importance of taking into account charge-sign dependent modulation effects when studying historic cosmic ray modulation. This can be seen in Fig. 5, taken from that study, which illustrates the clear charge-sign dependence in galactic CR proton intensities computed using extrapolations of the temporal behaviour of the tilt angle from sunspot numbers following the approach of Asvestari and Usoskin (2016), and historic heliospheric magnetic field values at Earth inferred by McCracken and Beer (2015) from beryllium-10 observations (for more detail on very longterm modulation effects, see the review by Usoskin 2013). Note that the grey and cyan lines in Fig. 5 correspond to computed CR intensities, as opposed to the black line, which corresponds to spacecraft observations reported by Gieseler et al. (2017). Engelbrecht and Moloto (2021) directly employed spacecraft observations of turbulence as well as large scale plasma quantities as inputs to study GCR antiproton modulation, estimating the intensities of these particles at the heliospheric termination shock, demonstrating that these differed considerably from what was expected from galactic propagation models. Shen et al. (2021) extended the model of Shen et al. (2019) to consider the influence of observationally motivated, latitudinally varying turbulence quantities on GCR latitude gradients, finding good agreement with Ulysses observations during the first and third fast latitude scans, as illustrated in Fig. 6 (taken from that study), and concluding that the assumption of an anisotropic perpendicular diffusion coefficient is not necessary to fit the observations. That conclusion agrees with Moloto et al. (2019), who employed diffusion coefficients derived from different scattering theories, namely QLT and NLGC.
Several GCR modulation studies used outputs from turbulence transport models as direct inputs for diffusion and turbulence-reduced drift coefficients. Florinski and Zank (2006) Shen et al. (2021), shown alongside with Ulysses observations taken during the first and third fast latitude scans, as reported by Heber et al. (1996Heber et al. ( , 2008, de Simone et al. (2011), Gieseler and Heber (2016), Vos and Potgieter (2016) (figure adapted from Shen et al. 2021) demonstrated the sensitivity of computed GCR intensities on turbulence conditions and the properties of the surrounding interstellar space, including those when the solar system was immersed in the environment of the Local Bubble, representing low-density, hightemperature interstellar conditions (see ). This study employed, for the inner heliosphere, turbulence quantities modelled using the Zank et al. (1996) TTM, in conjunction with a full (axisymmetric) MHD model yielding large-scale background plasma quantities (see Florinski et al. 2003), and modeled the transport of turbulence across the termination shock by scaling magnetic variances up with the compression ratio. The diffusion coefficients employed by Florinski and Zank (2006) were modelled following Zank et al. (2004) with drift effects omitted, using a 2D finite-difference GCR transport code. In a subsequent study, Florinski and Pogorelov (2009) investigated GCR transport in the heliosheath during the 2008-2009 solar minimum using a 3D, stochastic solver for the Parker TE, demonstrating that the heliosheath presents a considerable modulation barrier for these particles, and reporting that GCRs reside for considerably longer periods in the heliosheath than they do inside the termination shock. That study employed, for background plasma quantities, the MHD model of Pogorelov et al. (2006), and modelled turbulence quantities in the supersonic solar wind by scaling the magnetic variance with the solar wind density. Across the termination shock, Florinski and Pogorelov (2009) scale the variance up with the square of the compression ratio, and beyond that assume that this quantity scales as the square of the HMF magnitude, employing these quantities as inputs for the parallel and perpendicular MFPs proposed by le Roux et al. (1999), Zank et al. (2004). Note that these authors did not consider the influence of compressible turbulence on GCR transport coefficients. Engelbrecht and Burger (2013a,b) studied the transport of GCR protons and electrons, as well as their antiparticles, inside the termination shock, using turbulence quantities computed using the two-component Oughton et al. (2011) TTM as inputs for QLT parallel MFPs, as well as a perpendicular MFP expression derived from the extended NLGC theory proposed by Shalchi (2006). Their computed intensities of GCR protons and antiprotons were in reasonable to good agreement with spacecraft observations at Earth, as can be seen from Fig. 7. The authors reported a reduction in GCR proton latitude gradients due to a latitude dependence in particle diffusion coefficients resulting only from the latitude dependence of turbulence quantities, as opposed to parametric increases in meridional per- pendicular diffusion coefficients, or the influence of additional meridional transport due to a Fisk (1996)-type HMF (their study employed a purely Parkerian HMF model). Engelbrecht and Burger (2015a) demonstrated, using the same TTM and perpendicular diffusion coefficients derived from the unified nonlinear theory of Shalchi (2010), the extreme sensitivity of computed GCR intensities to the choice of theoretical models for the turbulence-reduced drift coefficient, highlighting the importance of this quantity to numerical modulation studies. Engelbrecht and Burger (2015b) demonstrated the sensitivity of GCR proton intensities to the assumed low-wavenumber behaviour of the 2D turbulence power spectrum, which strongly influences the perpendicular diffusion coefficients calculated using the NLGC-type theories (see also . Guo and Florinski (2016) investigated the transport of GCRs in the presence of corotating interaction regions, using the single-component Breech et al. (2008) TTM in conjunction with a MHD code (see Guo and Florinski 2014, for more detail) as inputs to calculate a QLT derived parallel diffusion coefficient and a simplified perpendicular diffusion coefficient (see le Roux et al. 1999). These authors reproduced the recurrent variations associated with these structures, and attributed them to changes in diffusion coefficients associated with stream interfaces, as opposed to the action of the heliospheric current sheet. Furthermore, Engelbrecht (2019b) showed, using the GCR modulation model of Engelbrecht and Burger (2015b), that computed low-energy GCR electron intensities display a sensitivity to the assumed low-wavenumber, as well as the high-wavenumber, behaviour of the turbulence power spectrum, emphasizing the need to accurately model a broad range of turbulence quantities in particle transport studies.

Beyond the Parker Formulation
The Parker transport equation, valid for nearly isotropic distributions, is of limited use at low particle speeds (v ∼ u, where u is the flow speed of magnetic fluctuations) and in situations where scattering is weak and the distribution function is driven away from isotropy by transport effects. By definition, a strongly anisotropic distribution function has f (p) − f (p) ∼ f (p), where the subscript denotes averaging over a solid angle in momentum space. The simplest mechanism for producing anisotropy is adiabatic focusing in a diverging magnetic field, which is a common phenomenon in SEP transport (Bieber et al. 1986;Dröge et al. 2006). Contrary to the nearly-isotropic Parker transport equation that contains no explicit dependence on the mean magnetic field B 0 , a SEP propagation model must include, at the lowest order, the effects of changing magnetic field magnitude and the corresponding change in the pitch angle to conserve the magnetic moment of the particle (Roelof 1969). Some common locations with large ion anisotropies are the upstream regions of interplanetary shocks (Pesses et al. 1982), the Earth's bow shock (Möbius et al. 2001), and the solar wind termination shock (Decker et al. 2005), where field aligned ion beams are often observed owing to reflection off the shock ramp.
An early formulation of the evolution of a nearly gyrotropic distribution function, f (z, μ, t) was captured by Roelof (1969) in his well-known equation: for which the third term on the left describes anisotropic behavior via adiabatic focusing, comprised of cosine of pitch angle μ, particle speed v, pitch-angle diffusion coefficient D μμ , and focusing length, L, where: Equation (40) is the precursor to the equation of focused transport (FT) widely used today. Earl (1976) obtained a solution to Eq. (40) by performing an eigenfunction expansion on the first two terms, and showed that the transport changes from diffusive to weakly focusing, to coherent, to strongly focusing and super-coherent. This formulation enabled the effects of pitch angle scattering and adiabatic focusing to be combined, an attribute that enabled, for example, Ma Sung and Earl (1978) to successfully fit over 30 SEP intensitytime profiles. Ng and Wong (1979) derived numerical solutions to these early equations and adopted even more realistic approximations of the IMF. They found that monopolar and Archimedean-spiral mean fields behaved very differently from the exponential profile often used at the time, and that the monopolar solution agreed fairly well with predictions from classical diffusion. Ruffolo (1991) used the Ng and Wong (1979) approach to simulate the decay of protons from solar flare neutrons. They were able to successfully reproduced the flare features observed by ISEE 3. Ruffolo (1995) added the effects of deceleration and convection in addition to pitch-angle scattering and magnetic fluctuations. By simulating the transport of solar flare protons, they demonstrated the importance of including adiabatic deceleration in FT models, thereby laying a more accurate framework for the modeling the evolution of pitch angle distributions in the expanding solar wind. Focused transport studies have also included the effects of shock acceleration (e.g. Kirk and Schneider 1987;Ruffolo 1999), which has led to an explanation of precursory anisotropies observed before Forbush decreases (Leerungnavarat et al. 2003), and to model particle acceleration in compression regions (Giacalone et al. 2002). For a review of solar energetic particle focused transport studies, the interested reader is invited to consult Klein and Dalla (2017), .
An extended FT equation (FTE) (Skilling 1971(Skilling , 1975Isenberg 1997) can be derived under the assumption of gyrotropy of the distribution in the frame associated with the moving plasma, where the electric drift is absent, and the advection term acquires the plasma flow velocity. Note that the pitch angle is measured relative to the magnetic field at the actual particle's position, so drift acceleration is retained by the FTE. Of course, the distribution is actually gyrotropic not in the plasma frame, but in the guiding center frame, by its definition. So an alternative approach would be to transform to the GC frame, so that the drift advective terms would reappear in the FTE. However, it can be argued that the gradient and curvature drifts are of a higher order in r g /L (where L here denotes a characteristic length scale of the electromagnetic field) compared to the electric drift, so the resulting equation would be order-mismatched. For a more in-depth discussion of this matter, see Appendix B of le Roux et al. (2007). Given that the FTE is derived under the assumption of gyrotropy in the plasma frame, the momentum coordinate in the FTE refers to momentum in the moving, non-inertial plasma frame. By transforming the Boltzmann equation into the plasma frame and averaging over the gyrophase, the resulting pitch-angle dependent transport equation is readily shown to be ∂f ∂t Here b = B/B, and d/dt = ∂/∂t + (u · ∇) is the convective derivative. The term dependent on ∇ · b describes the effects of focusing and mirroring due to the change in magnetic field strength. The terms proportional to b · du/dt account for the flow inertia; they are second order in u/v and are often neglected. The remaining terms containing the derivatives of the plasma velocity describe the effects of flow compression (∼ ∇ · u) and shear (∼ ∇u) on the velocity and pitch angle of the particles. The right hand side describes scattering in pitch angle by magnetic irregularities. It should be understood that Eq. (42) is not a direct generalization of the Parker TE. The reason is that the gyrotropic approximation eliminates spatial curvature, polarization, and gradient drift effects, as well as diffusion normal to the mean magnetic field. It does, however, include the changes in parallel and perpendicular velocities from magnetic focusing, temporal changes of the magnetic field strength (the so-called betatron effect), and the work of the motional electric field on the guiding center drift motions ).
Because the FTE has one extra velocity coordinate, its use is typically restricted to systems with reduced spatial dimensionality. A popular way to numerically solve (42) is by the Monte-Carlo method, where multiple guiding center trajectories are integrated in time from the point of entry and until a spatial or momentum boundary is crossed Fahr 1998, 2000;Gieseler et al. 1999). The FTE can also be solved in the plasma frame where only the streaming, focusing, and scattering terms are present; this requires transforming the distribution between fixed and moving frames at each step during the simulation .
The FTE can be transformed by expanding f in Legendre polynomials, which results in an infinite system of hyperbolic PDEs . Truncating the expansion at order n yields n + 1 characteristics that correspond to coherent propagation along the magnetic field at some fraction of particle's speed. Even truncations are more useful because they capture the non-propagating mode (particles with pitch angle of 90 • ). The n = 2 system consists of three modes, namely the isotropic intensity, the field-aligned streaming mode responsible for first-order anisotropy, and the second order mode corresponding to distribution compression or elongation along B. Solving a three-equation system is much less costly numerically than solving the full FTE. A comparison between the two methods was performed by Florinski et al. (2008a) in the context of particle acceleration at quasi-perpendicular shocks. The resulting omnidirectional intensities near the shock were qualitatively similar, but the Legendre model yielded an unexpectedly hard power law spectra at high energies (> 1 MeV).
The shock acceleration problem treated with FTE is quite different from the nearly isotropic version because of the possibility of particle mirroring from a kink in the magnetic field. Louiville's theorem then implies that the intensity of particles within the loss cone is continuous across the shock, while the mirroring particles only contribute to the intensity upstream. As a consequence, particle intensity is discontinuous across the shock (Ruffolo 1999), something that is not permitted by the diffusive theory. Monte-Carlo simulations of particle acceleration by shocks (Gieseler et al. 1999) and compression regions (Klappong et al. 2001) featured a mirroring peak where accelerated particles are concentrated in a precursor ahead of the shock front. The height of the peak is inversely proportional to the particle velocity and the pitch-angle scattering rate. The exponential shock precursor is narrower in the anisotropic case compared with the diffusive limit. Spectra of accelerated particles at highly oblique shocks can develop a pronounced dip followed by a hump as a consequence of efficient acceleration by adiabatic reflection where particles receive a large boost in energy from reflection Florinski et al. 2008b). Dröge et al. (2006) used the FTE to model the spectra, temporal profiles, and anisotropies of protons and iron nuclei following impulsive SEP events. They emphasized the role of adiabatic acceleration in the nearly scatter-free regime to shift the injected spectra toward lower energies. Subsequently,  added cross-field diffusion to the FTE in an adhoc manner as a random walk process of the guiding center in the plane normal to the mean magnetic field. They found that more realistic results were obtained with the perpendicular MFP λ ⊥ proportional to the perpendicular component of the particle's velocity, which appears to be contrary to the weakly nonlinear theory prescription for incompressible turbulence (Shalchi et al. 2004c;, where λ ⊥ depends on the parallel velocity. The model was able to reproduce the intensity dropouts often observed at 1 AU using the ratios of λ ⊥ /λ as small as 10 −5 . A similar model was used for SEP electrons in the energy range measured by the Wind spacecraft (Dröge et al. 2018) where it was found that QLT prescriptions for resonant pitch-angle scattering based on observed spectra of slab fluctuations were inconsistent with observed intensities and angular distributions at 1 AU. Low energy electrons appeared to experience nearly isotropic scattering that would require non-resonant effects such as transit time damping in the presence of compressive turbulence. The model was applied for particles with v u, so all plasma advection-related terms in Eq. (42) were ignored.
The full FTE was used by Chalov and Fahr (2003) to simulate the suprathermal tail often observed in the spectra of pickup protons. They argued that the tail could be explained by secondary charge exchange of solar-wind protons with energetic neutrals originating from the inner heliosheath. A time dependent model of pickup ion acceleration at the termination shock based on a solution to the FTE was presented by le Roux and Webb (2009). Unlike most studies that used a stochastic trajectory approach, this model solved the FTE on a mesh in a three-dimensional phase space (r, p, μ) that restricted the model to a single spatial dimension. In focused transport theory injection of particles into the shock acceleration process requires that v > u 1 cos θ bn , i.e., the particle speed must exceed the de Hoffmann-Teller speed, which is much larger than the shock speed for high obliquities. However, the authors calculated the probability distribution of magnetic field angles from Voyager observations and found that, while the termination shock was nearly perpendicular on average, for about 10% of the time the shock obliquity was below 60 • . The model varied the magnetic field in the upstream region in a stochastic fashion that both reduced the injection threshold and introduced effective cross field diffusion due to field line meandering, allowing ions from energies as low as 1 keV to be efficiently accelerated. A similar approach was used by Chalov et al. (2016), who implemented a local anisotropic model of pickup ion transport in the vicinity of the termination shock embedded in a global model describing the evolution of plasma and neutral atoms over the entire heliosphere. Because of the stochastic approach their model was fully threedimensional. The model included small scale variations of the magnetic field angle with the shock normal, as well as large scale variations from the magnetic sector structure. It was demonstrated that ion spectra measured by Voyager 1 and Voyager 2 up to a few MeV could be explained by shock drift acceleration where the shock normal angle was a random parameter selected based on the probability distribution deduced from Voyager magnetic field observations (Fig. 8). Zuo et al. (2011) presented a detailed description of solving the FTE in the vicinity of a shock with a forward-in-time stochastic approach. They found that the FTE produced the same power law spectrum near quasi-parallel shocks as predicted by diffusive shock acceleration theory. Conversely, at more oblique shocks the spectrum had an initially harder segment followed by the standard DSA spectrum at the highest energies. Voyager 1 observed highly anisotropic particle beams during the two years preceding its termination shock crossing (Krimigis et al. 2003). Kóta and Jokipii (2004) modeled acceleration and transport of low-energy "termination shock" particles as well as higher energy ACRs using the FTE in a form somewhat different from Eq. (42). Because in the absence of perpendicular transport the particles are tied to the magnetic field lines, one could solve the FTE by integrating the trajectories of the guiding centers along the field lines. The latter were assumed to cross the termination shock multiple times. The model produced a two-population spectrum where low energy ions were recently injected and did not have  Chalov et al. (2016) sufficient time to be accelerated to ACR energies. The calculated pitch-angle distribution exhibited beaming in the sunward direction, in agreement with observations. The same model was used by Kóta and Jokipii (2017) to explain the anisotropy events observed by Voyager 1 in the LISM that were thought to be related to a passage of compressive structures past the spacecraft (Gurnett et al. 2015). It was shown that GCR depletion near the 90 • pitch angle could be produced by energy loss of trapped ions in the expansion region following the compression. While the theory of Kóta and Jokipii (2017) explains several key aspects of the observations, many issues remain unresolved, particularly related to the episodes of the GCR intensity depletions near 90 • pitch angle (Rankin et al. 2019b). The reason for the longevity of the events, as well as their recovery mechanism remains unclear. Further, their timing is not always well correlated with the behavior of the local field, leading several authors to suggest that the large-scale draping of the field may influence the GCR trapping (Rankin et al. 2019b;Hill et al. 2020). Moreover, while the events are evident for protons, the electrons show no clear evidence of anisotropy, despite both species being clearly reflected and accelerated at the forefront of the shock (Rankin et al. 2020). This appears contrary to the measurements of the ambient turbulence spectra (Burlaga et al. 2018), which is higher in amplitude for wavelengths that interact with the protons as opposed to the electrons. Fraternale et al. (2020) suggested pickup ion instability as a candidate process for generating high-frequency turbulence capable of isotropizing the electrons, a topic that awaits further investigation. See the Rankin et al. chapter of this book for further details.
In recent years, observations of CR anisotropy have garnered much attention. Here we give a brief review of CR anisotropy simulations as they relate to observations in the solar wind and heliosphere, including those at TeV energies, near stream-interfaces, and encountered via in-situ measurements at the TS and HP. The large-scale compressive structures responsible for the transients seen by the Voyagers in the very local interstellar medium (VLISM) result from the pile-up of solar events that coalesce and form merged interaction regions capable of propagating through the heliosphere and transmitting pressure waves across the boundary. In mid-2012, such an event was observed to modulate cosmic rays first in the heliosheath at Voyager 2 and later in the interstellar medium at Voyager 1, as reported by Rankin et al. (2019a) (see the Rankin et al. 2022 paper for further details). The time profiles of the GCR intensity responses were remarkably similar (cross correlation of 91.2 ± 2.1%), which posed an intriguing question: why did the modulation appear isotropic in the heliosheath but highly anisotropic in the VLISM? This was investigated by , who employed a Vlasov-Fokker-Planck equation to simulate particle transport in a shock propagating through a simple paramagnetic shielding model of the heliosheath. They demonstrated that trapping and adiabatic cooling of particles in the rarefied fields downstream of the disturbances occurred in both places. However, the weak scattering of the interstellar medium enables pitch angles to be preserved, while the strongly-turbulent field in the heliosheath more effectively isotropizes the distribution.
GCRs at TeV to PeV energies have been observed to exhibit a small but significant sidereal diurnal anisotropy as they propagate through the heliosphere and arrive at Earth, with amplitudes of ∼ 10 −3 on large scales (> 60 • ) and ∼ 10 −4 on small scales (10 • to 30 • ) (see, e.g., Nagashima et al. 1998;Hall et al. 1999;Guillian et al. 2007;Munakata et al. 2010;de Jong 2011;Abbasi et al. 2012;Aartsen et al. 2016;Amenomori et al. 2017;Bartoli et al. 2018;Abeysekara et al. 2019). The origins of this anisotropy are actively debated, and various mechanisms have been explored by models such as the altering of particle transport by magnetic turbulence in the local interstellar medium (e.g. Giacinti Schwadron et al. (2014) showed that the direction of the global anisotropy at TeV energies is consistent with the mean field derived from the IBEX ribbon, indicating that the GCR propagation paths are mostly along the field, with very little cross-field diffusion. Analytical models by several authors gained some traction, showing that the small-scale anisotropies can arise from the interstellar magnetic field's configuration (e.g. Giacinti and Sigl 2012;Battaner et al. 2015;López-Barquero et al. 2016). Nonetheless, the origins of the global anisotropyboth its amplitude and shape -remain the subject of debate. A simple dipole anisotropy could potentially arrive from the heliosphere drifting through an isotropic distribution of GCRs caused by the Compton-Getting effect, but observations indicate that the anisotropy is far more complex (see, e.g., Aartsen et al. 2016).  account for these complexities through the method of Liouville Mapping -by phenomenologically representing the large-scale anisotropy as a sum of Legendre polynomials (see also , and references therein), a topic that is discussed in greater detail in the Rankin et al. (2022) paper.

Beyond the Diffusion Approximation
In addition to the phenomena of magnetic focusing and non-isotropic pitch-angle distributions, a second aspect of particle transport has received increased interest in recent years, namely the anomalous diffusion of energetic particles. The basic idea is that nonhomogeneous structures in the background plasma can lead to deviations in the classical (Gaussian) diffusion approach, which assumes a certain level of homogeneity of the background medium. On a macroscopic level, this can manifest itself as deviations of the temporal behavior of the particle's average moments. Most prominently, the linear time dependence of the second moment can be different. If x is the particle position and the mean square displacement dependence on time is described by the transport is called superdiffusive for α > 1 and subdiffusive for α < 1. The case where α = 1 corresponds to the classical Gaussian diffusion case where each step in the random walk is uncorrelated with the preceding steps. Anticorrelation ("negative memory") leads to subdiffusion and correlation ("positive memory") leads to superdiffusion. Note, however, that this is only a very coarse, high-level characterisation. The behaviour in other particle moments might be critical. Also, a combination of sub-and superdiffusive processes can lead to a macroscopic α = 1 behaviour in special cases, while the underlying processes are not necessarily Gaussian (e.g. Magdziarz and Weron 2007). A good overview of the application of anomalous diffusion processes in astrophysical plasmas (with connections to laboratory plasmas) is given in . The observational evidence for anomalous diffusion in the heliosphere and its interpretation are already discussed in some length in a chapter by Perri et al. in this journal. Here, we will focus on specific aspects of the theoretical description relevant to particle transport in the heliosphere.

Non-Gaussian Transport Processes Induced by Field-Line Turbulence
Some evidence for non-Gaussian transport of field lines and particles in the solar wind comes from observations of "dropouts" in which SEP fluxes undergo repeated, non-dispersive increases and decreases . These have been interpreted in terms of filamentary, non-diffusive connection via magnetic field lines back to a narrow source region at the Sun . The filamentation and dropout pattern have been explained in terms of temporary topological trapping Zimbardo et al. 2004;Tooprakai et al. 2016) by the flux-tube structure of the solar wind (McCracken and Ness 1966;Burlaga 1969), which can be modeled by a 2D component in the 2D+slab model of magnetic turbulence. It should be emphasized that such behavior arises even for mathematical models of homogeneous 2D+slab turbulence, and it can arise at any location in such models; in other words, no boundary condition is used to impose special structures such as magnetic flux ropes. Further work identified local trapping boundaries (even in homogeneous turbulence) across which the transport of magnetic field lines is suppressed Seripienlert et al. 2010).
Other observations of sudden changes in the energetic particle flux in interplanetary space could be explained by local trapping boundaries (Tessein et al. 2016;. Such non-Gaussian transport has also been proposed to explain extreme ultraviolet "moss" emission from the solar transition region, with a network pattern at the base of hot coronal loops (Kittinaradorn et al. 2009). Drift motions due to the 2D turbulence structure lead to local squeezing of the solar particle distribution. This process was invoked to explain eventto-event variability in the intensity of relativistic solar particles during GLEs .
Another example of non-Gaussian transport in a turbulent plasma is "compound" transport (Getmantsev 1963;Urch 1977;, which we call compound subdiffusion. This can result from the random walk of particles along field lines ( z 2 ∝ t ) that themselves undergo a random walk in space ( x 2 ∼ z), so that x 2 ∝ t 1/2 . Referring to Eq. (44), this corresponds to α = 1/2, a subdiffusive case. Such compound subdiffusion mainly applies to slab turbulence in which particles are tied to a single field line (see Sect. 5.1.1). Physically, subdiffusion results from the "negative memory" between steps in the random walk, as backtracking along field lines implies that later steps in the plane perpendicular to the large-scale magnetic field can be anticorrelated with the earlier steps. Subdiffusion in slab-like turbulence was demonstrated by numerical simulations of Qin et al. (2002b). However, once the magnetic field acquires a degree of transverse complexity (as in the 2D+slab model) particles undergo transient subdiffusion at early times followed by asymptotic diffusive behavior at later times Qin et al. (2002a). This is due to true cross-field diffusion, in which the particle no longer follows the same field line; consequently the random walk loses the memory of previous steps and becomes diffusive.
One should also note that an arguably trivial case of superdiffusive transport is ballistic transport, in which particles have not had enough time to change direction so that x ∝ t for individual particles and x 2 ∝ t 2 , corresponding to α = 2 in Eq. (44). In general, diffusive behavior (representing a random walk with no memory) is always preceded by such an initial superdiffusive phase while t τ , where τ is a mean free time, as seen in Fig. 3.

Modelling of Anomalous Transport
The generalization of transport equations to anomalous diffusion behavior leads to the notion of fractional Fokker-Planck equations (FFPEs) in space and time (see, e.g., Metzler and Klafter 2000). For the treatment of superdiffusion, the second-order Laplacian for the spatial diffusion is usually replaced by a diffusion operator of fractional order (for an alternative see, e.g., Sokolov and Metzler 2004). An example derivation of the fractional diffusion equation starting from the equation of continuity and generalizing Fick's law can be found in Chaves (1998). For a more detailed physics-based derivation of the resulting integral equations, see also Chukbar (1995). A proper generalization of the Laplacian is the so-called Riesz fractional derivative, which can be expressed by the Riemann-Liouville fractional derivatives, here given in a form with integration domain to the left of x: with m − 1 < μ < m, where m is an integer (note that in this section μ, following standard notation, denotes the order of the fractional derivative and not the pitch-angle cosine). The symmetric fractional Riesz derivative ∂ μ /∂|x| μ (Gorenflo et al. 1999;Metzler and Klafter 2000) is then given by and generalizes the Laplace operator. A corresponding (reduced) transport equation, or fractional Parker equation (Litvinenko and Effenberger 2014) using this definition, can be written as where k μ is a generalized diffusion constant with the appropriate physical dimensions . This type of equation can be solved numerically with grid-based discretization schemes (see, e.g., Stern et al. 2014, and references therein) and with stochastic methods, Fig. 9 Pseudo-particle paths of a two-dimensional Lévy process with μ = 1.5. The trajectories show alternating episodes of unperturbed propagation and intense local scattering. Note that the trajectories do not describe real particle motion but rather phase-space elements of the underlying distribution function similar to the Itô-equivalent SDEs. Analytical solutions that asymptotically approximate the full solution of such a diffusion-advection equation are discussed in Litvinenko and Effenberger (2014). The existence of moments of the solution to the equation depends on the parameter μ. In order for higher order moments to strictly exist, it may be necessary to move to tempered distributions or e.g. Lévy Walk processes (see, e.g., the discussion in Metzler and Klafter 2000). Particle conservation is, however, fulfilled, since the zero-order moment of the distribution is finite. The derivation of a fractional diffusion equation via a generalized Fick's law (Chaves 1998) further illustrates this notion. When solving the generalized problem of anomalous diffusion using the SDE approach the classical Wiener process (i.e. a Gaussian stochastic process) has to be replaced by heavytailed Lévy distributions (Jespersen et al. 1999;Magdziarz and Weron 2007), which results in distinctive patterns for the behaviour of phase-space elements with extended periods of scatter-free Lévy flights alternating with periods with local scattering. The SDE takes the form: where dL μ (t) is now a μ-stable Lévy motion with Fourier transform e ikLμ(t) = e −t|k| μ , wherek is a wavenumber. Figure 9 shows sample pseudo-particle trajectories of such a process in two dimensions, illustrating the numerical integration using appropriate random number generators for these generalized probability distributions (Chambers et al. 1976). Note that the recent work of le Roux and Zank (2021) derives such fractional kinetic equations from first principles for a scenario of particle trapping in small scale magnetic flux ropes. Together with evidence from solar flares that points to anomalous acceleration behaviour in fragmented current sheets (Bian and Browning 2008;Isliker et al. 2017), these approaches provide a strong argument to study anomalous transport and acceleration processes in the heliosphere in much more detail in the future.

Summary and Outlook
As discussed above, extraordinary advances, made in the past few decades, in our understanding of heliospheric turbulence, its transport, and the plasma physics of the greater heliosphere itself, coupled with great strides in our understanding of the scattering of energetic charged particles traversing turbulent magnetic fields, enable cosmic ray transport modelling at an unprecedented level of complexity and detail. Such models, in turn, provide valuable insights into the transport and modulation of cosmic rays. Relatively recently, these advances in CR transport theory experienced broader application in the studies of SEP transport. For example Tooprakai et al. (2016) studied the transport of particles originating in impulsive solar flares in the presence of composite turbulence, and Laitinen et al. (2016) investigated the longitudinal transport of SEPs due to field line meandering (see also studies by Chhiber et al. 2021b,a). Furthermore, studies that employ QLT and FLRW diffusion coefficients, such as , emphasised the importance of perpendicular diffusion in the transport of SEPs, while van den Berg et al. (2021) investigated the influence of turbulence on potential drift effects experienced by such particles. Given recent interest in particle transport in various asterospheres (as reviewed elsewhere by Herbst et al. 2022 in this journal), the study of CR transport in these novel scenarios may provide useful insights into, for example, the potential habitability of the exoplanets they contain (see, e.g., Mesquita et al. 2021). Furthermore, these advances can be of great benefit to particle transport studies in more exotic astrophysical contexts, such as pulsar wind nebulae (see, e.g., Vorster and Moraal 2013;Porth et al. 2016).
Many challenges still remain. The next generation of models will need to move beyond qualitative comparisons with spacecraft observations, and undertake careful comparisons with precision measurements of cosmic ray intensities such as those reported by PAMELA (for two such comparisons that have already been made, see Moloto et al. 2018;Shen and Qin 2018) and AMS. These new models will be taking into account our increasing knowledge of the complexities of the heliospheric plasma environment (e.g. ) and the corresponding theoretical developments (e.g. ). This will require further advances in turbulence transport modelling, and further improvements in our understanding of cosmic ray transport coefficients, given the differences in values for these quantities yielded by various theories discussed in Sect. 5.2. For example, the assumption of transverse turbulence is explicitly made in the various scattering theories discussed in this review. This limits their applicability in the heliosheath, where observations show that a significant compressible component is present in the HMF turbulence, as reviewed by Fraternale et al. elsewhere in this journal. This is a significant issue, considering that Voyager observations show that most of the modulation experienced by galactic cosmic rays occurs in this region. To date, relatively few theoretical studies have investigated particle transport coefficients in the presence of such turbulence (see, e.g., Schlickeiser 2002;Lazarian 2002, 2004;le Roux et al. 2005;Lazarian and Beresnyak 2006;Strauss et al. 2016), yielding relatively simple results for these quantities in specific, usually astrophysical, contexts. It remains to be seen what sort of influence their use in an ab initio CR transport model would have on GCR intensities and ACR modulation.
As MHD models of the entire heliosphere, including the tail region, become more refined (see, e.g., Opher et al. 2021, and references therein), future work will combine such models with state-of-the-art GCR transport models. One could hope, for example, that such modelling efforts would be able to reproduce the sharp transition of cosmic rays observed by the two Voyager spacecraft when transiting the heliopause boundary (within fractions of an astronomical unit), which in itself indicates a stark contrast between diffusion coefficients in the heliosphere compared to the very local interstellar medium. Several models have demonstrated that small diffusion coefficients could cause a gradual intensity change that extended tens to hundreds of au beyond the heliopause (e.g. Scherer et al. 2011;Strauss and Potgieter 2014). Several models show that a dramatic increase in the ratio of the parallel to perpendicular diffusion coefficients is required to explain the observations (e.g. Guo and Florinski 2014;, although the fine spatial resolution required and the need to include accurate models of both fields and particles at the two boundaries remain open challenges (see the Rankin et al. chapter for more details). Further theoretical insights would greatly enhance our understanding of the fundamental physics behind the transport of these particles in this region. Recent observations have also presented new challenges to modeling ACR and GCR transport near the Sun. The Parker Solar Probe reported much stronger ACR radial gradients than expected by models, revealing a significant deviation from the typical r −2 Parker-spiral field behavior near the Sun (see Rankin et al. 2021a,b, andthe Giacalone et al. 2022 paper for more detail). Many questions arise about what sort of assumptions should be made when modeling particle transport in this new regime. For example, the unknown behavior of drifts close to the Sun makes the treatment of the inner boundary condition uncertain. In addition, the enhanced turbulence levels near the Sun may considerably reduce the drift effects (see van den Berg et al. 2021). Moreover, large gradients in the magnetic field could lead to strong magnetic focusing or mirroring, rendering the traditional Parker transport models invalid.
Lastly, the connections between heliospheric and cosmic ray transport in the broader galaxy are yet to be fully explored. Such an exploration would be of great use to the study of the galactic propagation of cosmic rays (see, e.g., Moskalenko et al. 2003), specifically in terms of their transport parameters in the interstellar medium (see, e.g., Yan and Lazarian 2008), which in turn would inform heliospheric modulation studies with improved estimates of the very local interstellar spectra of these particles. Such information would also be of great use in the search for the signatures of dark matter annihilation (see, e.g., Cholis et al. 2019;von Doetinchem et al. 2020). Studies of such connections would benefit from any potential observational insights yielded by the planned interstellar probe discussed by Brandt et al. (2022) in this journal. Acknowledgements This work was made possible by the International Space Science Institute and its interdisciplinary workshop "The Heliosphere in the Local Interstellar Medium", www.issibern.ch/workshops/ heliosphere.

Funding Note
Open Access funding enabled and organized by Projekt DEAL. This work is based on the research supported in part by the National Research Foundation of South Africa (Grant Numbers: 111731, 119434, 137793), Thailand Science Research and Innovation (Grant Number: RTA6280002), and United States National Aeronautics and Space Administration (Grant Numbers: 80NSSC18K1209, 80NSSC20K0786). Opinions expressed and conclusions arrived at are those of the authors and are not necessarily to be attributed to the NRF or TSRI.

Data Availability Not applicable.
Code Availability Not applicable.

Conflicts of interest/Competing interests The authors declare no conflict of interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.