A Focused Transport-based Kinetic Fractional Diffusion-advection Equation for Energetic Particle Trapping and Reconnection-related Acceleration by Small-scale Magnetic Flux Ropes in the Solar Wind

Analysis of energetic particle inner heliospheric spacecraft data increasingly suggests the existence of anomalous diffusion phenomena that should be addressed to achieve a better understanding of energetic particle transport and acceleration in the expanding solar wind medium. Related to this is fast-growing observational evidence supporting the long-standing prediction from magnetohydrodynamic ( MHD ) theory and simulations of the presence of an inner heliospheric, dominant quasi-two-dimensional MHD turbulence component that contains coherent contracting and merging ( reconnecting ) small-scale magnetic ﬂ ux rope ( SMFR ) structures. This suggests that energetic particle trapping in SMFRs should play a role in anomalous diffusion in the solar wind that warrants further investigation. However, progress in studying such anomalous energetic particle transport phenomena in the solar wind is hampered by the lack of a fundamental derivation of a general fractional kinetic transport equation linking macroscopic energetic particle fractional transport to the microscopic physics of energetic particle interaction with SMFR structures. Here, we outline details of how one can derive a closed ensemble-averaged focused transport equation in the form of a general kinetic fractional diffusion-advection equation from ﬁ rst principles following the nonlinear Eulerian correlation function closure approach of Sanchez et al. With this equation one can model the anomalous diffusion of energetic particles in ordinary, momentum, and pitch-angle space in response to particle trapping in numerous SMFRs advected with the solar wind ﬂ ow.


Introduction
The number of studies finding potential anomalous transport features in energetic particle spacecraft data in the solar wind have been steadily increasing. For example, Mazur et al. (2000) and Chollet & Giacalone (2008) reported large abrupt drops in the counting rate of solar flare generated solar energetic ions (SEPs) at 1 au over small spatial scales of the order of the solar wind turbulence correlation length, suggesting that these sudden drops are connected to the properties of low-frequency magnetohydrodynamic (MHD) turbulence. Ruffolo et al. (2003) offered an explanation for these "dropout events" that involves SEPs from solar flare events populating two classes of magnetic flux tubes. SEPs that propagate to 1 au while being trapped in closed small-scale quasi-helical magnetic field structures or small-scale magnetic flux ropes (SMFRs) will retain a relatively high counting rate, while those SEPs propagating along adjacent open magnetic field lines that spread efficiently in latitude and longitude will arrive at 1 au with a strongly reduced counting rate. The trapped SEPs can be thought of as undergoing subdiffusive transport across the magnetic field.
The idea of particle trapping in SMFRs in the solar wind being responsible for observed potentially anomalous energetic particle transport behavior has received a considerable boost from the identification of unprecedented numbers (a few hundred per month) of SMFRs in the inner heliosphere from ∼0.3-7 au using Helios, Advanced Composition Explorer (ACE), Wind, Ulysses, and Voyager data (e.g., Hu et al. 2018;Chen et al. 2019;Zhao et al. 2019a;Chen & Hu 2020). SMFRs are defined as quasi-helical, coherent magnetic field structures advected with the solar wind flow. They are composed of a guide field component in the axial direction and a 2D twist or magnetic island component perpendicular to the guide field component (Cartwright & Moldwin 2010). SMFRs near Earth typically have cross sections -<0.1 0.01 au, which overlap largely with the inertial range scales of solar wind magnetic turbulence. These structures are different from pseudo-flux ropes such as torsional Alfvén waves propagating at the Alfvén speed ). There is a long-standing debate whether SMFRs mainly originate near the Sun (e.g., Borovsky 2008) or whether they are generated mainly locally in the solar wind. Evidence is increasing that in a space plasma with a significant guide/background magnetic field like the solar wind, SMFRs naturally form as part of a self-generated quasi-two-dimensional (quasi-2D) MHD turbulence component that might dominate other MHD wave turbulence modes. This is supported by observations in the slow solar wind near 1 au (e.g., Matthaeus et al. 1990;Bieber et al. 1996;Greco et al. 2009;, MHD simulations (e.g., Shebalin et al. 1983;Dmitruk et al. 2004) and nearly incompressible MHD turbulence transport theory (e.g., Zank & Matthaeus 1992Zank et al. 2017Zank et al. , 2018Zank et al. , 2020. Furthermore, observational evidence is accumulating that SMFRs also form locally in the solar wind near the heliospheric current sheet, large-scale current sheets associated with interplanetary coronal mass ejections, and corotating interaction regions through turbulent magnetic reconnection, and that neighboring SMFRs undergo merging through small-scale magnetic reconnection and large-scale compression (Khabarova et al. , 2016Hu et al. 2018;; Chen et al. 2019; Chen & Hu 2020).
Increasingly energetic particle acceleration events are observed behind interplanetary traveling shocks, and downstream of the heliospheric termination shock as well, that coincide with and are interpreted to be caused by contracting and merging (reconnecting) SMFRs. Typically, such events display an intensity amplification factor that increases with particle energy (e.g., Zank et al. 2015;Zhao et al. 2018;Adhikari et al. 2019;Malandraki et al. 2019;Zhao et al. 2019b), an effect that can be modeled with kinetic transport equations for energetic particle interaction with SMFRs that we developed (Zank et al. 2014;le Roux et al. 2015le Roux et al. , 2018. In these transport equations particle acceleration is modeled in terms of normal advection and diffusion in momentum space. However, two acceleration events were found where the normal intensity amplification scaling with energy holds at higher energies, but at lower suprathermal energies an inverse scaling of the amplification factor with energy applies, i.e., the amplification factor decreases with particle energy Khabarova et al. 2016). In Voyager 2 data from behind the heliospheric termination shock, Zank et al. (2015) found that the anomalous decrease in the intensity amplification factor for ions with particle energy at lower energies ranged from ∼0.028-3.5 MeV and the normal increase in the amplification factor at higher energies ranged from ∼1. 8-22.3 MeV. The same trend appeared in SEP data at 1 au from Solar Terrestrial Relations Observatory A and B and in ACE behind the traveling shock event of 2006 December 14 (Khabarova et al. 2016) where the intensity amplification factor for ions was observed to decrease with particle energy between ∼0.047 and 1.05 MeV nuc −1 and increase with particle energy between ∼13.6 and 100 MeV nuc −1 . These authors suggested that this anomaly in the particle acceleration occurs below a certain threshold energy where energetic particles have gyroradii much smaller than the SMFR cross-section scale so that they experience trapping effects. One would expect also that these particles experience subdiffusive perpendicular transport and it might be that trapped particles experience anomalous diffusive acceleration as argued below from the viewpoint of simulations.
In addition, in a series of papers Zimbardo and coworkers reported spacecraft observations of SEP intensity-time profiles ahead of interplanetary shocks that are power laws instead of the exponential profiles predicted by standard diffusive shock acceleration. They interpreted this as a sign of superdiffusive SEP transport that modifies the efficiency of the shock acceleration of SEPs (Perri & Zimbardo 2007Zimbardo et al. 2017.  detected superdiffusive transport upstream of three traveling shocks in ACE data with shock obliquities that are quasi-parallel, oblique, and quasi-perpendicular. In all these events, upstream wave amplification was negligible and thus not responsible for the power-law behavior. The results showed that superdiffusive transport was most pronounced for parallel transport upstream, while perpendicular transport upstream was less superdiffusive with diffusive behavior in the lower energy channels. Recently,  explained such superdiffusive transport as a product of subdiffusion in pitch-angle space arising from the observation that pitch-angle scattering times in the solar wind have a power-law distribution instead of a single value.
Besides observational evidence, there is ample evidence from simulations that anomalous diffusion occurs in space and astrophysical plasmas. Anomalous diffusion that can be either superdiffusive or subdiffusive has been detected in many simulations of energetic charged particles in prescribed magnetic turbulence that depends on parameters such as turbulence levels, geometry, kurtosis, and the ratio of the particle gyroradius over the turbulence correlation scale (e.g., Qin et al. 2002;Zimbardo et al. 2006;Pommois et al. 2007;Shalchi & Kourakis 2007;Tautz 2010;Beresnyak 2013;Alouani-Bibi 2014;Servidio et al. 2016). Furthermore, energetic particle trapping in magnetic islands produced by primary current sheet magnetic reconnection has been reported in many simulations (e.g., Drake et al. 2006;Guidoni et al. 2016;Xia & Zharkova 2018). In addition, the simulations of Drake et al. (2006Drake et al. ( , 2010, e.g., exhibited systematic acceleration of particles trapped in a contracting magnetic island, which may be a sign of superdiffusive energization. Superdiffusive energization is supported by the results of tracing test particle trajectories in the fragmented reconnection electric fields at reconnection sites and in electric fields produced by dynamic SMFRs as a consequence of a large-scale current sheets undergoing turbulent magnetic reconnection as expected to occur during jet and eruption events at the Sun (e.g., Isliker et al. 2017Isliker et al. , 2019. In many research fields it has been fruitful to model anomalous transport phenomena with the standard fractional diffusion equation. For an overview, see, e.g., Metzler & Klafter (2000. This equation is derived by taking the late-time and large-spatial-scale (fluid) limit of the Montroll-Weiss equation (Montroll & Weiss 1965). The Montroll-Weiss equation represents a model for microscopic continuous-time random walk (CTRW) that is in turn derived from a so-called master equation. The CTRW model generalizes the Brownian random walk model typically in terms of two independent general probability distribution functions for the waiting time between spatial steps and the step size that are specified to be Lévy stable distributions that exhibit powerlaw tails. The general anomalous diffusion equation in one dimension for the charged particle number density in a plasma is given by , where A is a constant, σ is a measure of the particle displacement, and τ is a reference waiting time. The implication of Equation (1) is that the average particle displacement is determined by | | ( ) á D ñ µ D a b x t . The nature of the fractional transport can be classified in terms of the Hurst fractional index b a = H , where subdiffusive transport holds when < < H 0 1 2 and superdiffusive transport occurs for < < H 1 2 1. In the special case where a = 2, subdiffusion requires b < 1 and superdiffusion is valid for b < < 1 2. However, if a < 2 one can have superdiffusion even when b < 1, because subdiffusion requires a b > 2 , while superdiffusion demands that a b < 2 . According to the literature survey by Sanchez et al. (2006) the fractional index constraints a < 2, b < 1 of Equation (1) combined with the superdiffusive requirement a b < 2 appear to be very much applicable to plasma transport in fusion devices.
It should also be noted that in certain problems, such as cosmic-ray reacceleration in astrophysics or solar flare associated acceleration of SEPs in the solar corona, researchers saw the need for modeling the acceleration as an anomalous diffusion process in momentum (velocity or energy) based on microscopic CTRW in momentum space (Bian & Browning 2008;Uchaikin 2010;Isliker et al. 2017Isliker et al. , 2019. For this purpose, Equation (1) was proposed with a fractional diffusion term in momentum space instead of ordinary space. However, in energetic particle transport problems in the solar wind advective transport can be important. For example, the diffusive shock acceleration of SEPs cannot be modeled without solar wind advection, while advection of energetic particles in momentum space is considered to be important in the transport theory of systematic particle acceleration by SMFRs in the solar wind (e.g., Zank et al. 2014;Zhao et al. 2018;le Roux et al. 2015le Roux et al. , 2019. Thus, a more appropriate model for fractional energetic particle transport in the solar wind is the so-called fractional diffusion-advection equation, which follows from extending the CTRW-based general anomalous diffusion equation (Equation (1)) to include spatial advection. Such fractional diffusion-advection equations are discussed in both the subdiffusive and superdiffusive limits (Metzler & Klafter 2000. Combining the two limits yields a general fractional diffusion-advection equation based on the CTRW model valid for b < 1 and a < 2 given by is the left Riemann-Liouville fractional time derivative. In the case of energetic particle transport in the solar wind ahead of interplanetary shocks, observed and interpreted by Perri & Zimbardo (2007 to be superdiffusive, a fractional generalization of the Parker transport equation (e.g., Litvinenko & Effenberger 2014;Zimbardo et al. 2017) has been put forward to model these observations. This equation in one spatial dimension is where ( ) f x p t , , is the energetic particle distribution function that also depends on the particle momentum p, and U is the solar wind speed. In this equation, a < < 1 2 to ensure superdiffusive transport ( < < H 1 2 1). This equation can be viewed as the superdiffusive limit of the general fractional diffusive-advection equation given by Equation (2). Such superdiffusive transport, viewed by  as largely a parallel transport phenomenon, is interpreted by  as a consequence of subdiffusion in pitch-angle space arising from the observation that pitch-angle scattering times in the solar wind have a power-law distribution that emphasizes long scattering times (long waiting times). To model such pitch-angle transport,  suggested the fractional pitch-angle transport equation where μ is the cosine of the particle pitch angle.

Motivation
As discussed above, energetic particle trapping in SMFRs in the solar wind is expected to play a key role in energetic particle perpendicular transport and acceleration being anomalous diffusion processes in ordinary and momentum space, respectively. However, progress in studying such anomalous energetic particle transport phenomena in the solar wind is hampered by the lack of a fundamental derivation of a general fractional kinetic transport equation linking macroscopic energetic particle fractional transport to the microscopic physics of energetic particle interaction with SMFR structures.
In previous publications we presented derivations of closed ensemble-averaged focused transport equations in the form of multidimensional kinetic diffusion-advection or Fokker-Planck equations with normal diffusion terms in pitch-angle and momentum space. These equations were derived by fundamentally linking macroscopic energetic particle transport to the microscopic physics of particle interaction with SMFRs (Zank et al. 2014;le Roux et al. 2015le Roux et al. , 2018. However, in these efforts no attempt was made to extend these derivations to include fractional transport terms for modeling particle trapping effects. Our normal diffusion-advection equations were derived in two limits. In the limit of particle interaction with SMFRs with a sufficiently weak twist or magnetic island component compared to a locally uniform guide magnetic field component, the standard quasi-linear theory (QLT) method was used where particle trajectories were assumed to be approximately undisturbed by SMFRs on short timescales in the derivation. In the second limit a stronger magnetic island component was assumed for which a nonlinear extension of the QLT method was employed by modeling particle trajectories to be disturbed on short timescales by SMFRs. This approach can be viewed as a simplified way of achieving the end result of a complicated renormalization treatment of the nonlinear terms (e.g., Owens 1974;Bykov & Toptygin 1993). The QLT limit yielded a normal diffusion-advection equation after closure because the disturbed particle trajectory Lagrangian correlation functions, arising from the microscopic physics of particle interaction with SMFR structures, were explicitly assumed to undergo statistically a rapid exponential decay (Poisson statistics), while the correlation functions were implicitly specified statistically to decay locally spatially. The nonlinear limit differs noticeably from the QLT limit because in the former case the disturbed particle trajectory correlation functions are Eulerian because the disturbed particle trajectory information did not enter the argument of the correlation functions as in the QLT limit. Consequently, it was possible to specify explicitly how the correlation functions decay in both time and space. Whereas the Eulerian correlation functions were assumed to decay exponentially according to Poisson statistics on short timescales as before, the correlation functions were specified to decay locally according to a Gaussian distribution. In the end the nonlinear Eulerian correlation function approach also yielded a normal diffusion-advection equation because in both cases the disturbed particle trajectories were assumed to be statistically Markovian (without memory effects) and spatially local.
From these focused transport-based kinetic diffusion-advection equations for energetic particle interaction with SMFRs we derived Parker-type diffusion-advection transport equations by imposing the limit of nearly isotropic distributions functions on large spatial scales (Zank et al. 2014;le Roux et al. 2015le Roux et al. , 2018. The challenge is to generalize these fundamental derivations that start with the underlying microscopic physics of particle interaction with SMFRs to also include anomalous diffusion. In other words, can we derive general kinetic fractional diffusion-advection equations similar to Equation (2) from first principles?
In a pioneering publication, Sanchez et al. (2006) showed how one can derive fractional transport equations from the basic underlying microscopic physics of particle interaction with the turbulent flow to provide a general description of anomalous diffusive transport of the tracer density in incompressible fluid turbulence using the continuity equation as a starting point. Sanchez et al. (2006) discussed how one can, by staying within the framework of the standard Lagrangian correlation function closure approach, transition from normal diffusion as a result of QLT to anomalous diffusion using a nonlinear renormalization theory. In this approach the explicit QLT assumption of the tracer particles perceiving decorrelated fluid turbulence on short timescales (Markovian assumption) was relaxed to bring in memory effects. However, the resulting closed transport equation was limited to fractional partial time derivatives with normal partial spatial derivatives because the limitation of the classical Lagrangian approach, that tracer particles perceive decorrelated fluid turbulence on local spatial scales, has not been addressed. Sanchez et al. (2006) showed how one can overcome this limited description of anomalous transport by not integrating over the tracer particle Green's function propagator, thus avoiding a Lagrangian closure approach where trajectory information enters the turbulence correlation function. In this essentially Eulerian correlation function closure approach the ensemble-averaged triple product of the second order turbulence correlation function and the Green's function propagator was dealt with by applying an integration technique to the disturbed tracer particle trajectories along the turbulent flow, resulting in a general closed transport equation with both fractional time and spatial derivatives. In this approach a b < = < H 0 2 (the index of the fractional time derivative) and a < < 0 2(the index of the fractional spatial derivative). For b < 1 and a < 2 the derived fractional transport equation is consistent with the general anomalous diffusion equation based on the separable CTRW model (see Equation (1)). For b < < 1 2 the fractional transport equation appears to go beyond the limit of a CTRW-based model (although see Luchko (2013)) and the equation was classified as being a fractional generalization of the wave equation if a < < 1 2 (e.g., Momani 2006).
In this paper, using the methods of Sanchez et al. (2006) briefly discussed above as a guide, we present some of the details of how one can derive a general multidimensional fractional diffusion-advection equation that links the macroscopic fractional transport statistics of energetic particles in ordinary, pitch-angle, and momentum space to the microscopic physics of energetic particle interaction with dynamic SMFRs using the standard focused transport equation as a basis. The outline of the paper is as follows: In Section 3 we present the application of the usual QLT-type perturbation analysis method to the standard focused transport equation. This is followed by a solution of the nonlinear evolution equation of the disturbed particle distribution function (including the nonlinear quadratic terms in fluctuating quantities) on relatively short timescales using a Green's function approach. In Section 4, we discuss three approaches to achieve a closed ensemble-averaged focused transport equation in the form of a kinetic diffusionadvection equation. In Section 4.1, we show how the standard QLT Lagrangian correlation function closure approach results in a normal kinetic diffusion-advection equation, in Section 4.2 how a nonlinear extension of the QLT Lagrangian correlation function closure approach yields a kinetic fractional diffusionadvection equation with anomalous diffusion terms where only partial time derivatives are fractional, and in Section 4.3 how a different nonlinear Eulerian-based correlation function closure approach produces a more general kinetic fractional diffusionadvection equation that allows for anomalous diffusion terms where both the time and the spatial, momentum, and cosine of the pitch-angle partial derivatives are fractional. This is followed by Section 4.4, which includes a discussion on how to turn the general fractional diffusion-advection equation into a general fractional Fokker-Planck equation, and by the Summary in Section 5.

Perturbation Analysis of the Focused Transport Equation
Consider the standard focused transport equation (Skilling 1975 for the evolution of the nearly gyrotropic gyrophase-averaged in a high conductivity plasma as a function of particle position x, time t (both measured in the plasma flow frame), particle momentum p, and the cosine of the particle pitch angle μ (both determined in the stationary (observer) frame). In Equation (5), á ñ f V g is the fixed-frame gyrophase-averaged guiding center velocity given by The first term models the change in particle pitch angle when particles experiences the magnetic mirroring/focusing force. The term on the right-hand side of Equation (5) represents energetic particle scattering in pitch-angle space due to interactions with low-frequency Alfvén wave turbulence. The pitch-angle diffusion coefficient mm D A is determined by standard QLT for gyro-resonant particle interaction with parallel propagating Alfvén waves (Schlickeiser 1989). This implies interaction with relatively small-scale wavelengths of the order of the particle gyroradius r g , which at least for energetic ions in the inner heliosphere, belong partly to the smaller scales in the inertial range of plasma wave turbulence. It could be argued that one should also take into account MHD waves that propagate at oblique angles relative to the background magnetic field. Here, we follow a simplified model for lowfrequency turbulence in the solar wind based on nearly incompressible MHD turbulence theory (e.g., Zank et al. 2017). The theory predicts that turbulence in the solar wind is mainly composed of a dominating quasi-2D MHD turbulence component containing SMFRs and a minor parallel propagating Alfvén wave component. This is supported by observations in the slow solar wind near Earth (e.g., Matthaeus et al. 1990;Bieber et al. 1996;Greco et al. 2009; and MHD turbulence simulations (e.g., Shebalin et al. 1983;Dmitruk et al. 2004). To model statistically energetic particle interaction with a collection of intermediate-scale SMFRs, with cross-sectional scales belonging to the large-scale end of the turbulence inertial range, the variables in Equation (5) are decomposed into mean (ensemble-averaged) and random fluctuating parts following the standard quasi-linear perturbation analysis technique. Accordingly, Note that we specify random fluctuations d mm D A generated by the intermediate-scale SMFR medium in the pitch-angle diffusion coefficient due to interaction with gyroscale Alfvén waves. We do allow for a distinction between the mean spatial, momentum, and pitch-angle rates of change caused by particle interaction with the background nonuniform plasma flow and magnetic field, and the average rates of change induced by the mean nonuniform plasma flow and magnetic field (mean flow and field gradients) that potentially exist in a collection of coherent SMFRs, by defining where the subscript "•" denotes the background plasma flow and magnetic field contribution and subscript "I" the contribution from the average nonuniform plasma flow and magnetic field from SMFR structures. It is through this distinction that our perturbation analysis of energetic particle interaction with coherent SMFR structures differs from that of particle interaction with classical random wave turbulence because in the latter case the mean nonuniform wave turbulence flow and magnetic field (mean wave flow and field gradients) are assumed to be zero so that there is no wave contribution to the particle's mean spatial, momentum, and pitch-angle rates of change. Thus, for example, although the average SMFR flow velocity d á ñ = U 0 I , the average SMFR shear flow gradient is not zero in our approach. This enables us to model a situation where multiple magnetic islands are simultaneously contracting, such as expected during early-time magnetic reconnection at a large-scale current sheet.
After taking the ensemble average of Equation (5) to model particle interaction with SMFRs on large spatial scales and long + á ñ To simplify the analysis, we assume that particle transport on intermediate scales is predominantly affected by SMFR structures so that the effect of the gyroscale Alfvén waves on particles can be neglected. Thus, gyroscale Alfvén waves are considered to be so weak that their contribution to particle transport becomes only detectable on large spatial scales and long timescales in Equation (11). Furthermore, we neglect the ensemble-averaged quadratic (nonlinear) terms in fluctuating quantities (last four terms in Equation (12)). By later iterating the solution for df to include the neglected terms and inserting the iterated solution back into Equation (11) to achieve closure, one can determine if their neglect was justified (Balescu 1995;Sanchez et al. 2006). The simplified df equation reads as 13) contains the exact propagator of particles through SMFR structures because the nonlinear terms containing information about the effect of the SMFR on particle transport have been retained (the transport coefficients are complete as defined in Equation (9).
In order to solve Equation (13), we first need to solve the equation for the Green's function that provides the exact particle trajectory in ordinary, momentum, and pitch-angle space inside an SMFR structure. This equation is The solution for the Green's function can be expressed as . The Green's function describes the integral solution of the exact particle trajectory in ordinary, momentum, and pitch-angle space backwards in time from t to ¢ t ( > ¢ t t ) during a particle's interaction with an SMFR structure in the plasma. This enables us to write the solution of Equation (13) in terms of the Green's function, which is where G and Q are expressed by Equations (16) and (14), respectively. In the solution for -¢ to follow particle trajectories backward in time from the final time t to the initial time t 0 . By taking into account the range of initial values m x p , , 0 0 0 relative to finite final values x, p, and μ, and by choosing = t 0 0 , we can express the solution for df with expanded boundaries as The upper boundary for the momentum integral is the result of specifying the minimum physically possible initial momentum = p 0 0 . Since the probability of such particles reaching a final momentum p of energetic particles is low, one can extend the upper boundary to infinity. No such extension of the integral boundaries of the μ integral can be imposed because there is a significant probability of finding energetic particles at m = 1. Thus, finally We can do the next iteration of the df solution by inserting the current df solution (Equation (17)) back into the evolution equation (Equation (12)) for df , where we keep the neglected ensemble-averaged nonlinear terms. After inserting the new solution for df into the evolution equation (Equation (11)) for á ñ f , we find that the previously neglected nonlinear terms in Equation (12) feature in a series of terms that are proportional to in the lowest order iteration of the df solution, one can justify having neglected the ensemble-averaged nonlinear terms in Equation (12) (see discussions in Balescu 1995;Sanchez et al. 2006). Finally, note that after inserting the solution for df in the nonlinear closure terms (quadratic terms in fluctuating quantities) of Equation (11), one gets correlation function terms such as ( · )( ) ( , , which can be neglected at late times  t t cor , where t cor is the characteristic timescale on which these correlation functions decay.

Different Approaches to the Closure Terms
The solution (Equation (17)) for df is inserted into the nonlinear closure terms of Equation (11) to accomplish a closed ensemble-averaged focused transport equation for the evolution of the average energetic particle distribution function á ñ f on large spatial scales and long timescales after interaction with numerous SMFRs. How one approaches the closure terms will determine whether one ends up modeling energetic particle interaction with SMFRs in terms of Fokker-Planck terms for classical diffusion or for anomalous diffusion involving fractional derivatives (Sanchez et al. 2006).

Standard Quasi-linear Lagrangian Correlation Function Closure
To illustrate the standard QLT closure approach, we consider a term from the first closure term in the last line of Equation (11) after the integrals containing the delta functions from the Green' function propagator (Equation (16)) in the df solution were executed:  In classical QLT one assumes weak SMFR turbulence by dropping the effect of SMFRs on the particle trajectory in x, p, and μ space (the nonlinear terms in fluctuating quantities on the left-hand side of Equation (13) are neglected). Consequently, ¶á ñ ¶ f x j is independent of SMFR effects and can be removed from the ensemble average in Equation (18), which enables formation of the disturbed guiding center trajectory correlation function d d á ñ V V gi gj induced by interaction with SMFR structures. This correlation function's decay occurs on a correlation timescale t cor along the particle's trajectory in x, p, and μ space through SMFR structures, which in the QLT limit is predominantly determined by the remaining background plasma flow and fields. Our assumption is that when a particle's guiding center exits an SMFR structure, it will see a decorrelated magnetic field and plasma flow that can abruptly change the particle's guiding center trajectory. Thus, in QLT the energetic particle guiding center's SMFR crossing time, based on a trajectory along the background/ guide magnetic field undisturbed by the SMFR structure's weak magnetic island component, can be viewed as a measure of t cor (see Equation (A1) in Appendix A for an estimate of t cor ). If the background plasma flow and magnetic field are taken to be weakly nonuniform and the background parallel electric field is assumed to be negligible, variation of particle momentum and pitch angle is weak. Then, in Equations (18) and (19) the particle trajectory expressions in the arguments of the correlation function and in ¶á ñ ¶ f is the unit vector along the background magnetic field and U 0 is the background plasma flow velocity). Finally, it is specified that in Equation (18) gi gj g 0 decays implicitly along the particle's guiding center trajectories through SMFR structures on both short time and local spatial scales. This implies that the disturbed particle guiding center trajectories and related SMFR turbulence are statistically Markovian (without memory effects) and spatially local. This is captured in the approximation that ¶á ñ ¶ f , , , j , which removes any possibility for modeling anomalous fractional transport. Accordingly, the closure term (Equation (18) where C ij L is the Lagrangian correlation function defined as in terms of a volume integral of the locally decaying disturbed guiding center trajectory correlation function along the particle's approximately undisturbed guiding center trajectory, thus removing any further explicit reference to the spatial trajectory. In , which refers to the guiding center velocity disturbed weakly by an SMFR structure with a weak twist or magnetic island component dB I perpendicular to a relatively strong, weakly nonuniform axial or guide ). The guide field component is assumed to be predominantly determined by the background magnetic field (see Hu et al. 2018) for some observational support). It is then assumed that SMFRs are statistically homogeneous and stationary (Ruffolo et al. 2003). Furthermore, consistent with the assumed short characteristic correlation (SMFR crossing) time t cor along the undisturbed guiding center trajectory, the disturbed guiding center trajectory correlation function, and by implication the correlation function of the associated SMFR flow and magnetic field generating the disturbed trajectories, are modeled explicitly to decay rapidly as an exponential function in time. Accordingly, we specify This expression refers to the probability of a certain correlation time (undisturbed energetic guiding center's SMFR crossing time), which can be related to a Poisson probability distribution for the waiting time (escape time from SMFR) given by There is only one value for the characteristic correlation time t cor because the twist component of SMFRs in QLT is simply too weak to disturb the guiding center trajectory significantly and thus to transmit information about the intermittent nature of the plasma medium in which SMFRs occur to the guiding center trajectory. After execution of the time integral over a timescale of  t t cor , the final result is a diffusion term for particle transport in the 2D magnetic island plane perpendicular to the background/guide magnetic field where the expression of the perpendicular diffusion coefficient is In Equation (26), s c I is the normalized cross helicity and r A I is the Alfvén ratio of the magnetic island or twist component, whereas V A0 is the Alfvén speed associated with the background magnetic field B 0 assumed to be the guide field of SMFRs. Equation (26) was derived assuming that statistically the SMFR magnetic island and associated flow component in the 2D plane perpendicular to B 0 is axisymmetric around the background/guide magnetic field. The main point of this section is that the standard QLT Lagrangian correlation function closure approach results in normal diffusion terms in x, p, and μ space because of the assumption that the disturbed particle trajectory correlation functions decay on short time and local spatial scales along the guiding center trajectory, the latter specified to be undisturbed by SMFRs on the correlation function decay timescale t cor .

A Classical Diffusion-advection or Fokker-Planck Equation
After closing the rest of the nonlinear terms on the right-hand side in the ensemble-averaged focused transport equation (Equation (11)) in the same fashion as was done for the term · d d - á ñ V f g above, and for simplicity ignoring the nonlinear closure terms associated with fluctuations in the Alfvén wave scattering frequency n sc A , the closed ensemble-averaged focused transport in conservation form reads as where locally the magnetic field coordinate system with unit vectors ( ) = e e e b , , 1 2 3 0 is defined so that the background/ guide magnetic field unit vector  b e 0 3 . Since the background magnetic field is assumed to be weakly nonuniform, the magnetic field coordinate system is assumed to be locally equivalent to a Cartesian coordinate system so that e e e e e e , , , , , , x y z 1 2 0 1 2 3 . Accordingly, the position coordinates ( . Conservation form was achieved by limiting ourselves to fast but nonrelativistic particle speeds, and by assuming that the average magnetic mirroring force and the mean divergence of the plasma flow in SMFRs are small. In Equation (27), the combination of transport terms is recognizable as the standard focused transport equation (Skilling 1975;Isenberg 1987) for energetic particle advective transport in the nonuniform background plasma flow U 0 and magnetic field B 0 in conservation form with a scattering term on the right-hand side due to particle gyro-resonant interaction with parallel propagating Alfvén waves. In Equation (28), In Equation (29), · E b 0 0 refers to the large-scale parallel electric field, which can be put to zero because on large scales in a high conductivity plasma the dominating electric field component The expressions for á ñ dp dt I and m á ñ d dt I , referring to coherent advective energetic particle transport in p and μ space in response to the average nonuniform flow and magnetic field of SMFRs, can be found in Appendix B (Equations (B2)-(B5)). Furthermore, the diffusion coefficients can be expressed as  (25) and (26). The detailed expressions for these diffusion coefficients can be constructed by inserting the expressions for the ensemble averages of the associated variances of the disturbed particle trajectories in response to fluctuations in the nonuniform SMFR flow and magnetic field presented in Appendix B (Equations (B6)-(B13)), and the expression for the energetic particle correlation time t cor (the characteristic time over which particles see coherent SMFR magnetic fields) in the QLT limit presented in Appendix A (Equation (A1)).
In conclusion, the standard QLT approach results in an ensemble-averaged focused transport equation that is a multidimensional diffusion-advection equation with nonuniform diffusion coefficients for normal diffusion in x, p, and μ space that is at the same time a classical Fokker-Planck equation. The key assumptions for this outcome are that the disturbed particle trajectory and associated SMFR correlations functions decay explicitly exponentially on short timescales, and implicitly on local spatial scales. In other words, the disturbed particle trajectories in x, p, and μ space and associated SMFRs are assumed to be statistically (in an ensemble-averaged sense) Markovian (no memory effects) and local spatially.

Nonlinear Lagrangian Correlation Function Closure
Following the approach by Sanchez et al. (2006), we consider again the closure term expressed in Lagrangian form in Equations (18) and (19) (the exact particle trajectories in x, p, and μ space for particle propagation through SMFR structures are specified in the arguments of the variables). Different from the QLT approach, these exact expressions are now retained. Thus, individual SMFRs affect particle transport in x, p, and μ space via the magnetic island (twist) component on short timescales through the quadratic (nonlinear) terms in fluctuating quantities present in the evolution equation for df (Equation (13)). A second difference compared to QLT is that the disturbed energetic particle guiding center correlation function , is allowed to decay on long timescales with a significant probability. However, the QLT restriction of the local decay of the correlation function along the particle trajectories in x, p, and μ space is retained. This is expressed implicitly by assuming that , , j , where, different from the QLT approach, the now non-negligible time lag in the average particle distribution is retained. It is this feature that brings into play fractional time derivatives as the agent for anomalous particle transport (Sanchez et al. 2006).
Once again, ¶á ñ ¶ f x j can be removed from the ensemble average in Equation (18) because its argument is independent of the effects of the fluctuating SMFR magnetic island component. Then, where C ij L is the Lagrangian correlation function defined as in terms of integrals of the disturbed guiding center trajectory correlation function that decays locally along the exact particle trajectories , through SMFRs (Equation (19)), thus removing any explicit reference to these trajectories.
After assuming that SMFRs are statistically homogeneous and stationary (see discussion by Ruffolo et al. (2003)), and having previously specified implicitly that there is a significant probability of the decay of the disturbed guiding center trajectory Lagrangian correlation function over longer time intervals due to particle trapping in SMFRs, it makes sense to replace the former rapid exponential decay in QLT with a more gradual power-law decay according to the function 1 . This power-law dependence can be derived by assuming that SMFRs, defined to have cross-sectional scales coinciding with scales associated with the inertial range of MHD turbulence, are always imbedded in a nonuniform solar wind plasma medium containing, for example, intermittently larger-scale magnetic flux ropes, magnetic clouds, streams, interaction regions, and shocks. Consequently, SMFRs have intermittently intense magnetic fields which efficiently trap energetic particles, resulting in relatively long correlation times for the decay of the SMFR-disturbed guiding center trajectory correlation function. Thus, instead of there being only one characteristic correlation time t cor for the decay of the correlation function when particles interact with SMFRs in a uniform plasma medium, particle interaction with SMFRs in an intermittent plasma medium generates a distribution of characteristic t cor values with an average value t á ñ cor . It can be shown that the probability of particle escape from SMFRs over a time interval Dt, when averaged over a nonuniform plasma medium generating a distribution of characteristic correlation frequencies times ( ) n P cor , can be expressed as (e.g., Lepreti et al. 2001;le Roux et al. 2010) , 34 t cor 0 cor cor cor 2 cor cor where n t = 1 cor cor . By specifying a general distribution of characteristic correlation frequencies , 35 cor cor cor cor cor cor we find that the probability of particle escape from SMFRs distribution after a time interval Dt is described by the distribution function The total probability of particle escape over a time interval cor . Consequently, the total probability of particle trapping in this time interval is  (31), we get

Subdiffusion
Following the approach of Sanchez et al. (2006), we apply integration by parts to the time integral in Equation (38) and assume late times (  t á ñ t 1 cor ) before expressing the time integral partly in terms of a fractional partial time derivative for subdiffusive transport. Then, , the normal perpendicular coefficient kÎ is expressed in Equations (25) and (26), and the perpendicular transport coefficient for fractional diffusion has the expression where a = 2 is the fractional index of the spatial derivative to indicate a normal second order derivative, β is the fractional index of the time derivative, and b - is the left Riemann-Liouville fractional differential operator, which can be expressed as The first term in Equation (39) reproduces the normal perpendicular diffusion term that we derived above for the case of standard QLT Lagrangian closure (Equations (24)- (26)) that is in competition with the second term, the anomalous diffusion term with a fractional time derivative. The problem with this result is that the fractional diffusion term has a negative sign that needs to be corrected. Sanchez et al. (2006) suggested a way to overcome this problem by modifying the Lagrangian correlation function (Equation (33)) to be , which reverses the sign of the fractional term. The result is The fractional diffusion term now has a positive sign, while the normal diffusion term will have a positive sign when 3 2 , only fractional diffusion occurs. As shown by Sanchez et al. (2006), one can use dimensional analysis to estimate on what timescale the fractional diffusion term dominates the normal diffusion term when it is positive. By specifying . This indicates that for | | d b b values sufficiently large, fractional diffusion can dominate normal diffusion for a long time period before normal diffusion takes over, which is the assumption we will make for the purposes of deriving the complete ensemble-averaged focused transport equation for anomalous diffusion.

Superdiffusion
To model superdiffusion in terms of a fractional time derivative requires modifying Equation (38) to apply for b < < 1 2. This can be accomplished by adjusting the Lagrangian correlation function in Equation (33) to be positive , , , , 46 where the perpendicular fractional transport coefficient is given by Equation (40), and b - represents the left Riemann-Liouville fractional integral (left Riemann-Liouville fractional operator with a negative fractional index because b -< 1 0 ), which we express as

A Fractional Diffusion-advection Equation Limited to Fractional Time Derivatives
After closing the rest of the nonlinear terms on the right-hand side in the ensemble-averaged focused transport equation (Equation (11)) in the same fashion as was done for the term · d d - á ñ V f g above, and for simplicity ignoring the nonlinear terms associated with fluctuations in the Alfvén wave scattering frequency n sc A , the closed ensemble-averaged focused transport equation in conservation form, valid for b < < 0 2 and a = 2, becomes  The terms on the left-hand side and the first term on the righthand side are in exact agreement with those in the classical Fokker-Planck equation (Equation (27)). The rest of the terms on the right-hand side refer to new scattering terms for anomalous diffusion containing fractional time derivatives but normal x i⊥ , p, and μ derivatives. In these terms, b - represents the left Riemann-Liouville fractional operator given by Equation (41) for subdiffusion ( b < < 0 1 ) and denotes the left Riemann-Liouville integral of Equation (47) in the case of superdiffusion ( b < < 1 2). Furthermore, the anomalous diffusion coefficients, which now depend on the fractional index β, can be expressed compactly as  The expressions for the variance in the disturbed particle trajectories in response to fluctuations in SMFR properties listed in Equation (49) can be found in Appendix B (Equations (B6)-(B13)), while the expression for the average energetic particle correlation time t á ñ cor related to particle trapping in SMFRs is presented in Appendix A (Equations (A5) and (A7)).
The closed ensemble-averaged focused transport equation (Equation (48)), being valid for b < < 0 2 and for a = 2, yields subdiffusive transport for b < < 0 1(Hurst fractional index b < = < H 0 2 1 2), and superdiffusive transport for b < < 1 2 (Hurst fractional index b < = < H 1 2 2 1). This transport equation can be related to simplified fractional transport equations in the literature by assuming constant diffusion coefficients, and neglecting mixed derivative terms as well as the term for particle scattering on Alfvén waves (first term on the right-hand side). Then, we get By restricting ourselves to a fractional time index of b < < 0 1 (subdiffusion), this equation can be identified as a threedimensional version of the one-dimensional fractional diffusion-advection equation for modeling subdiffusive transport discussed in Metzler & Klafter (2000) (see also the Introduction in Section 1). The latter equation can be interpreted as a general anomalous diffusion equation based on the separable CTRW model (see Equation (1) 1and a = 2) and extended to include advection. This connection becomes clear if one neglects the advective terms in Equation (50) and multiplies the equation with the fractional differential operator  (1)) in the subdiffusive limit (e.g., Metzler & Klafter 2000;Sanchez et al. 2006). The lefthand side of Equation (51) features the left Caputo fractional differential operator, which we express as If one neglects the advection terms in Equation (50) and restricts oneself to superdiffusive transport ( b < < 1 2, a = 2), the fractional transport equation becomes  (1)) if in the latter we put a = 2 and extend its range to b < < 1 2. There are indications that this extension is possible (Luchko 2013). However, in our case it is not possible to cast Equation (53) in the form of the typical time-fractional wave equation to make potentially the connection to a CTRW model (Sanchez et al. 2006).
In this nonlinear Lagrangian correlation function closure limit, it is the combination of particle trajectories disturbed significantly by SMFRs inx and μ space on short timescales with the decay of the disturbed particle trajectory correlation functions with enhanced probability on long timescales that enables one to model particle trapping inside SFMR structures. The trapped particle's guiding center, and its associated pitch angle, can thus be modeled to be strongly but coherently disturbed when it experiences magnetic mirroring between the stronger magnetic fields at magnetic island midpoints (Guidoni et al. 2016). When viewed over a bouncing period timescale, however, both the particle's guiding center position and the particle pitch-angle change little during trapping. When the guiding center finally escapes across the SMFR boundary, its trajectory in both ordinary and pitch-angle space changes abruptly as shown in the simulations of Drake et al. (2006), e.g., resulting in the decay of the disturbed particle guiding center and pitch-angle trajectory correlation functions. From the viewpoint of a CTRW process we have a relatively long waiting time before the particle guiding center executes a spatial and pitch-angle step to exit the SMFR with the spatial step having a maximum size capped by the size of the SMFR. One can argue that SMFRs are by definition small structures, so the maximum spatial step size should be small. If the magnetic island component of the SMFR is not too strong, the step size in pitch angle might also be relatively small. This is consistent with the nonlinear Lagrangian closure approach limitation of local decay of the disturbed particle correlation function inx and μ space. From this perspective it makes sense to model anomalous spatial diffusion perpendicular to the magnetic field and anomalous diffusion in pitch angle as subdiffusive processes with a fractional time derivative (b < 1), while retaining normal spatial and pitch derivatives.
The key assumption that made fractional transport possible in the nonlinear Lagrangian closure limit was that the disturbed particle trajectory correlation functions, generated by particle interaction with SMFRs, decay as power laws with time to include long decay timescales appropriate for modeling particle trapping in SMFR structures, as discussed above. However, only the partial time derivative became fractional. All other fractional derivatives remained normal. This limitation can be traced to retaining the QLT assumption that the decay of the correlation functions occurs locally inx , p, and μ space. Thus, the disturbed particle trajectories are assumed to be statistically non-Markovian (with memory effects), but local inx , p, and μ space. The main drawback of the latter is that all the different transport processes have the same fractional index of b < < 0 2 (the same Hurst fractional index b < = < H 0 2 2). The needed flexibility to model transport processes with different fractional indexes does not exist. For instance, if one wants to model energetic particle trapping in SMFRs as a subdiffusive spatial transport process across the background/guide magnetic field ( ) b < < 0 1, as discussed above, by default particle transport in p and μ space will also be subdiffusive with the same fractional index. Then, acceleration of trapped particles inside strong dynamic SMFRs, which one might expect to have a more systematic character to be superdiffusive in an ensemble-averaged sense (Drake et al. 2006;Isliker et al. 2017Isliker et al. , 2019, will be underestimated. More flexibility can only be achieved if the limitation of local decay of the correlation functions inx , p, and μ space is relaxed, which is the subject of the next section.

Nonlinear Eulerian Closure
Using the approach by Sanchez et al. (2006) as a guideline, we will now show how one can approach the nonlinear closure problem to derive a more general closed ensemble-averaged fractional focused transport equation in which the anomalous diffusion terms contain both fractional time and fractionalx , p, and μ derivatives. A key element in this approach is a Eulerian treatment of the disturbed particle trajectory correlation functions in the nonlinear closure terms in Equation (11). This is accomplished when, after inserting the solution (Equation (17)) for df (that includes the exact disturbed trajectories for particle interaction with SMFRs) into the closure terms in Equation (11), we do not execute the integrals over the space, momentum, and pitch-angle intervals by applying the delta functions from the Green's function propagator (Equation (16)). In this way we avoid the situation where only a time interval integral remains in which the disturbed particle trajectory information entered the argument of the disturbed particle trajectory correlation functions (Lagrangian correlation functions).
Consider as an example a term from the second closure term in the last line of Equation (11) after having inserted the general solution for df as stated in Equation (17): where the Green's function (Equation (16)) has been simplified to In the simplified Green's function the particle trajectory integrals of Equation (16) have been simplified by introducing the averages of á ñ á ñ á ñ f f f V dp dt dp dt , , g over the time interval t cor , the characteristic timescale over which the disturbed particle trajectory correlation function in Equation (54) decays, which is assumed to correlate with the characteristic timescale for particle escape across SMFR boundaries. Thus, average rates of change in the particle guiding center position, momentum, and cosine of the pitch angle are used to characterize the particle's coherently disturbed trajectory inside each coherent SMFR structure. This averaging procedure can be viewed as simplified version of the approach followed by Sanchez et al. (2006), which appears to involve numerical integration to deal with the particle trajectory integrals as given in Equation (16). In what follows the time average notation áñ t cor will be neglected.
The ensemble average in Equation (54) should be understood as an integral over many realizations of disturbed trajectories of particles in x, p, and μ space induced by interaction with SMFRs structures so that ( ) ( ) áñ = áááñ ñ ñ d d m d dp dt d dt V g . In the closure term (Equation (54)), we need only to specify explicitly the ensemble average over all possible disturbed trajectories in momentum space, which can be expressed as x x x dp dt t p dp dt t t p p dp dt p d dp dt dp dt t p dp dt t t p p P dp dt dp dt t p dp dt t d dp dt P dp dt , , , where ( ( )) d P dp dt is the probability of a certain disturbed particle trajectory in momentum space during an SMFR encounter, which is defined so that the total probability for all disturbed trajectories in momentum space is ( ) ( ( )) ò d d = -¥ ¥ d dp dt P dp dt 1 I I . The delta function in Equation (56) reflects the fact that we did the replacement ( ) d á ñ = á ñ + f dp dt dp dt dp dt in the corresponding delta function in Equation (55) to distinguish between the particle's mean rate of momentum change á ñ dp dt generated by the combination of the nonuniform background plasma flow and magnetic field with the nonuniform SMFR plasma flow and magnetic field, and the fluctuation in the rate of momentum change ( ) d dp dt produced by disturbances in the nonuniform SMFR plasma flow and magnetic field (see Equations (9) and (10)).
Assuming that on macroscopic scales anomalous transport of particles occur in p space when interacting microscopically with numerous SMFRs, we define is the Hurst exponent for fractional transport and ( ) D ¢ = D -á ñ D p p dp dt t . Therefore, | á á ñ -á f dp dt dp ñ µ Ddp dt t H 1 . One can rescale | ( )| d á ñ dp dt to become statistically independent of Dt (Sanchez et al. 2006 Accordingly, we define the rescaled quantities (˜) d = dp dt ( )( ) d Ddp dt t H 1 and ( ( )) ( ) ( (˜)) d d = D -P dp dt t P dp dt H 1 , which will be introduced when in Equation (56) the transformation ( ) (˜) d d  dp dt dp dt is performed further below. The next step is to do a Taylor series expansion of ) m x t p , , , , which compactly expressed, is given by Then, x x dp dt t p dp dt t t p p dp dt t p H H dp dt dp dt t p dp dt t p H H dp dt t p , , , where we introduced the new function ( (˜)) d F H dp dt , defined in terms of the last expression in square brackets in the second line of Equation (59). The introduction of H in the second line of Equation (59) is done to ensure consistency with the more complicated approach involving numerical integration of particle trajectories in Sanchez et al. (2006). By doing the transformation ( ) (˜) d d  dp dt dp dt I I in the ensemble-averaged expression in Equation (56) using the rescaled variables defined just below Equation (57) and inserting the expression into Equation (59), we find that x dp dt t p dp dt t t p p dp dt p H t d dp dt H dp dt dp dt P dp dt dp dt p t dp dt H H dp dt t p dp dt t P p t dp dt H H t p dp dt t P p t dp dt , , , where the average of the function Φ with respect to (˜) d dp dt I was introduced before executing the integral after the first equal sign. In the expression after the second equal sign we assumed that ( H dp dt H dp dt , , dp dt . Taking into account that | (˜)| (| | ( ) ) d á ñ = á D ¢ D ñ = a a dp dt p t H constant (independent of Dt) as discussed above (see Equation (57) and its discussion), it follows that ( (˜)) (˜) d áF ñ = d H dp dt , dp dt ( ) a F H, . With this approach, in the last line of Equation (60) we arrived at the combination of constants ( ) a F H H, , thus reproducing the combination that Sanchez et al. (2006) derived with a more sophisticated approach involving particle trajectory integration.
Upon inserting Equation (60) back into the closure term (Equation (54) The goal of expressing Equation (61) solely in terms of a fractional time and fractional momentum derivative is hampered by the presence of integrals over Dx and m D . The particle trajectories are disturbed simultaneously in x, μ, and p space during interaction with individual SMFRs, suggesting that fractional derivatives with respect to x and μ should also be introduced. Further assumptions are needed to avoid this complication. Since d 0, the ensemble averages of the delta functions can be further simplified to become t . However, after executing the integrals with the delta functions in Equation (61), the argument of ¶á ñ ¶ f p contains the trajectory information t , making it partially Lagrangian. To avoid this from occurring, the approximation that over a correlation timescale is made. This implies that the particle guiding center SMFR crossing distance parallel to guide field is much less than the distance to the SMFR structure (particles are interacting with SMFRs not located close to the origin of the spatial coordinate system). It is also assumed that m m t m -á ñ » d dt cor , indicating that the systematic variation in μ caused by the nonuniform background magnetic field and flow, and by the nonuniform SMFR magnetic field and flow, over a correlation time is overall much smaller than μ ( m á ñ d dt is small). This implies that systematic advection of particles in μ space in the final transport equation will be relatively weak because the various average gradients of the magnetic field and flow in a collection of SMFRs are relatively small. On this basis, Equation (61) simplifies to At this point we introduce the asymptotic power-law tail of a symmetric Lévy stable distribution to model the probability of a specific rescaled momentum rate of change when particle trajectories are modified by numerous encounters with individual SMFR structures. In this way, we link the basic physics of particle acceleration by SMFRs according to focused transport theory on microscopic scales to macroscopic fractional transport statistics in momentum space characterized by the indexes H and α. This approach enables specification of enhanced probabilities of long waiting times and large steps in momentum space (Lévy flights) from a CTRW perspective, which can be used to model particle trapping effects in dynamic SMFRs. Accordingly, we specify where the measure of the particle displacement in momentum space σ to the fractional power α is defined as a a a a -dp dt dp dt . 65 H cor 1 After inserting the Lévy distribution (Equation (64)), introducing the transformation ( ) D ¢ = D -á ñ D p p dp dt t , and assuming  t á ñ dp dt p 1 cor in the argument of ¶á ñ ¶ f p (weak systematic advection of energetic particles in momentum space generated in the nonuniform background and SMFR plasma flow and magnetic field so that á ñ dp dt is small), In Equation (66), the accent symbol of D ¢ p was discarded. Before expressing the time and momentum integrals as fractional derivatives, one can make some progress in specifying the unknown constant ( ) a F H, by taking the limit of classical diffusion. This can be accomplished by letting a  -2 and a b inside the integrals, and using ( ) G » -  1 for   1. From this it can be concluded that the correct standard diffusive limit will be reached by specifying beforehand that thus avoiding a blowup that appears during the limit-taking process. Then, the updated expression for ( ) a C H, is , 1 1 2 2 cos 2 . 6 9 When the standard diffusive limit is taken, one finds that ( ) = g 1 2, 2 1, but ( ) a g H, remains undetermined for anomalous transport. One can express the momentum integral naturally in terms of left and right Caputo fractional derivatives with respect to momentum by first decomposing the integral as follows: Then, 0. Switching to Riemann-Liouville fractional derivatives is a necessary step because only then can one take the momentum derivative of the fractional derivative terms as specified in the closure term (Equation (66)). The Riemann-Liouville fractional derivatives in Equation (71) can be expressed as As demonstrated in Sanchez et al. (2006), in the limit a < < 0 1 the two momentum integrals in Equation (70) must first be integrated by parts before expressing them in fractional form. After integration by parts, we get The two momentum integrals can be directly expressed in terms of left and right Riemann-Liouville fractional derivatives, which leads to which formally is the same as Equation (71), but because a < < 0 1, the Riemann-Liouville fractional derivatives have the expressions fractional derivative operator. Since the fractional index is a -< 1 0, these are really Riemann-Liouville fractional integrals. (66) with the one in Equation (38) in the section on the nonlinear Lagrangian correlation function closure approach reveals that at late times we have the same power-law integrand for Dt because

Comparison of the Dt integral in Equation
. This leads us to the conclusion that the time integral in Equation (66) can be expressed in terms of left Riemann-Liouville fractional time derivatives according to and for  To verify that this more general fractional result for the closure term is consistent with the result for the closure term in the case of the nonlinear Lagrangian approach to closure discussed in Section 4.2, where the time derivative is fractional but the momentum derivative is normal (a = 2), we let a  -2 in Equation (79)

A General Fractional Diffusion-advection Equation
After treating all the closure terms in the ensemble-averaged focused transport equation in the manner discussed above, ignoring the nonlinear terms associated with fluctuations in the Alfvén wave scattering frequency n sc A , and leaving out mixed derivative terms because they require a different approach that still needs to be determined, we arrive at a closed ensembleaveraged focused transport equation in conservation form with a general description of anomalous transport given by  2 . This equation yields subdiffusive transport for b a < 2 (Hurst fractional index b a < = < H 0 1 2) and superdiffusive transport for a b < 2 (Hurst fractional index < = H 1 2 b a < 1), keeping in mind that a < 2 and b a = < H 2. Thus, different from the nonlinear Lagrangian closure approach where b < 1 implies subdiffusion because a = 2, superdiffusion is now possible even when b < 1 by specifying a b < 2 , which provides greater flexibility. Keep in mind that the advection terms for transport in p and μ space on the lefthand side have been limited in the derivation to weak advection affects compared to the fractional diffusion terms on the righthand side. This implies, e.g., that systematic acceleration by SMFRs in this equation is restricted to be weaker than anomalous diffusive acceleration by SMFRs. The left Riemann-Liouville fractional derivative operator with respect to time is specified in Equations (77) and (78)  for a < < 1 2. Since the particle distribution is not necessarily zero at m = 1, this led to extra boundary-value terms for the particle distribution (see the bottom line in Equation (84)) after the Caputo fractional derivatives with respect to μ were transformed to Riemann-Liouville fractional derivatives for a < < 1 2 (see the discussion below Equation (71)) and when integration by parts was performed in the limit a < < 0 1 (see the beginning of Section 4.3.2). (84) can be expressed compactly as where the function ( ) a C H, is given by Equation (83), the expressions for the ensemble-averaged disturbed particle trajectories in response to statistical fluctuations in SMFR properties are presented in Appendix B, Equations (B14) and (B15), and the expressions for the average energetic particle correlation time t cor related to particle trapping in SMFRs are given in Appendix A, Equations (A5)-(A7).

The anomalous diffusion coefficients in Equation
The transport equation (Equation (84)) can be compared to simplified fractional transport equations in the literature by assuming constant transport coefficients, and neglecting the term for particle scattering on Alfvén waves (first term on the right-hand side). Then, we get x i with a fractional index α refers to a standard symmetric Riesz fractional derivative  To save space, the rest of the coordinates in the argument of á ñ f were omitted. Since exactly the same expressions hold for the Riesz fractional momentum derivative in Equation (88), we do not repeat them here. The fractional μ derivative is expressed in terms of a symmetric Riesz fractional derivative with finite left and right boundaries (Kelly et al. 2019), according to    (87), given by Equation (83), is retained.
By restricting ourselves to a fractional time index of a b < = < H 0 1 , while retaining the range of a < < 0 2 , Equation (88) can be identified as a three-dimensional version of the general one-dimensional fractional diffusion-advection equation. The latter is a unification of the one-dimensional superdiffusion-advection and subdiffusion-advection equations found in Metzler & Klafter (2000). The general fractional diffusion-advection equation itself can be interpreted as a general anomalous diffusion equation based on the separable CTRW model (see Equation (1)) extended to include advection (Metzler & Klafter 2000). This becomes clear by neglecting the advective terms in Equation (88) and then multiplying the equation with the fractional differential operator  (1)) based on the separable CTRW model (e.g., Metzler & Klafter 2000;Sanchez et al. 2006), where the left-hand side of Equation (94) features the left Caputo fractional differential operator given by Equation (52). The only difference is the presence of additional boundary-value terms that arose during the derivation of the anomalous diffusion term for transport in μ space, as discussed above. If one neglects the advection terms in Equation (88) and restricts oneself to a b < = < H 1 2 and a < < 1 2, we have a transport equation only devoted to superdiffusive transport (Hurst index b a < = < H 1 2 1) given by This equation can be thought of as a three-dimensional version of the general one-dimensional spacetime-fractional wave equation (Momani 2006;Sanchez et al. 2006), which in the limit where a b =  H 2 and a  2, reproduces the classical wave equation where all the fractional derivatives are normal second order derivatives. In the limit where a b = , Luchko (2012), e.g., refers to it as the fractional wave equation, while it has been dubbed the neutral-fractional diffusion equation by Gorenflo et al. (2000). However, different from Equation (95) (1)) if the latter's range is extended to b < < 1 2. There are indications that a connection to a CTRW model can be made in the extended range of b < < 1 2 (Luchko 2013). However, in our case it is not possible to cast Equation (95) for b < < 1 2 in the form of the typical spacetime-fractional wave equation to potentially make the connection to a CTRW model (Sanchez et al. 2006).
To summarize, we demonstrated, by using the pioneering method discussed by Sanchez et al. (2006) as a guide, how a nonlinear closure approach based on an Eulerian treatment of the disturbed particle trajectory correlation functions produced a closed ensemble-averaged focused transport equation that enabled a more flexible modeling of anomalous diffusion in x i⊥ , p, and μ space in the case of energetic particle interaction with SMFRs. Whereas the nonlinear Lagrangian correlation function closure approach yielded only fractional partial time derivatives, this approach generated both fractional partial time derivatives and fractional x i⊥ , p, and μ derivatives. The key assumption that made this possible is that, while the disturbed particle trajectory correlation functions generated by particle interaction with SMFRs are assumed to decay as a power law with time as in the nonlinear Lagrangian approach, these correlation functions are specified to also decay as power laws in x i⊥ , p, and μ space. Thus, the disturbed particle trajectories are assumed to be statistically non-Markovian (with memory effects) as in the nonlinear Lagrangian approach, but nonlocal in x i⊥ , p, and μ space too. Thus, we succeeded in deriving a closed ensemble-averaged focused transport equation in the form of a general, three-dimensional fractional diffusionadvection equation that links the macroscopic energetic particle fractional transport statistics to the microscopic physics of energetic particle interaction with SMFRs in focused transport theory.

Avoiding the Boundary-value Terms
A drawback of the ensemble-averaged focused transport equation (Equation (84)) that we derived in the form of a general fractional diffusion-advection equation is the presence of the boundary-value terms involving ( ) m á ñ  f 1 that result from the process of deriving the fractional diffusion terms in μ space. We argue that when energetic particles are trapped in SMFRs, they experience magnetic mirroring repeatedly inside SFMRs or they execute closed guiding center orbits around the magnetic island in the 2D plane perpendicular to the guide field (e.g., Ambrosiano et al. 1988;Drake et al. 2006;Guidoni et al. 2016;Du et al. 2018). We expect such trapping processes to yield small changes in μ over a bouncing or orbiting timescale (the particle tends to always return to the same μ value) until they exit the SMFR structure by making an abrupt change in pitch angle (e.g., Drake et al. 2006). This is effectively akin to a CTRW process in μ space with long waiting times between random steps in μ. Thus, particle transport in μ space can be effectively treated as a subdiffusive process. Such subdiffusive behavior can be modeled with a fractional time derivative alone so that there is no need to have a fractional μ derivative. A normal μ derivative can be recovered in Equation (84) by letting a  -2 in the m D integral before fractional μ derivatives are introduced (see the Dp integral equation (Equation (66) and its discussion below, which is relevant for the analog case of the m D integral discussed here) and we get rid of the boundary-value terms. Accordingly, we suggest the revised fractional diffusion-advection equation for b < 2 and a a < , 2  where we introduced different fractional indices for anomalous perpendicular diffusion (a 1 ) and anomalous momentum diffusion (a 2 ) to have more flexibility to model different anomalous diffusion processes with different fractional transport characteristics. We have argued above that for the problem of particle trapping and associated acceleration by SMFRs one may expect effective particle transport as subdiffusive perpendicular to the background/guide field direction, whereas anomalous momentum diffusion might be superdiffusive in an ensemble-averaged sense to reflect the more systematic character of the acceleration as reported in a number of simulations (Drake et al. 2006;Isliker et al. 2017;Du et al. 2018;Isliker et al. 2019). As before, we can get a simple general fractional diffusionadvection equation for b < 2 and a a < , 2 1 2 by assuming constant anomalous diffusion coefficients, and by neglecting the Alfvén wave scattering term. Then, , we recover the general anomalous diffusion equation

A Fractional Fokker-Planck Equation
An interesting question is whether one can relate Equation (84) to general fractional Fokker-Planck equations reported in the literature. By multiplying Equation (84) with the fractional differential operator , applying dimensional analysis to the advection terms to replace the fractional differential operator with a timescale T −( ɑH−1) , and finally applying the inverse fractional operator to the equation, one finds a general fractional Fokker-Planck equation for a b = < 2 and a < 2 given by In this form we can interpret our equation as a multidimensional general fractional Fokker-Planck equation that in one spatial dimension is similar to the equation discussed by Effenberger (2012). In one spatial dimension this equation can also be viewed as unifying the time-fractional Fokker-Planck equation (e.g., Metzler & Klafter 2000) with the space fractional Fokker-Planck equation limit (e.g., Liu et al. 2004). The equation is not as general as the fractional Fokker-Planck equation discussed by Baumann & Stenger (2017) where also the drag term features a fractional spatial derivative. Finally, by restricting a b = < H 1, and by specifying a normal second order μ derivative as above, we get

Summary
There is growing evidence from spacecraft observations that modeling of SEP acceleration and spatial transport in the solar wind needs further refinement in terms of an anomalous diffusion description. For example, it has been suggested that SEP trapping in SMFRs, resulting in perpendicular subdiffusion, is responsible for the large abrupt drops in intensity observed at 1 au for SEPs of solar flare origin (Mazur et al. 2000;Ruffolo et al. 2003;Chollet & Giacalone 2008). SEP trapping in SMFRs has also been invoked to explain an unexpected SMFR acceleration feature in the solar wind where the maximum SEP intensity amplification factor was observed to decrease with increasing particle energy instead of the opposite predicted by the SMFR transport theory for normal diffusive transport Khabarova et al. 2016). Furthermore, SEP transport ahead of interplanetary shocks has been interpreted to be superdiffusive because of observations of power-law intensity-time profiles upstream, implying superdiffusive shock acceleration of SEPs (Perri & Zimbardo 2007;Zimbardo et al. 2017). To model the latter kinetic transport phenomenon in the solar wind, it was proposed that a fractional generalization of the Parker transport equation (Litvinenko & Effenberger 2014;Zimbardo et al. 2017) be used. The fractional Parker transport equation can be viewed as a kinetic extension of the classical anomalous diffusion equation (Equation (1)) based on the separable microscopic CTRW model (Metzler & Klafter 2000) to include essential advection transport effects (Equations (2) and (3)). However, the proposed fractional Parker transport equation was not derived from first principles linking the microscopic physics of particle interaction with MHD turbulence in the solar wind to macroscopic fractional transport, and second, Parker transport equations are limited to nearly isotropic energetic particle distributions.
Consequently, we embarked on deriving a general kinetic fractional transport equation for energetic particle anomalous diffusion using the standard focused transport equation as a basis for modeling the microphysics of energetic particle interaction with SMFRs. In this way, we aimed to generalize our previous efforts in developing a kinetic transport theory for energetic particle interaction with SMFRs (Zank et al. 2014;le Roux et al. 2015le Roux et al. , 2018 to include anomalous diffusion effects, thus enabling us to address the issue of particle trapping in SMFRs. To connect the microscale dynamics of energetic particle interaction with SMFRs to macroscale fractional transport, we largely adopted and expanded the pioneering approach of Sanchez et al. (2006) to derive as a final and main result a closed ensemble-averaged focused transport equation in the form of a general, multidimensional, fractional diffusion-advection equation for anomalous diffusion in ordinary, momentum, and pitch-angle space.
In the derivation we followed the example of Sanchez et al. (2006) by treating the nonlinear closure terms, which link the microscopic physics to large-scale anomalous diffusion, in three different ways: (1) By applying the standard QLT closure approach of Lagrangian disturbed particle trajectory and associated SMFR turbulence correlation functions that decay explicitly exponentially on short timescales and implicitly on local spatial scales, we recovered the ensemble-averaged focused transport equation in the standard form of a normal diffusion-advection or a Fokker-Planck equation, which different from our previous derivations, also includes a term for pitch-angle-dependent perpendicular diffusion for particles interacting with SMFRs. In other words, the disturbed particle trajectories and associated SMFR turbulence were assumed statistically to be Markovian (without memory) and spatially local, which is a doubtful choice given the coherent nature of SMFRs, and how large they can be in the solar wind.
(2) This was followed by a nonlinear extension of the QLT method for the Lagrangian disturbed particle trajectory and associated SMFR turbulence correlation functions (Sanchez et al. 2006). The difference from QLT is that the correlation functions were specified to decay explicitly as power laws with time, thus allowing for long decay timescales appropriate for modeling particle trapping in SMFRs. However, in this mixed approach the QLT assumption of the implicit local decay of the correlation functions in ordinary, momentum, and pitch-angle space was retained (Sanchez et al. 2006). Thus, a kinetic fractional diffusion-advection equation was derived in which only the partial time derivatives were fractional, which is a direct consequence of specifying disturbed particle trajectories that are statistically non-Markovian (with memory), while remaining local in ordinary, momentum, and pitch-angle space. Since the fractional time index is b < 2, the Hurst fractional index is b < = < H 0 2 1 , thus enabling both subdiffusive transport ( ) b < 1 and superdiffusive transport ( ) b < < 1 2. The main drawback of this approach is that it is not possible to model different anomalous diffusion processes with different Hurst fractional indexes.
(3) Finally, we introduced a totally different nonlinear closure approach in which the Lagrangian correlation function formalism was replaced by a Eulerian approach (see Sanchez et al. 2006). With this approach it was possible to specify disturbed particle trajectory correlation functions that explicitly decay asymptotically as power laws in both time and in ordinary, momentum, and pitch-angle space by specifying Lévy stable distribution functions for the statistics of the disturbed particle trajectories in response to the associated SMFR turbulence. In this way, the most general case of disturbed particle trajectories that are statistically both non-Markovian and nonlocal in ordinary, momentum, and pitch-angle space were achieved, resulting in a general kinetic fractional diffusion-advection equation with anomalous diffusion terms that contained both fractional partial time derivatives with a fractional index of b < 2 and fractional partial derivatives in ordinary, momentum, and pitch-angle space with a fractional index of a < 2. This is the main result of this paper. This equation, which yields subdiffusion when b a < 2 and superdiffusion for a b < 2 , provides more flexibility in modeling anomalous diffusion phenomena occurring during energetic particle interaction with SMFRs because superdiffusion is now possible even when, e.g., b < 1.
By assuming constant transport coefficients, and by restricting b < 1, the most general fractional transport equation that we derived is clearly recognizable as a three-dimensional version of the general one-dimensional fractional diffusionadvection equation (see Equation (2)) that unifies sub-and superdiffusive advection equations found in the literature (Metzler & Klafter 2000). The fractional diffusion-advection equation itself is equivalent to the familiar general anomalous diffusion equation (Equation (1)) based on the standard microscopic separable CTRW model extended to include advection (Metzler & Klafter 2000). When our general fractional diffusion-advection equation is limited to superdiffusion by restricting both b < < 1 2 and a < < 1 2, it can be thought of as a general fractional wave equation with advection terms. However, after simplifying the equation by assuming constant transport coefficients and dropping the advection terms, our equation does not become a threedimensional version of the general fractional wave equation in one dimension in the literature (Gorenflo et al. 2000;Momani 2006;Luchko 2013) that can be linked to a separable CTRW model as noted by Sanchez et al. (2006). We also discovered that our general fractional diffusion-advection equation can be transformed into a general three-dimensional fractional Fokker-Planck equation similar to that found in one dimension in the literature (Effenberger 2012). The latter can be recognized as a unification of the time and space fractional Fokker-Planck equations (Metzler & Klafter 2000;Liu et al. 2004).
Our most general fractional diffusion-advection equation has some drawbacks and require further development in the future. The derivation of the fractional pitch-angle diffusion term resulted in extra distribution function boundary-value terms at m = 1 that cannot be dropped. The extra boundary terms can be sidestepped by expressing this term as a normal second order pitch-angle derivative while retaining a fractional time derivative. This limits the term to subdiffusive transport in μ space, which might be sufficient to account for SMFR-induced particle trapping effects on particle transport in μ space. Furthermore, mixed derivative transport terms were neglected because it is not clear yet to us how to cast them in a fractional derivative form. However, there are limits to which the neglect of mixed derivatives can be justified. In the case of a mixed derivative fractional transport term with respect to p and μ, one can argue that, if the particle distribution is nearly isotropic with respect to μ, particle acceleration by SMFRs will be predominantly described by the fractional momentum diffusion term. In addition, it has to be noted that the advection terms in μ and p space were assumed to be relatively weak to achieve a tractable fractional transport equation. This implies, e.g., that systematic acceleration of energetic particles by SMFRs due to average SMFR field and flow quantities is less important than acceleration by SMFRs involving anomalous diffusion in response to fluctuations in SMFR field and flow quantities. Finally, the expressions for the anomalous diffusion coefficients include an undetermined constant that needs to be determined from particle simulations in future work. J.A.l.R. acknowledges support from NASA grant 80NSSC19K027, and NSF-DOE grant PHY-1707247. G.P.Z. acknowledges the partial support of a NASA Parker Solar Probe contract SV4-84017, a NASA IMAP subaward under NASA contract 80GSFC19C0027, and a NASA award 80NSSC20K1783. Both authors acknowledge support from an NSF EPSCoR RII-Track-1 Cooperative Agreement OIA-1655280.

Appendix A The Energetic Particle Correlation Time
The energetic particle correlation time t cor is defined as the characteristic time it takes for the disturbed energetic particle trajectory correlation functions in x, p, and μ space, arising from particle interaction with coherent SMFR structures, to become decorrelated along guiding center trajectories through these structures. We assume that decorrelation sets in when energetic particles exit SMFRs so that the decorrelation time is determined by the residence time in or the escape time from SMFRs, which in turn depends on the nature of the guiding center trajectory through the SMFR structure. In the QLT limit of our transport theory (Section 4.1), when the twist or magnetic island component of SMFRs is weak compared to the background/guide field component ( 2 ), we assume that t cor is determined mainly by the characteristic particle crossing time of SMFRs along the background/guide field because the particle guiding center trajectory is approximately unaffected by the magnetic island component on short timescales. Accordingly, | | || t m = L v I cor , where || L I is the length of the SMFR along the background/guide field direction. To avoid a blowup at m = 0, we introduce the SMFR turbulence correlation time t c as a competing contribution to the particle correlation time. Its expression is d = á ñ t L U 2 c I I 2 1 2 , where dU I is the plasma flow speed in the 2D plane perpendicular to the guide field in an SMFR, and L I is the characteristic cross section of the magnetic island component of SMFRs. The particle correlation time t cor is then modeled as a competition between the characteristic crossing and turbulence correlation frequencies according to ) the turbulence correlation frequency acts as a small correction term. Thus, the flow speed in the SMFR structure weakly modifies the particle guiding center trajectory consistent with the QLT approach.
In the nonlinear limit of our transport theory, where the magnetic island component of the SMFR is allowed to be considerably stronger, the particle guiding center trajectory is assumed to be strongly affected by the magnetic island so that particle trapping can occur. It appears that particles can be magnetically mirrored at the magnetic island midpoints where the island magnetic field is strongest compared to the endpoints where it is weakest (e.g., Guidoni et al. 2016). We worked out a magnetic mirroring condition where d á ñ B I 2 refers to the average magnetic energy density at magnetic island midpoints. We also estimated a characteristic bouncing time t b between the magnetic island midpoints based on our closed ensemble-averaged focused transport theory, where In this case, the average particle trapping or escape time is determined by the mean total bouncing time, which depends on the average number of bounces n b before escape. It should be noted that this estimate might be inaccurate in situations where the average increase in particle speed is strong and/or where the average size of magnetic islands vary significantly during the trapping period. A final expression for t cor follows when it is modeled as a competition between the inverse of the mean total bouncing time and the inverse of the turbulence correlation time Particles with m m > m 2 2 do not experience bouncing between SMFR midpoints, but can still execute closed guiding center orbits in the magnetic island (e.g., Drake et al. 2006;Guidoni et al. 2016). Assuming that SMFR turbulence is on average axisymmetric around the background/guide field, one can model the average SMFR as a simple cylindrical SMFR structure with a finite axial length || L I and cross-sectional scale L I . We find that the path length along a helical magnetic field line is approximately

Appendix B The Transport Coefficients
We used perturbation theory and dimensional analysis to work out approximate, tractable expressions for the transport coefficients in the closed ensemble-averaged focused transport equations we derived that separate the effect of the background plasma flow and magnetic field from the turbulent SMFR plasma flow and magnetic field on energetic particle transport and acceleration. This was accomplished by employing the following assumptions: (1) SMFRs are quasi-helical magnetic field structures that are advected with the background plasma flow. They can be thought of as a manifestation of the perpendicular wavenumber component of an MHD shear Alfvén wave, which has a zero phase speed (Zank et al. 2017), and belong to the quasi-2D MHD turbulence component, which is thought to be the dominant MHD turbulence component in the large-scale inner heliospheric solar wind at lower helio latitudes (Matthaeus et al. 1990;Zank & Matthaeus 1992Bieber et al. 1996;Hunana & Zank 2010;Zank et al. 2017). Their magnetic field consists of a magnetic island or twist component dB I in the 2D plane perpendicular the axial or guide field component B 0 (Cartwright & Moldwin 2010). In the background plasma flow frame there is an SMFR flow also constrained to the 2D plane perpendicular to B 0 , thus limiting SMFR dynamics such as contraction and merging with neighboring SMFRs to the 2D plane of the magnetic island component (Birn et al. 1989;Dmitruk et al. 2004). The magnetic island component has statistical properties, whereas the guide field component is assumed to be dominated by the background magnetic field , which is nearly uniform on the scale of SMFRs.
(2) A strong background/guide magnetic field B 0 was assumed so that on average  d á ñ B B 1 , indicating that the background plasma flow speed is much larger than average flow speed in SMFRs. Thus, the model is best suited to the super fast magnetosonic solar wind.
(4) SMFR structures overall are elongated in the axial or guide field direction so that they depend more weakly on the spatial coordinate along B 0 than on spatial coordinates in the 2D plane perpendicular to B 0 (Weygand et al. 2011) (they are quasi-2D structures). (5) Statistically, both the magnetic field and the plasma flow of the magnetic island component of SMFRs, constrained to the 2D plane perpendicular to B 0 , manifest as axisymmetric turbulence around B 0 (the average magnetic island can be viewed as a circular structure perpendicular to B 0 . where r A I is the Alfvén ratio of the magnetic island component of SMFRs, V A0 is the Alvén speed associated with the background/ guide field component of SMFRs, L I is the characteristic crosssectional scale of SMFRs, y sin is the sine of the spiral angle of the large-scale magnetic field, | | | | = r p qB g 0 , with q the net particle charge, is the maximum size of the energetic particle gyroradius (m = 0), Z is the atomic number, A is the mass number, ( ) g v is the Lorentz factor, and d p e , is the proton(electron) inertial scale. In Equations (B2) and (B3), n á ñ I COM determines, respectively, the p and μ rate of change of energetic particles in response to the mean divergence of the SMFR flow (disregarding the μ dependence), n á ñ I SH the average p and μ rate of change due to the average component of the shear gradient in the SMFR flow along the magnetic field of the magnetic island component of SMFRs dB I , n á ñ I ACC the average p and μ rate of change in response to the mean component of the acceleration of the SMFR flow along dB I (noninertial force), n á ñ E I the average p and μ rate of change due to the mean component of turbulent motional electric field force parallel to the background/guide magnetic field B 0 , and n á ñ I REF the μ rate of change connected to the average magnetic mirroring force in the magnetic field of SMFRs.
Constants were introduced in the expressions in Equation (B4) to control the mean energetic particle acceleration rates. For example, we specify The average rate of the particle's guiding center displacement in response to the nonuniform plasma flow and magnetic field in SMFRs, determined by the average of the guiding center velocity V g I , is zero because SMFR quantities dU I and dB I without gradients are modeled as random variables, as discussed above. Therefore, The normal diffusion or Fokker-Planck scattering terms that we derived in the cases of the quasi-linear and nonlinear Lagrangian closure approaches contain the variances of the momentum and μ rate of changes for energetic particle interaction with SMFRs based on the following expressions for the random fluctuations in these rates:  For simplicity, most terms in the expressions in Equations (B9), (B11), and (B12) mixing different transport mechanisms were omitted. The neglect of some can be justified in the fast particle limit  v V 1

A0
. The variance in the rate of the particle spatial displacement is determined by fluctuations in the particle's guiding center velocity dV g I , according to the expression Simplified expressions for ensemble averages of fluctuations in the energetic particle spatial, momentum, and pitch-angle rates of change to some fractional index power that appear in the anomalous diffusion coefficients in Equation (87) of the general fractional diffusion-advection equation (Equation (84)) were derived. Consistent with the fact that Equation (84) was derived in the limit of weak coherent advection effects on energetic particles, we assumed weak advection effects in deriving the expressions for the ensemble averages. In these derivations, the inequalities | | | | | |