Investigating Superdiffusive Shock Acceleration at a Parallel Shock with a Fractional Parker Equation for Energetic-particle Interaction with Small-scale Magnetic Flux Ropes

It has been suggested before that small-scale magnetic flux rope (SMFR) structures in the solar wind can temporarily trap energetic charged particles. We present the derivation of a new fractional Parker equation for energetic-particle interaction with SMFRs from our pitch-angle-dependent fractional diffusion-advection equation that can account for such trapping effects. The latter was derived previously in le Roux & Zank from the first principles starting with the standard focused transport equation. The new equation features anomalous advection and diffusion terms. It suggests that energetic-particle parallel transport occurs with a decaying efficiency of advection effects as parallel superdiffusion becomes more dominant at late times. Parallel superdiffusion can be linked back to underlying anomalous pitch-angle transport, which might be subdiffusive during interaction with quasi-helical coherent SMFRs. We apply the new equation to time-dependent superdiffusive shock acceleration at a parallel shock. The results show that the superdiffusive-shock-acceleration timescale is fractional, the net fractional differential particle flux is conserved across the shock ignoring particle injection at the shock, and the accelerated particle spectrum at the shock converges to the familiar power-law spectrum predicted by standard steady-state diffusive-shock-acceleration theory at late times. Upstream, as parallel superdiffusion progressively dominates the advection of energetic particles, their spatial distributions decay on spatial scales that grow with time. Furthermore, superdiffusive parallel shock acceleration is found to be less efficient if parallel anomalous diffusion is more superdiffusive, while perpendicular particle escape from the shock, thought to be subdiffusive during SMFR interaction, is reduced when increasingly subdiffusive.


Introduction
The well-known Parker transport equation, for energetic charged particle (cosmic-ray) propagation and acceleration in a turbulent, high-conductivity plasma flow (Parker 1965), has been widely applied to numerous problems in space and astrophysical plasmas, including the solar wind that is our main focus. The Parker equation is appropriate for situations in which a nearly isotropic energetic-particle distribution in momentum space is maintained due to particle scattering on magnetic fluctuations in the solar wind flow. In this equation, large-scale spatial transport is considered as a competition between energetic-particle advection by the solar wind plasma flow and normal diffusive transport from particle scattering on magnetic fluctuations of low-frequency (magnetohydrodynamic, hereafter MHD) solar wind wave turbulence. Similarly, the equation sometimes include a competition between advection and normal diffusion in momentum space. The advection usually refers to coherent momentum gain or loss generated by the adiabatic compression or expansion of the large-scale solar wind flow, respectively, which is a first-order Fermi acceleration or deceleration process. Normal momentum diffusion, modeled usually as a result of particle interaction with the random-solarwind-motional electric field or velocity fluctuations of MHD wave turbulence, is also known as second-order Fermi acceleration (see review by, e.g., Zhang & Lee 2013). The Parker equation also contains the physics of standard diffusive shock acceleration (DSA) theory (see reviews by Drury 1983;Jones & Ellison 1991) used extensively for modeling energeticparticle acceleration at collisionless heliospheric shocks (see reviews by Verkhoglyadova et al. 2015;Desai & Giacalone 2016).
However, there is an increasing realization that a complete understanding of the impact of coherent, small-scale magnetic flux rope (SMFR) structures on energetic-particle transport and acceleration in the large-scale solar wind and at heliospheric shocks cannot be achieved within the frame work of normal advection and diffusion. This suggests the need for the development of a new Parker transport equation that includes anomalous advection and diffusion. For example, Mazur et al. (2000) and Chollet & Giacalone (2008) reported large abrupt drops in the counting rate of solar flare generated solar energetic particles (SEPs) at 1 astronomical unit (au) over small spatial scales of the order of the solar wind turbulence correlation length, proposing that these sudden drops are connected to the properties of low-frequency (MHD) turbulence. Ruffolo et al. (2003) suggested that these abrupt drops can be understood considering that SEPs temporarily trapped in these closed, quasi-helical SMFRs structures will maintain 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. Furthermore, more recently unusual energetic ion and electron acceleration events were observed in the solar wind that exhibit energetic-particle flux enhancements too far away from interplanetary shocks to be attributed to standard DSA theory, and that could not be explained in terms of second-order Fermi acceleration by Alfvén-wave turbulence Khabarova & Zank 2017;Malandraki et al. 2019;Zhao et al. 2019). These events were found to coincide with numerous dynamic SMFR structures and have the characteristic feature that the observed flux amplification factors increase with particle energy. New Parker-type transport equations with normal advection and diffusion transport terms were derived more recently to model the energetic-particle transport through and acceleration by dynamic SMFR structures in the largescale solar wind (Zank et al. 2014;le Roux et al. 2015le Roux et al. , 2018. The new equations were applied successfully to reproduce these flux enhancements for energetic ions (e.g., Zhao et al. 2018;Adhikari et al. 2019;le Roux et al. 2019;Zhao et al. 2019;Van Eck et al. 2021). However, SMFR ion acceleration events were identified with flux amplification factors decreasing with particle energy Khabarova et al. 2016;Van Eck et al. 2021). This anomaly might be attributed to the trapping of lower-energy particles in SMFRs  with less turbulent magnetic fields (limited scattering inside SMFRs reduces particle escape across the magnetic field) that may not escape before the SMFRs become dynamically inactive. In contrast, higher-energy particles with larger gyroradii will escape more easily to interact with more dynamic SMFRs elsewhere to experience relatively more efficient acceleration. Thus, also SMFR acceleration tentatively points to the existence of trapped suprathermal particles experiencing the subdiffusive perpendicular transport that would require a new Parker equation, which includes anomalous transport processes.
To this we can add that Zimbardo and coworkers reported in a series of papers that spacecraft observations of SEP intensitytime profiles ahead of interplanetary shocks yielded power laws instead of the exponential profiles predicted by standard DSA. They interpreted this as a sign of superdiffusive spatial SEP transport in the shock vicinity resulting in superdiffusive shock acceleration. They predicted the superdiffusive shock acceleration to be more efficient, producing a harder accelerated powerlaw spectrum for SEPs, than standard DSA (Zimbardo et al. 2006;Perri & Zimbardo 2007, 2008Zimbardo & Perri 2013;Zimbardo et al. 2017.  explained such superdiffusive transport as a product of subdiffusion in pitch-angle space arising from the observation that Alfvén-wave pitch-angle scattering times in the solar wind have a power-law distribution instead of a single value. Another possibility is that there is a strong presence of SMFRs in the vicinity of heliospheric shocks so that particle trapping effects could result in anomalous diffusive transport and shock acceleration. 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 1992, 1993Zank et al. 2017Zank et al. , 2018Zank et al. , 2020. It is plausible that SMFRs can have a strong presence ahead of interplanetary shocks. Furthermore, there is evidence that the occurrence of SMFRs peaks at large-scale current sheets in the solar wind where magnetic reconnection produces additional SMFRs that can trap energetic particles (Khabarova et al. , 2016Hu et al. 2018). Many of these current sheets occur behind traveling shocks (Khabarova & Zank 2017). Recently, a theoretical investigation was launched to determine the transmission of quasi-2D SMFRs through a perpendicular shock . The results show a strong enhancement in the magnetic energy density of the SMFR magnetic island component behind the shock, thus confirming a strong presence of SMFRs downstream of the shock as well.

Motivation
In le Roux & Zank (2021), we derived from first principles a new pitch-angle-dependent fractional diffusion-advection equation for energetic-particle interaction with dynamic SMFR structures in the solar wind on the basis of the standard focused transport equation (Skilling 1975). In this follow-up work, our first goal is to show how the fractional diffusion-advection equation can be converted into fractional Parker and Gleeson-Axford transport equations. This will be based on the assumption that during anomalous transport, even though large step sizes and long waiting times between scattering events have an enhanced probability, the probability for scattering events associated with small step sizes and short waiting times remains high so that nearly isotropic pitch-angle distributions can be maintained on large spatial scales (e.g., Effenberger 2012; Zimbardo et al. 2017). Two important features of the new fractional Parker transport equation discussed below are that parallel transport is superdiffusive when the fractional time index β < 1, and that energetic-particle advection processes are anomalous too. Consequently, we find that anomalous advection processes are increasingly dominated by superdiffusion during late times. Second, our goal is to investigate the consequences of the new fractional Parker equation for superdiffusive shock acceleration at a parallel heliospheric shock when dynamic SMFRs have a strong presence in its vicinity. It will be analyzed how anomalous DSA differs from standard DSA in regards to the shock-acceleration timescale, the accelerated particle spectrum that forms at the shock at late times, and the spatial variation of these spectra both sides of the shock.

The Existing Fractional Diffusion-Advection Equation
for Energetic-particle Interaction with SMFRs

The Assumptions
Using the standard focused transport equation (Skilling 1975) as the starting point, we used perturbation theory to derive a pitch-angle-dependent fractional diffusion-advection equation for energetic-particle interaction with SMFR structures in the largescale solar wind (for more details, see le Roux & Zank 2021). Approximate, tractable expressions for the transport coefficients in the new equation were derived that separate the effect of the background solar wind flow and magnetic field from the turbulent SMFR plasma flow and magnetic field on energeticparticle transport and acceleration. This was accomplished by employing the following assumptions: (1) SMFRs in the solar wind are quasi-helical magnetic field structures that are advected with the background plasma flow (e.g., Cartwright & Moldwin 2010;Zank et al. 2017). These structures can be thought of as a manifestation of the perpendicular wavenumber component of a MHD shear Alfvén wave that has a zero phase speed ) and belong to the quasi-2D MHD turbulence component that is considered 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 1992, 1993Bieber et al. 1996;Hunana & Zank 2010;Zank et al. 2017). Their magnetic field consists of a magnetic island or twist component δB I in the 2D plane perpendicular to the axial or guide field component B 0 (Cartwright & Moldwin 2010). The SMFR flow δU I is 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 δB I (and the associated flow δU I ) has statistical properties whereas the guide field component B 0 is assumed to be dominated by the background magnetic field B 0 , 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 Smith et al. 2006) (magnetic island or twist magnetic field component δB I with a relatively weak average magnetic energy density d á ñ B I 2 ). This does not preclude δB I from being strong outside the center of the SMFR where the guide or axial field dominates (Cartwright & Moldwin 2010) and particle trapping can occur (Guidoni et al. 2016). (3) The plasma flow of SMFRs δU I is such that d á ñ U U 1 I 2 1 2 0  indicating that the background solar wind flow speed U 0 is much larger than the average flow speed in SMFRs. Thus our SMFR 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 as hinted by the observations of 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 ). More recent state-of-the-art multispacecraft observations using four Cluster spacecraft data are consistent with the existing view of slow solar wind turbulence being dominated by a 2D turbulence (SMFR) component confined to the 2D plane perpendicular to the mean magnetic field (e.g., Narita et al. 2011). Multispacecraft data analysis has the advantage of not explicitly employing Taylor's hypothesis. However, different from our assumption, the observations of Narita et al. (2011) indicate that 2D turbulence (SMFR turbulence) is not axisymmetric in this 2D plane. However, Turner et al. (2012) concluded that the nonaxisymmetry detected in Cluster turbulence data might be due to a "preferential filtering of the modes with wavevectors mainly parallel to the solar wind flow," because data filtering in the spacecraft frame cannot be avoided, while Wicks et al. (2012) thought that the Cluster spacecraft did not spend enough time in the solar wind to address this issue. Given this uncertainty, we preferred not to deviate from the standard theoretical assumption of axisymmetric turbulence around B 0 in the derivation of transport coefficients in our transport theory because of the benefit of greater simplicity. (6) On intermediate scales (MHD scales), the magnetic reconnection and SMFR merging are determined predominantly by the turbulent component of the ideal or motional electric field δE I = δU I × δB I that is parallel to B 0 in our SMFR model. In this sense only, the electric field is interpreted as a parallel electric field. (7) SMFR magnetic fields and plasma flows are random variables with zero averages on large scales, but gradients of the SMFR flow and magnetic field are not necessarily zero on average, allowing us to model the particle acceleration by both the mean and the variance of the divergence, shear-flow gradient, and acceleration of the nonuniform SMFR flow, where appropriate. We also assume that the turbulent motional electric field might not be zero on large scales, so we distinguish between acceleration by its average and variance. In this way, we distinguish the coherent structures, such as SMFRs, statistically from traditional random wave turbulence.

The Equation
The pitch-angle-dependent fractional diffusion-advection equation for energetic-particle interaction with SMFR structures in the large-scale solar wind that we developed (le Roux & Zank 2021) on the basis of the standard focused transport equation (Skilling 1975) describes the evolution of the ensemble-averaged distribution function 〈f〉(x, p, μ, t) as a function of particle position x, particle momentum p, cosine of particle pitch-angle μ, and time t. In conservation form, the equation can be compactly presented according to where 〈f〉 is nearly independent of the particle gyrophase, and its evolution is expressed in terms of a mixed coordinate system where x and t are measured in the fixed (observer) frame whereas p and μ are determined in the noninertial moving frame at rest in the solar wind flow frame.
On the left-hand side of Equation (1), the transport terms two to four models of the divergence of the average differential energetic-particle flux in x, p, and μ-space, respectively, as a result of particle response to the nonuniform large-scale magnetic field B 0 (b 0 is the unit vector along B 0 ), plasma flow U 0 , parallel electric field E 0 · b 0 (terms with subscript "0") of the large-scale solar wind, and the average effects associated with the fluctuating magnetic field δB I , plasma flow δU I , parallel electric field δE I · b 0 belonging to dynamic SMFR structures advected with the large-scale solar wind flow (terms with subscript "I"). In these transport terms, 〈dx/dt〉 0 , 〈dp/dt〉 0 , and 〈dμ/dt〉 0 , which denote particle advection velocities in x, p, and μ-space from interaction with the large-scale solar wind, have the expressions where q is the net ion or electron charge. In Equation (2), E 0 · b 0 refers to the large-scale parallel electric field that is neglected because, in the high-conductivity solar wind plasma, the electric field is dominated by the perpendicular ideal or motional electric field E 0 ≈ − U 0 × B 0 ⊥b 0 . The total time derivative is defined as d/dt = ∂/∂t + U 0 · ∇. The expressions in Equation (2) are recognizable as those belonging to the standard focused transport equation (e.g., Skilling 1975;Isenberg 1997). Furthermore, the expressions of 〈dx/dt〉 I , 〈dp/dt〉 I , and 〈dμ/dt〉 I in Equation (1)  and r g = p/(|q|B 0 ) is the maximum charged particle gyroradius (μ = 0). In Equations (3) and (4), the mean relative p and μ rates of change of energetic particles depend on n á ñ I COM in response to the mean divergence of the SMFR flow in the 2D plane perpendicular to B 0 , n á ñ I SH due to the effect of the average component of the shear gradient in the SMFR flow in this 2D plane parallel to δB I , n á ñ I ACC in response to the mean component of the acceleration of the SMFR flow in the same 2D plane (noninertial force) parallel to δB I , and n á ñ E I due to the mean component of the turbulent motional electric field force (qδE I = − q〈δU I × δB I 〉) parallel to B 0 . These expressions, associated with the various aspects of the nonuniform plasma flow and the turbulent motional electric field in dynamic SMFRs, determine the four SMFR acceleration mechanisms in focused transport theory that were applied before to model the acceleration of energetic particles in the large-scale solar wind within the context of normal diffusion and advection in p and μ-space (e.g., le Roux et al. 2019). Furthermore, the mean rate of change in μ also depends on n á ñ I REF induced by the average SMFR magnetic mirroring force in the 2D plane perpendicular to B 0 . Note that the average particle spatial advection velocity 〈dx/dt〉 I = 0 because SMFR fluctuations δU I and δB I are treated as random random variables, as discussed above.
On the right-hand side of Equation (1) are all the scattering terms. Term one models the classical diffusion of particles in μspace due to their gyro-resonant interaction with Alfvén waves propagating along the large-scale magnetic field. This is followed by anomalous diffusion terms for particle propagation perpendicular to the large-scale magnetic field (term two), in pspace (term three), and in μ-space (last term), all due to nonresonant energetic-particle interaction with SMFRs. The expressions of the anomalous diffusion coefficients in the fractional transport terms are given by ] is the anomalous diffusion coefficient for fractional transport in μ-space. In Equation ( In the limit of classical diffusion where α = 2, h(β, α = 2) = 1 so that C(β, α = 2) = 1/2, but for α ≠ 2, h(β, α) is an unknown parameter that needs to be determined from turbulence simulations (Sanchez et al. 2006). Since in our application of this transport theory discussed further below we specify α = 2, this problem does not arise. Furthermore, in Equation (5) I I I I energetic particle p and μ-rates changes in terms of fractional averages.
In Equation (5), 〈τ cor 〉 represents the average time it takes for energetic particles inside SMFRs to encounter a decorrelated SMFR magnetic field δB I . We assume that decorrelation sets in when energetic particles exit through the sharp boundaries of coherent SMFR structures (Greco et al. 2009). Thus, 〈τ cor 〉 is an estimate of the average residence time of energetic particles in SMFR structures (escape time), based on a simple model of the guiding center trajectory through SMFRs. Our fractional transport theory is nonlinear in the sense that the particles' guiding trajectory is modified by the SMFR magnetic field (see le Roux & Zank 2021). Accordingly, assuming that the statistically averaged SMFR structure can be thought of as a simple helical magnetic field (the twist component of SMFR magnetic field is considered to be, on average, axisymmetric around the background/guide magnetic field) with an average axial length L I∥ , and crosssectional width L I , we estimated the average path length along a helical magnetic field line to be Assuming that a particle will escape the average SMFR after propagating through the structure (guiding center completing the distance s I ), one can express 〈τ cor 〉 as Our propagation model for the particle decorrelation time, where energetic particles propagate through a SMFR structure by following the associated helical magnetic field lines, can be interpreted as a weak particle trapping model. This can be understood by considering that, based on the expression from Equation (9), the path length along a helical field line for the average SMFR structure in the solar wind near 1 au is s I ≈ 1.4L I|| since L I|| /L I ≈ 2 (Weygand et al. 2011). Thus, the time for an energetic particle to propagate through a helical SFMR structure will be less than a factor of two longer compared to when the particle guiding center propagates the length of a SMFR along the guide or axial field component alone. We chose the weak particle trapping model based on the estimate that for most particle pitch angles; strong particle trapping through magnetic mirroring and the bouncing of particles inside SMFRs are unlikely in the strong guide field limit of our theory. This is apparent from our mirroring condition as suggested by observations at 1 au (Smith et al. 2006), only particles with pitch angles close to 90 o (|θ|  78 o ) can be expected to be reflected inside SMFRs at 1 au, on average. For a nearly isotropic energetic-particle pitchangle distribution, this implies that only a small fraction of particles, ∼20%, is subject to magnetic mirroring (strong trapping) inside SMFRs.
In Equation (1), each of the anomalous diffusion terms includes a left Riemann-Liouville fractional time derivative that is defined as 1 is the right Caputo fractional derivative. The strength of the anomalous diffusion is captured by the fractional indices α and β. In le Roux & Zank (2021), 0 < α < 2 and 0 < β < 2, but in this paper, we restrict ourselves for simplicity to the traditional range for superdiffusive transport 1 < α < 2 and subdiffusive transport 0 < β < 1.

A New Fractional Parker Transport Equation for
Energetic-particle Interaction with SMFRs To convert Equation (1) into a fractional Parker-type transport equation, we assume that enough scattering occurs during anomalous transport for particles to develop a nearly isotropic pitch-angle distribution on large scales as discussed above. We follow the standard procedure of doing a moment expansion of the particle distribution in μ-space where the μdependence of the distribution is modeled in terms of Legendre polynomials P n (μ) according to is the nth moment of 〈f〉. Consistent with the assumption of nearisotropy, the expansion is terminated at the second moment f 2 . First, we execute moment zero (n = 0) of Equation (1) by multiplying each each term with ) and integrating over all μ-values. The zeroth moment of Equation (1) is In Equation (14), 〈〉 μ indicates an average over all μ-values ( The expressions in Equation (15) clarifies how the divergence of the plasma flow velocity, the parallel electric field in the plasma, the acceleration of the plasma flow along the magnetic field, and the shear tensor of the plasma flow along the magnetic field in Equation (14) are decomposed into background solar wind components (terms with subscript "0") and mean SMFR components (terms with superscript "I"). The expressions for the latter is given in Equation (4). As stated above, the parallel large-scale electric field E 0 · b 0 = 0 because, appropriate for a high-conductivity collisionless plasma like the solar wind, (15), the parallel shearflow tensor of the large-scale solar wind is given by Since Equation (14) contains higher moments of the distribution function f 1 and f 2 , closure in terms of f 0 can only be achieved by finding solutions for f 1 and f 2 in terms of f 0 . For this purpose, we need to calculate the first (n = 1) moment of Equation (1) by multiplying each term with ò mm d 1 2 1 1 and integrating over μ. We find that

[ ]
In Equation (17), ν A is the scattering frequency of particles on Alfvén waves, having assumed that the energetic-particle pitch- where ν A is specified to be independent of μ to mimic in a simple way the effect of resonance broadening (Lee 2005). Furthermore, ν I1 is the scattering frequency of particles interacting with SMFRs associated with the first moment of Equation (1) and defined in terms of a fractional index β < 1 given by By imposing the ordering f 0 ? f 1 ? f 2 consistent with a nearly isotropic pitch-angle particle distribution, retaining only f 0 -terms on the left-hand side of Equation (17), and considering large scattering frequencies ν A and ν I1 on the right-hand side, a simple solution for f 1 in terms of f 0 can be obtained given by the following: is the effective fractional scattering frequency for particle scattering on both Alfvén waves and SMFRs defined as To find a solution of f 2 in terms of f 0 , we derive the second (n = 2) moment of Equation (1) by multiplying each term by  where the fractional scattering frequency for particle scattering on SMFRs associated with the second moment of Equation (1) is given by By again applying the ordering f 0 ? f 1 ? f 2 , retaining only f 0 -terms on the left-hand side of Equation (21), and by considering large scattering frequencies ν A and ν I on the right-hand side, a simplified solution for f 2 in terms of f 0 can be obtained given by s is the effective scattering frequency for particle scattering on both Alfvén waves and SMFRs defined as Upon inserting the solutions for f 1 (Equation (19)) and f 2 (Equation (23)) into Equation (14), the closed equation for the evolution of f 0 (x, p, t) (the direction-averaged or isotropic part of the energetic-particle distribution 〈f〉(x, p, μ, t)) can be expressed as a fractional Parker-type transport equation in a conservation form according to or in a fractional diffusion-advection form as Alternatively, the closed equation can be cast in the form of a Gleeson-Axford-type fractional transport equation given by is known as the Compton-Getting factor for transforming the differential particle flux between inertial frames (e.g., Gleeson & Axford 1968). In Equations (25)-(27), It is insightful to write Equation (27) compactly as a continuity equation for the differential density n p (x, where the energetic-particle differential flux in ordinary space is and the differential flux in momentum space is In the definitions for S p and J p , 〈〉 Ω indicates an average over all directions in momentum space, and

The Fractional Parker Transport Equation to Be Solved
Let us consider a simplified version of the fractional Parker transport equation (Equation (25)) in one spatial dimension appropriate for studying anomalous DSA at a parallel shock when energetic particles interact on both sides of the shock solely with SMFRs. Interactions with Alfvén waves are neglected (ν A = 0), allowing us to focus on a shock where SMFRs form the dominant low-frequency turbulence component. We do not think that the inclusion of Alfvén waves will change the superdiffusive nature of energetic-particle parallel transport in the vicinity of the shock for the following reasons: (1) Ahead of interplanetary shocks in the inner-heliosphericbackground slow solar wind, the evidence suggest that lowfrequency quasi-2D turbulence (SMFR turbulence) is energywise by far the dominant component, followed by Alfvén-wave turbulence as a minor secondary component (see Section 3.1). In a recent analysis of observations of the transmission of turbulence through three heliospheric shocks at 1, 5, and 84 au with a turbulence transmission model, Zank et al. (2020) concluded that both upstream and downstream of these shocks SMFR turbulence carry far more magnetic and kinetic energy compared to Alfvén-wave turbulence.
(2) Even if Alfvén waves can significantly impact energetic-particle transport in the vicinity of shocks, we still expect parallel transport to be superdiffusive because, as discussed by , the superdiffusion is a product of subdiffusion in pitchangle space arising from the observation that the Alfvén-wave pitch-angle scattering times in the solar wind have a power-law distribution instead of a single value.
The simplified equation is given by where f 0 (z, p, t), and z is the spatial coordinate perpendicular to a planar shock surface. The large-scale solar wind flow U 0 (z), now measured in the shock frame, and magnetic field also points in the direction perpendicular to the shock surface. By multiplying Equation (35) with the left Riemann-Liouville where an escape term with an escape time τ esc (z, p) and a source term Q(z, p) were added on the right-hand side. The first term was transformed to a left Caputo fractional derivative by applying , , having assumed that the initial condition ∂f 0 (z, p, 0)/∂t = 0. Since we are injecting particles at the shock at all times with a constant source term Q(z, p), there is no need for the initial condition.
Note that at later times the time-dependent fractional parallel transport (Equation (36)) simplifies to a time-dependent anomalous diffusion equation for β < 1 given by since the advection, source, and loss terms decay as 1/t 1− β relative to the anomalous diffusion term in Equation ( in those terms. Equivalently, Equation (35) converges to a time-dependent anomalous diffusion equation at late times since the anomalous diffusion term grows as t 1− β relative to the other terms in Equation (35) because of the presence of the fractional time derivative in the anomalous diffusion term. Either way, this suggests that ahead of the shock the competition between the advection of energetic particles toward the shock by the solar wind flow and the anomalous parallel diffusion away from the shock (anomalous modulation) weakens with time. Our solution for time-dependent anomalous parallel DSA presented further below supports this interpretation. Further inspection of Equation (37) reveals that the anomalous diffusion term has normal spatial derivatives (fractional index α = 2), but that the time derivative term includes a fractional time derivative with a fractional index in the range 1 < 2 − β < 2 because of the restriction 0 < β < 1. Thus, the time-dependent anomalous parallel diffusion is a transport process with a Hurst fractional index in the range 1/2 < H = (2 − β)/2 < 1, which puts the transport in the superdiffusive domain (e.g., Sanchez et al. 2006). To put Equation (37) in context, we can think of a diffusion equation in which the time derivative is replaced by a fractional time derivative with a fractional index varying between zero and two. Equations of this type are classified as time-fractional diffusionwave equations where, for a fractional index between zero and one, they are called time-fractional diffusion equations and, for a fractional index between one and two, they are known as timefractional wave equations. Since Equation (37) holds for 1 < 2 − β < 2, it is a time-fractional wave equation that varies between a diffusion equation (2 − β = 1) to a wave equation (2 − β = 2) (e.g., Sanchez et al. 2006;Povstenko 2008;Luchko 2013). Solutions of the time-fractional wave equation equation (Povstenko 2008;Luchko 2013) suggest that in our application an initial particle pulse will spread out fast diffusively when the time-fractional index is close to one so that unphysically some particles in the tail of the pulse will propagate faster than the speed of light. In contrast, the particle pulse will propagate without dissipation at a fixed velocity, thus ensuring that all particles adhere to causal propagation obeying the relativistic constraint when the fractional index reaches two. For a fractional index between one and two, the particle pulse will evolve as a propagating wave that increasingly experiences diffusive spreading with time (wave damping). Evidence of wave-like features in our solutions for parallel superdiffusive shock acceleration are discussed further below.
The coefficient for parallel superdiffusion in Equation (37) where 1 < 2 − β < 2 because of the constraint 0 < β < 1. Thus, 〈|Δz| 2 〉 ∝ (Δt) 2−β , confirming that we are in the superdiffusive regime of parallel anomalous diffusion. Interestingly, the timedependent parallel superdiffusion equation (Equation (37)), which is part of the fractional Parker transport Equation (36), can be traced back to the underlying equation for anomalous pitch-angle diffusion that forms part of the pitch-angledependent diffusion-advection Equation (1). The equation for anomalous pitch-angle diffusion is ) . For this equation, 1 < α < 2 and 0 < β < 1 so that the fractional Hurst index varies in the range 0 < H < 1, thus covering the whole traditional subdiffusive (0 < H < 1/2) and superdiffusive (1/2 < H < 1) regimes. Thus, we see that this equation for general time-dependent anomalous diffusion in μ-space, which includes both fractional time and μderivatives, has transformed in the Parker limit to a timedependent parallel superdiffusive equation with a fractional time but normal spatial derivatives, irrespective of the fractional index α of the μ-derivatives. This narrowing in the range of anomalous transport in the Parker limit appears to be the result of enforcing a nearly isotropic pitch-angle distribution in deriving the fractional Parker transport equation from Equation (1). Nonetheless, the magnitude of k b a -= I 2 , 2 || [ ] is sensitive to the fractional index α of the μ-derivatives of the anomalous diffusion term in μ-space.
A scenario that could potentially support superdiffusive parallel transport for energetic ions interacting with SMFRs in the large-scale solar wind is when (1) these structures can be viewed as coherent, quasi-helical magnetic field structures with a typical width L I of intermediate size and sharp boundaries (Greco et al. 2009), and (2) energetic particles have gyroradii r g that are small compared to the typical SMFR width L I . There is observational support at 1 au for all these assumptions (e.g., Greco et al. 2009;Cartwright & Moldwin 2010;le Roux et al. 2015). Thus, it is reasonable to envisage that the particle guiding centers follow the helical magnetic field lines (negligible cross-field drifts) to propagate through SMFR structures in a coherent fashion for a considerable distance before scattering during exit of the structure through the narrow magnetic boundaries. To be more specific, SMFRs identified at 1 au typically have widths L I ≈ 0.01 − 0.001 au (Cartwright & Moldwin 2010;Khabarova et al. 2015), while data analysis of low-frequency 2D turbulence suggests elongated SMFR structures with an average length L I|| so that L I|| /L I ≈ 2 − 3 (Weygand et al. 2011). Based on the expression for the path length s I of a SMFR helical field line (Equation (9)), we estimate that the distance a particle guiding center travel through a SMFR structure by following magnetic field lines can be considerable because s I ≈ 0.004 − 0.04 au. Furthermore, the protons accelerated by SMFRs reach kinetic energies of 5 MeV (Khabarova & Zank 2017) in the solar wind at 1 au, thus easily fulfilling the requirement r g = L I . Consequently, the transport through SMFRs can be characterized as involving an enhanced probability of relatively long transport distances or step sizes along the background/guide field perpendicular to the shock surface (in the z-direction) between scattering events indicating Lévy walks as discussed by  and Zimbardo et al. (2021).
This begs the question, what kind of anomalous transport in μ-space does this correspond to? One can envisage particles gyrating around the helical magnetic field through the SMFR experiencing little variation in their pitch angles. In our approach, the random changes in μ are reserved for when particles enter or exit the SFMR structure through the narrow SMFR boundary. We estimate that energetic protons will traverse SMFR structures during a time interval in the range Δt ≈ 0.42-4.2 hours for 1 keV protons, or 25-2.5 minutes for 100 keV protons. Thus, the timescale for traversing such large magnetic structures can be considerable at lower suprathermal particle energies so that the waiting times between random changes in μ (SMFR residence times) can be relatively long. We conclude that, in this scenario, the parallel superdiffusion corresponds best to subdiffusive transport in μ-space (see discussions by and Zimbardo et al. 2021 in a different context). To conform to this interpretation, Equation (43) for general time-dependent anomalous pitchangle diffusion needs to be restricted to subdiffusion. This can be accomplished by specifying α = 2 in Equation (43) resulting in normal μ-derivatives. Then, the underlying equation for anomalous diffusion in μ-space that generates the parallel superdiffusion is given by implying that 〈|Δμ| 2 〉 ∝ (Δt) β (0 < β < 1), or that the Hurst fractional index is in the range 0 < H = β/2 < 1/2, whereby anomalous transport in μ-space is restricted to the subdiffusive regime. Consequently, the expression for parallel superdiffusion simplifies to , and τ esc are spatially uniform. Thus, we have to solve The source term Q(z, p) is also neglected because it is reserved for particle injection at the shock. After introducing the Laplace transform to calculate the Laplace transform of Equation (46), we find that the solution of Equation (46) for f z p s , , 0 ( ) ahead of the shock z > 0 (shock is located at z = 0) is  Next, in the equation for parallel transport given by Equation (36), we assume that f 0 (z, p, t) is continuous across the shock and treat the shock as a discontinuity from the perspective of energetic particles. After retaining the dominant transport terms at the shock, and integrating those terms across the shock from z = −ò to z = +ò letting ò → 0, one finds a condition for parallel superdiffusive shock acceleration given by after having specified a steady-state particle source Q(z, p) = S (z, p)δ(z) at the shock for injection purposes into superdiffusive shock acceleration, and having introduced the Compton-Getting factor C given by Equation (28). This condition suggests that the anomalous differential particle flux is conserved across the shock in the absence of particle injection at the shock (S(0, p) = 0). Upon taking the Laplace transform of the acceleration condition (Equation (52)), one can express this condition as   (54) and (55), s sh is the shock compression ratio, and q = 3s sh /(s sh − 1) is the familiar power-law index for the accelerated particle distribution f 0 (0, p) ∝ p − q as predicted by the standard steady-state DSA theory (Drury 1983). Variables with a subscript including "1" are upstream variables whereas variables with a subscript including "2" indicate downstream variables. In deriving Equation (53)   The solution to Equation (57) can be simplified by specifying the following: (i) τ esc1,2 → ∞ to ensure negligible anomalous diffusive particle escape from the shock during acceleration, (ii) , and (iii) g 1 (p, s) = g 2 (p, s), implying that the parallel superdiffusion timescale is equal on both sides of the shock (τ ad1 (p) = τ ad2 (p) = τ ad (p)). Then, the solution for parallel superdiffusive shock acceleration can be expressed as refers to the fractional parallel superdiffusive-shock-acceleration timescale for particles with momentum p (for confirmation, see the next section). The solution to Equation (59) is partly the result of having specified a mono-energetic steady-state particle source at the shock where p 0 is the momentum of the source particles, and n is the total number density of the source particles in the associated thin spherical momentum shell with radius p 0 advected into the shock from upstream at the upstream solar wind speed U 01 as measured in the shock frame. Furthermore, the solution to Equation (59) includes the specification of a general power-law relationship for the momentum-dependence of the parallel superdiffusion coefficient modeled according to The inverse Laplace transform in Equation (59) can be exactly evaluated using the convolution theorem for Laplace transforms so that the solution for parallel superdiffusive shock acceleration at the shock becomes where the time variable t is normalized according to b = t t T p, ADSA ( ) with T ADSA (p, β) determined by Equation (60). The solution is valid for 0 < β < 1. For the special case of superdiffusive parallel transport when β = 1/2, resulting in 〈|Δz| 2 〉 ∝ (Δt) 3/2 or a Hurst fractional index H = 3/4, a much simpler solution can be found. After evaluation the inverse Laplace transform in Equation (59) for β = 1/2, we find that ) , and erfc is the complimentary error function. Although the solution does not apply to early times (t/τ ad = 1), it still produces a sensible result, because then  t erfc 1 2 0 1 2

( ( ))
, and the solution converges to the well-known power-law spectrum result predicted by the standard steady-state DSA theory at a shock (e.g., Drury 1983).
Consider the superdiffusive case where β = 1/3. Then, 〈|Δz| 2 〉 ∝ (Δt) 5/3 or H = 5/6, representing more extreme parallel superdiffusion than β = 1/2. Evaluation of the inverse Laplace transform in Equation (59)  in Equation (65), and we recover again the classical solution for steady-state DSA. Third, when β = 1, |Δz| 2 〉 ∝ (Δt) or H = 1/2, we recover normal parallel diffusion. Upon evaluating the inverse Laplace transform in Equation (59) for this case, we find that ) , and the Heaviside step function -H t 1 ( ) indicates a time-dependent cutoff in the accelerated power-law spectrum that shifts to higher particle momenta with increasing time. When t → ∞ , the step function can be discarded, and the solution to Equation (66) agrees exactly with the standard steady-state DSA solution.
Thus, for both β = 1/2 and β = 1/3, the solutions for timedependent parallel superdiffusive shock acceleration converge to the standard steady-state DSA solution at late times. This indicates that our fractional Parker transport theory predicts that time-dependent parallel superdiffusive shock acceleration, even though involving anomalous diffusion operating on different different timescales than normal diffusion in classical timedependent parallel DSA, will eventually reach the same familiar steady-state accelerated spectrum µ - ) at the shock as predicted by the standard steady-state DSA theory. We suspect that this outcome is the result of having assumed in the derivation of the fractional Parker transport equation above that the anomalous pitch-angle diffusion allows for enough particle scattering on large spatial scales to maintain a nearly isotropic pitch-angle distribution. Inspection of simulations of anomalous transport of individual particles indeed shows that, despite the presence of long waiting times or step sizes between scattering events, there are also many occurrences of short waiting times and step sizes, which might be sufficient to maintain a nearisotropic distribution (Effenberger 2012;Zimbardo et al. 2017).

The Timescale for Parallel Superdiffusive Shock Acceleration
In this section, we derive an expression for the timescale over which parallel superdiffusive shock acceleration generates particles of momentum p to determine how this timescale features in our solutions for shock acceleration presented above. Assuming a balance between fractional advection and parallel superdiffusion away from the shock in the limit of long shock escape times (τ esc → ∞ ), and applying dimensional analysis, the parallel fractional Parker transport Equation (46)  This yields an estimate of the effective parallel superdiffusion spatial scale L ad over which the particle distribution decays away from the shock given by where τ ad is the effective parallel superdiffusion timescale, and 0 < β < 1. The parallel superdiffusion timescale for the upstream particles crossing the shock from downstream is Solving for τ ad1 , we get an expression for the parallel superdiffusion timescale to a fractional power β < 1 given by In the same way, the parallel superdiffusion timescale to the fractional power β for the downstream particles crossing the shock from upstream is Therefore, the total parallel superdiffusion timescale to a fractional power β for the particles to complete one acceleration cycle back and forth across the shock is . 7 2 The characteristic parallel superdiffusive-shock-acceleration timescale to the fractional power β, based on the particlemomentum increase during a shock-crossing cycle, is defined as where q = 3s sh /(s sh − 1) as defined above, and s sh is the shock compression ratio. We assume that this expression can be generalized to yield the time for accelerating particle from an initial momentum p 0 to a final momentum p according to In the special limit of normal diffusion (β = 1), Equation (74) yields the the correct acceleration time for normal time-dependent DSA (Drury 1983). To relate b T ADSA to the expression from Equation (60) for b T ADSA in our solutions for the accelerated spectrum at the shock, we assume   (60) as the parallel superdiffusive-shock-acceleration time.
We compare the efficiency of parallel superdiffusive shock acceleration with normal DSA. The ratio of the acceleration times for superdiffusive shock acceleration over classical DSA is , li is the logarithm integral function, and v 0 is the initial particle speed. Furthermore, ò = s I /τ c where p = + s L L L 1 || || is the path length that a particle guiding center follows along a helical SFMR field line, L I|| is the mean length of SMFRs in the axial or guide field direction, and L I is the average SMFR width. Also, in ò, in which, assuming that the anomalous pitch-angle scattering of energetic particles is dominated by abrupt magnetic mirroring/focusing when particles cross SMFR boundaries (see Equations (5)-(8)), we can approximate the anomalous diffusion coefficient for pitch-angle scattering mm b a= ( | | ) is the mean effective particle correlation time modeled as a competition between the energetic-particle SMFR inverse crossing time and the SMFR inverse turbulence correlation time.
Upon inserting the expression for mm b a= , we get in the fast particle limit (v/V A0 ? 1) that  (76) is a || = 1. The normal paralleldiffusion coefficient has a slightly weaker momentum-dependence because of the natural logarithmic function. To estimate the ratio of the parallel shock-acceleration times T ADSA /T DSA using the expression from Equation (76), we specify the following parameters that include the reasonable parameters related to traveling shocks at 1 au: β = 1/2 so that the fractional index for superdiffusion is 2 − β = 3/2, p/p 0 = 10, the shock compression ratio s sh = 2 yielding a steady-state spectral power-law index q = 6, the upstream flow speed in the shock frame U 0 = U 01 = 200 km s −1 , the source particle injection  Smith et al. 2006). We find that T ADSA /T DSA ≈ 39, suggesting that the parallel superdiffusive shock acceleration is significantly less efficient than normal parallel DSA at a traveling shock in the solar wind at 1 au. We note that Zimbardo & Perri (2013) and Perri & Zimbardo (2015) came to the opposite conclusion, based on taking into account that the particle step-size probability distribution is a power law (Lévy distribution) only beyond a certain minimum step size and finding that the minimum step size is small at heliospheric shocks. This is an issue that needs further study.
having applied the late time condition is the characteristic fractional superdiffusive timescale upstream) and the no escape limit τ esc1 → ∞ . Further insight can be gained by realizing that k containing inverse Laplace transforms that can be evaluated analytically. After the evaluation of the inverse Laplace transforms, the upstream solution becomes ). The solution is valid for 0 < β < 1.
. By pulling out the erfc-function as a common factor, the solution can be more compactly reexpressed as where E 1/2 is a Mittag-Leffler function valid for z L p t , 1 ad1 ( )  , and the effective superdiffusion spatial scale upstream of the shock is defined as The function g(x) was introduced as a consequence of approximating the erfc-function as which is especially accurate for x 0 (see the website https:// mathworld.wolfram.com/Erfc.html).
In the limit  ¥ t ,  t erfc 1 2 1 1 2 ( ( )) . Furthermore, x(p, t) → 0 so that g(x) → 1, and  ¥ L ad1 whereby E 1/2 ( − z/L ad1 (p, t)) → 1. This suggests that the accelerated particle distribution, even though it converges to the standard solution for steady-state DSA at the shock as discussed above, behaves differently from the standard theory upstream of the shock by becoming spatially uniform, instead of decaying exponentially with increasing distance from the shock. This can be understood by applying dimensional analysis to the equation for parallel fractional transport (Equation (46)), which can then be expressed as where we neglected the escape term by assuming a long escape time consistent with our upstream solution. In Equation (89), L is the spatial scale over which the particle distribution decays ahead of the shock, and T is timescale over which the distribution function evolves upstream. When T → ∞ in this equation, the equation simplifies to a steady-state equation which implies that L → ∞ assuming a finite value for f 0 (the steady-state particle source at shock ensures that f 0 stays finite in the vicinity of the shock). This indicates that, as the energetic-particle advection term in Equation (89) decays with time relative to parallel superdiffusion, the advection of energetic particles toward the shock becomes increasingly inefficient compared to parallel superdiffusion away from the shock. As a result, particles at late times are freely undergoing parallel superdiffusion away from the shock without any modulation from the incoming solar wind flow, thus forming a plateau distribution function upstream. It is unlikely that such a uniform upstream particle distribution will occur. There are indications that in complex systems, anomalous diffusion processes do cross over into normal diffusion states for system evolution times longer than a certain correlation time (e.g., Ito & Miyazaki 2003). In the solar wind, a crossover from anomalous diffusion to normal diffusion is also expected to occur when energetic particles temporarily trapped in SMFR structures escape and propagate along diffusing magnetic field lines (Ruffolo et al. 2003), or when solar energetic particles escaping nearly scatter free in front of traveling shocks scatter on magnetic turbulence at later times (e.g., Effenberger & Litvinenko 2014).

The Downstream Solution for Parallel Superdiffusive Shock Acceleration
According to Equation (50), we can express the downstream solution for parallel superdiffusive shock acceleration as where λ 2 (p, s) is defined by Equation (51). By taking the Laplace transform of the solution for f (0, p, t) in Equation (59) and inserting this result for f p s 0, , 0 ( ) in the downstream solution, we get having applied the late time condition t/τ ad2 ? 1 (τ ad2 is the characteristic superdiffusive timescale downstream) and the no escape limit τ esc2 → ∞ . After assuming for simplicity that |z|/(|U 02 |t) = 1 so that » + z U s z U s ). To gain more insight, we again specify β = 1/2 to simplify the general solution. Evaluation of the inverse Laplace transforms in Equation (93) for β = 1/2 yields the solution adv2 adv2 adv2 02 and x(p, t) and g(x(p, t)) are specified in Equation (87). If  ¥ t , x(p, t) → 0, g(p, t) → 1, L adv2 → ∞ , so that E 1/2 → 1. Accordingly, we recover at late times the result of the familiar accelerated particle power-law spectrum f 0 (p) ∝ p − q that is uniform, spatially behind the shock as predicted by the standard steady-state DSA theory. Note that the evolution toward a plateau downstream is controlled not only by the increase in the time-dependent advection scale length L adv2 = |U 02 |t but also by the effective, fractional parallel-diffusion-length scale downstream through the parameter x(p, t) because the latter depends on the parallel superdiffusive-shock-acceleration timescale T ADSA (p, β = 1/2) (see the first line below Equations (85) and (87)). Thus, for t/T ADSA ? 1, µ L t adv2 3 2 , indicating that the formation of the downstream plateau for f 0 appears to develop more rapidly compared to upstream where the characteristic spatial scale µ L t ad1 1 2 .

More General Solutions for Parallel Superdiffusive Shock Acceleration
So far, we concentrated in producing the analytical solutions for parallel superdiffusive shock acceleration at the shock as well as upstream and downstream. The results were limited to having no particle escape from the shock and tractable solutions were presented for mainly one fractional index value for parallel superdiffusion. To facilitate more general solutions that include the particle escape from the parallel shock due anomalous perpendicular diffusion, we retain particle escape effects in the solution for the accelerated particle spectrum at the shock as presented in Equations (57) and (58) by not simplifying the functions g 1,2 (p, s) in Equation (58). Then, as before, we assume times later than the effective parallel superdiffusive timescale so that the functions g 1,2 (p, s) = 1 in Equation (58) for expansion purposes, and for simplicity specify that g 1 (p, s) = g 2 (p, s) = g(p, s), implying that both the effective parallel superdiffusive timescale t = and behind the shock reads where the acceleration time for parallel superdiffusive shock acceleration unaffected by energetic-particle escape (T ADSA ) is defined as before by In Equations (98)-(100), the particle escape time τ esc is modeled assuming that energetic particles can escape a parallel shock with a finite size by transport in the perpendicular direction. It assumed that, because particles tend to be trapped SMFRs, perpendicular particle escape occurs as an anomalous diffusion process that is time fractional, as is discussed further below. Accordingly, the characteristic escape time is defined as a fractional time with an expression given by where L sh is the size of the shock related to its curvature radius, is the μ-averaged diffusion coefficient for perpendicular anomalous diffusion of energetic particles interacting with SMFRs (see definition in Equations (5)- (7)). Also, α is the fractional index of the spatial derivatives, and β ⊥ is the fractional index of the time derivative of time-dependent anomalous perpendicular diffusion (β ⊥ should be distinguished from β, which refers to the time-fractional index of anomalous parallel diffusion). In this simplified model of energetic-particle escape from a parallel shock, the particles escape when they have diffused anomalously a distance of the order of the shock size L sh across the magnetic field (e.g., Forman & Webb 1985).
in the fast (super Alfvénic) particle limit (V A0 /v = 1). Thus, can be thought of as modeling the anomalous perpendicular diffusion of energetic particles as a consequence of these particles following the helical magnetic field lines of numerous coherent SMFR structures (temporary trapping). In Equation (104), 0 < α < 2, 0 < β ⊥ < 1, and b a = C , , the power-law index a ⊥ of the momentum-dependence of k á ñ b a m^p I , ( ) [ ] listed in the escape time expression (Equation (103)) is a ⊥ = β ⊥ .
In the results, which we present below, we specify that α = 2 (C(β ⊥ , α = 2) = Γ(β ⊥ )/2). Thus, in the 2D plane perpendicular to the background/guide field, the particle displacement is statistically characterized by á D ñ = á D ñ µ D bx y t 2 2 | | | | ( ) , or the Hurst index 0 < H < 1/2. This indicates that we are restricting anomalous perpendicular diffusion for energetic particles interacting with numerous SMFRs to the subdiffusive regime because 0 < β ⊥ < 1. This is consistent with our discussion above (see discussion below Equation (43) in Section 5.1) where we envisaged the energetic particles largely following the quasi-helical magnetic field lines of coherent SMFR structures before scattering when crossing the sharp SMFR boundaries as a way to explain parallel superdiffusion. Viewed in the 2D plane perpendicular to the background/guide field, the particle guiding trajectories can be thought of as approximately circular, indicating little spatial displacement in the perpendicular direction when energetic particles traverse SMFR structures. Random perpendicular spatial displacement is reserved for the particle entry into and the exit from SMFRs. The waiting times between random displacements (SMFR residence times) of lower-energy suprathermal energetic particles in SMFRs were estimated above in the solar wind at 1 au to be potentially considerable for protons. On this basis, we think it is reasonable to model anomalous perpendicular diffusion through SMFRs as a subdiffusive process.
Having evaluated the inverse Laplace transforms in the solutions to Equations (98)-(100) numerically with MATLAB, we present time-dependent solutions showing how superdiffusive shock acceleration of energetic particles at a parallel shock is modified as a function of the level of parallel superdiffusion (2 − β-value) and the level of perpendicular subdiffusive particle escape from the parallel shock (β ⊥ -value). Information about the parameters that we specified for these calculations can be found in Section 5.3 in the discussion below Equation (80). In Figure 1, we present solutions for the accelerated, direction-averaged particle distribution function f 0 (p) at the shock located at z = 0 as a function of particle momentum normalized to the particle source momentum p 0 . The distribution function curves are normalized to a value of 1 at the injection momentum (p/p 0 = 1). The solutions represent the results calculated for ∼9 days of shock acceleration for different fractional index values (β-values) associated with parallel superdiffusion characterized statistically by 〈|Δz| 2 〉 ∝ Δt 2− β . Particle escape from the shock is ignored in Figure 1. We display solutions for β = 0.9 when 〈|Δz| 2 〉 ∝ Δt 1.1 (black), β = 0.7 when 〈|Δz| 2 〉 ∝ Δt 1.3 (red), β = 0.5 when 〈|Δz| 2 〉 ∝ Δt 1.5 (green), β = 0.3 when 〈|Δz| 2 〉 ∝ Δt 1.7 (blue), and β = 0.1 when 〈|Δz| 2 〉 ∝ Δt 1.9 (purple). The curve with open circles indicates the shock acceleration when parallel diffusion is normal. The results, calculated for a shock compression ratio of 2 (s sh = 2), show that the more superdiffusive parallel transport across the shock is, the less efficient the shock acceleration becomes. This is apparent in the shift of the spectral rollover to lower momenta and the steepening of the spectra. As expected, when parallel superdiffusion is approaching normal diffusion (β → 1), the spectral slope for time-dependent superdiffusive shock acceleration (black curve) converges to the spectral slope for normal time-dependent DSA (curve with open circles). Since shock acceleration is more efficient when β → 1, it so happened that at lower particle energies the black and open-circle curves both reproduced the power-law spectrum f 0 (p) ∝ p −6 predicted by normal steadystate DSA theory within nine days of acceleration. The fact that the black curve, which still represents parallel superdiffusion, reproduces the spectral prediction for normal steady-state DSA theory, and it is consistent with the results discussed above that the steady-state limit for superdiffusive time-dependent shock acceleration reproduces the familiar result of normal steadystate DSA theory, namely that µ f p p s s ) . The abrupt cutoff in the accelerated energetic-particle spectrum of normal time-dependent DSA at a parallel shock (curve with open circles) is in agreement with the analytical solution for regular time-dependent DSA presented above (see Equation (66) and its discussion). This suggests a limitation in the solution when it comes to modeling a more realistic gradual spectral rollover for normal time-dependent DSA.
In Figure 2, we illustrate how the time-dependent solutions of the accelerated particle distribution f 0 (z) at a specific particle momentum (p/p 0 = 3) vary as a function of distance z perpendicular to the parallel shock surface in astronomical units (au) for different levels of parallel superdiffusion after ∼9 days of shock acceleration. The shock is located at z = 0, z > 0 indicates upstream of the shock, and z < 0 means downstream. The parameters specified to calculate results of this figure and in Figure 1 are identical, and just like in Figure 1, the particle escape effects are neglected. As in Figure 1, the solid curves represent β = 0.9 when 〈|Δz| 2 〉 ∝ Δt 1.1 (black), β = 0.7 when 〈|Δz| 2 〉 ∝ Δt 1.3 (red), β = 0.5 when 〈|Δz| 2 〉 ∝ Δt 1.5 (green), β = 0.3 when 〈|Δz| 2 〉 ∝ Δt 1.7 (blue), and β = 0.1 when 〈|Δz| 2 〉 ∝ Δt 1.9 (purple). The results show that the characteristic scale over which f 0 (z) decay ahead of the shock increases the more superdiffusive that the parallel transport across the shock becomes. The curve associated with maximum parallel superdiffusion (purple) suggests that an increasingly exponential decay is being established with time followed by a rollover at higher particle momenta. In the less superdiffusive cases, it appears that more time is needed before those curves will also become exponential at lower particle momenta. Thus, more superdiffusive parallel transport ahead of the shock is more efficient in transporting energetic particles against the incoming solar wind flow away from the shock. The results behind the shock are more complicated. The scale over which the distribution function decays downstream first decreases and then increases with an increasing level of parallel superdiffusion (increasing 2 − β-values). Overall, our impression is that the characteristic distance of effective parallel superdiffusion back and forth across the shock is increasing with increasing superdiffusiveness. This is consistent with the trend in Figure 1 that superdiffusive shock acceleration at a parallel shock becomes less efficient with increasing superdiffusiveness over a finite shock-acceleration time interval, because the diffusive distance that the particles travel back and forth across the shock to complete one acceleration cycle increases.
It is also interesting to see how the solutions ahead of the shock in Figure 2 vary from a diffusive solution to a wave solution as the fractional index 2 − β changes from close to one to close to two, as discussed above. The black curve implies almost purely a diffusive spreading of particles (noncasual propagation not subject to the relativistic constraint) because the spatial decay of the distribution is close to an exponential ahead of the shock (close to a straight line). The purple curve at the bottom, on the other hand, indicates a solution upstream that approximates a wave solution. This can be seen in the prominent rollover in the distribution function that shows that the leading particle pulse ahead of the shock exhibits a finite propagation speed, thus restoring causal propagation obeying the relativistic constraint to the theory. Careful inspection of the rest of the upstream curves reveals how the upstream solutions gradually evolve away from a purely diffusive solution to solutions that combine diffusive and wave characteristics (increasingly strong rollover) as the fractional index 2 − β varies from ∼1 to ∼2. Viewing the different upstream curves from the perspective of a decreasing fractional index, one sees how the leading wave-like particle pulse increasingly gets damped out as the diffusive aspect of the solution becomes more prominent.
The results that we present in Figures 3 and 4 illustrate how varying the level of subdiffusiveness of perpendicular subdiffusive particle escape from the parallel shock affects the solutions of parallel superdiffusive shock acceleration for a given level of parallel superdiffusion. In Figure 3, all the solutions shown for the accelerated particle spectra at the shock are for the case of parallel superdiffusion across the shock when β = 0.5 (〈|Δz| 2 〉 ∝ Δt 1.5 ). What is varied in these figures is the fractional index β ⊥ of perpendicular subdiffusive energeticparticle escape from the shock that determines the perpendicular average displacement á D ñ µ D bx t 2 | | . As before, the shock compression ratio is 2, and the results are shown for ∼9 days of shock acceleration calculated with the same parameters as in the previous figures. We present results for no escape (dashed curve) as a reference, for escape determined by β ⊥ = 0.9 when 〈|Δx| 2 〉 ∝ Δt 0.9 (solid black curve), β ⊥ = 0.7 when 〈|Δx| 2 〉 ∝ Δt 0.7 (red curve), and β ⊥ = 0.5 when 〈|Δx| 2 〉 ∝ Δt 0.5 (green curve). The main effect is that the more subdiffusive perpendicular particle escape from the shock is, the less the accelerated particle spectrum is modified (steepened) by the escape. Sensibly, this shows that the particle escape from the shock becomes less efficient because the mean perpendicular displacement squared of energetic particles across the background magnetic field is reduced.
In Figure 4, we display the solutions showing how varying the level of subdiffusiveness of perpendicular subdiffusive particle escape from the parallel shock affects spatial profiles of the energetic-particle distribution f 0 (z) across the shock for a given level of parallel superdiffusion, i.e., for β = 0.5 when 〈|Δz| 2 〉 ∝ Δt 1.5 ). As in Figure 2, the results are shown for p/p 0 = 3, and for a shock compression ratio of 2 after ∼9 days of shock acceleration. As in Figure 3, the dashed curve forms the reference solution when there is no particle escape, the solid black curve includes escape effects determined by β ⊥ = 0.9 when 〈|Δx| 2 〉 ∝ Δt 0.9 , the red curve indicates the effect of escape for β ⊥ = 0.7 when 〈|Δx| 2 〉 ∝ Δt 0.7 , and the green curve represents the effect of escape for β ⊥ = 0.5 when 〈|Δx| 2 〉 ∝ Δt 0.5 . All the curves were normalized to a value of 1 at the shock position (z = 0). The results show that, ahead of the shock (z > 0), the spatial scale over which f 0 (z) decays increases with an increase in the level of subdiffusiveness in the perpendicular particle escape term (decrease in the β ⊥ -value). This is expected, because the perpendicular subdiffusive transport becomes less efficient so that energetic-particle escape from the shock is reduced. In contrast, behind the shock (z < 0), we found no significant effect from varying the level of subdiffusiveness in the perpendicular particle escape.

Summary and Discussion
Previously, we presented a pitch-angle-dependent fractional diffusion-advection equation for energetic-particle interaction with SMFRs in the solar wind on large spatial scales (MHD scales) derived from first principles from the standard focused transport equation (le Roux & Zank 2021). Since fractional transport, despite the enhanced probability for long waiting times and large step sizes between scattering events, still includes a high probability for small step sizes and short waiting times (Effenberger 2012;Zimbardo et al. 2017), we assumed that it was reasonable to specify a near-isotropic pitchangle distribution for energetic particles. This enabled us to derive a new fractional Parker transport equation that can also be cast in the form of a fractional Gleeson-Axford transport equation, thus opening up the opportunity of finding the analytical solutions for anomalous transport problems that can serve as a bench mark for future numerical solutions. As a first application, we decided to investigate the time-dependent anomalous DSA at a parallel shock.
Simplified for parallel transport, our fractional Parker transport theory for energetic-particle interaction with SMFRs has the following main features: (1) Anomalous parallel diffusion of energetic particles is superdiffusive when the time-fractional index β is restricted to 0 < β < 1. This fractional index β has its origin in the fractional time derivative that appears in the underlying pitch-angle-dependent fractional diffusion-advection equation from which the fractional Parker transport equation was derived above. Time-dependent anomalous parallel diffusion in the fractional Parker transport equation is described in terms of a fractional time derivative with fractional index 2 − β, and a term for parallel anomalous diffusion with normal spatial derivatives (fractional index α = 2). Accordingly, the anomalous parallel diffusion in the Parker limit is predicted to be superdiffusive for β < 1 (〈|Δz| 2 〉 ∝ (Δt) 2−β with 1 < 2 − β < 2).
(2) Energetic-particle advection becomes negligible compared to parallel superdiffusion at late times. Transport terms in the fractional Parker transport equation for advecting energetic particles with the solar wind flow, and for generating adiabatic energy changes of energetic particles induced by the divergence of the solar wind flow, are preceded by a fractional time derivative with a fractional index 1 − β, or equivalently, one can view the parallel superdiffusion term as being preceded by a fractional time derivative with a fractional index β − 1. Thus, the advection terms decay with time as t −(1− β) relative to the term for parallel superdiffusion, or from an equivalent perspective, the parallel superdiffusion term grows with time as t 1− β relative to the advection terms. From either perspective, the outcome at late times is that parallel superdiffusion is the dominating transport process.
Superdiffusive parallel transport of energetic particles through coherent SMFR structures has some physical justification. A scenario with some observational justification that Figure 1. Solutions for time-dependent superdiffusive shock acceleration at a parallel shock with a compression ratio of 2 after ∼9 days of acceleration. Energeticparticle escape from the shock is neglected. Plotted is the direction-averaged particle distribution f 0 (p) as a function of particle momentum p normalized to the source particle momentum p 0 . The curves are normalized to a value of 1 at p/p 0 = 1. Results are shown for parallel superdiffusion when 〈|Δz| 2 〉 ∝ Δt 1.1 (black solid curve), 〈|Δz| 2 〉 ∝ Δt 1.3 (red curve), 〈|Δz| 2 〉 ∝ Δt 1.5 (green curve), 〈|Δz| 2 〉 ∝ Δt 1.7 (blue curve), 〈|Δz| 2 〉 ∝ Δt 1.9 (purple curve), and 〈|Δz| 2 〉 ∝ Δt (open circles). would support such transport is the fact that energetic particles interact with SMFRs characterized by (1) a width L I ≈ l c  0.01 au in the solar wind at 1 au as indicated, e.g., by Cartwright & Moldwin (2010) and Khabarova et al. (2015) (l c is the solar wind turbulence correlation scale) so that their gyroradii r g = L I , (2) coherent (relatively smooth) quasi-helical magnetic fields (Cartwright & Moldwin 2010) with sharp magnetic boundaries (Greco et al. 2009), and (3) small electric fields induced by sub-Alfvénic flows (e.g., . Then it is reasonable to think of particle propagation occurring through SMFRs approximately as undisturbed gyromotion along the quasi-helical field lines of SMFRs with limited variation in pitch angle. Consequently, the parallel transport (transport along the SMFR guide field) will be coherent for relatively long distances before scattering occurs when the particles exit these structures through narrow SMFR boundaries. This suggests the large parallel spatial step sizes of up to ∼0.01 au at 1 au between scattering events as expected for superdiffusive transport.
Subdiffusive transport in μ-space best explains our result for superdiffusive parallel spatial transport through SMFRs. The equation for time-dependent anomalous transport in μ-space (see Equation (43)), which forms part of the complete pitchangle-dependent diffusion-advection (Equation (1)), underlies the time-dependent parallel superdiffusion equation that dominates in the fractional Parker transport equation at late times, allowing for subdiffusion or superdiffusion in μ-space. However, given the proposed characteristics of parallel superdiffusive transport through SMFRs as discussed above, the parallel superdiffusion would best correspond to particle transport in pitch-angle space inside SMFRs being subject to relatively long waiting times (SMFR residence times) before scattering in pitch-angle occurs upon SMFR exit. This implies that the particle transport in μ-space that generates parallel superdiffusion is a subdiffusive process. We estimated that the waiting times can be a number of hours for lower-energy suprathermal protons at 1 au. Accordingly, we proposed that the correct equation for time-dependent anomalous diffusion in μ-space that causes parallel superdiffusion is one where the time derivative is fractional with index β, but the term for anomalous pitch-angle diffusion contains normal pitch-angle derivatives (fractional index α = 2) resulting in 〈|Δμ| 2 〉 ∝ (Δt) β . Given the constraint 0 < β < 1, the transport is subdiffusive. We note that a similar connection between subdiffusive pitch-angle transport and superdiffusive parallel transport from a different perspective was made in the case of particle interaction with Alfvén waves propagating in a nonuniform, random plasma medium .
Theoretically, the superdiffusive parallel transport described by the fractional Parker equation is facilitated by the fact that our underlying pitch-angle-dependent fractional diffusion-advection equation is nonlinear, meaning the latter model particle guiding Figure 2. Solutions for time-dependent superdiffusive shock acceleration of energetic particles with p/p 0 = 3 across a parallel shock with a compression ratio of 2 after ∼9 days of acceleration. Particle escape is neglected. Plotted is the normalized direction-averaged particle distribution f 0 (z) as a function of distance z in astronomical units (au) perpendicular to the shock surface. The shock is located at z = 0, z > 0 means ahead of the shock, while z < 0 implies behind the shock. Just as in Figure 1, results are shown for parallel superdiffusion when 〈|Δz| 2 〉 ∝ Δt 1.1 (black solid curve), 〈|Δz| 2 〉 ∝ Δt 1.3 (red curve), 〈|Δz| 2 〉 ∝ Δt 1.5 (green curve), 〈|Δz| 2 〉 ∝ Δt 1.7 (blue curve), and 〈|Δz| 2 〉 ∝ Δt 1.9 (purple curve). center trajectories are coherently disturbed when propagating through intermediate-size coherent SMFR structures at 1 au in the solar wind before scattering in pitch-angle space when exiting SMFR boundaries. Our estimate that the guiding center trajectories can be coherently disturbed on relatively long timescales inside such intermediate-scale SMFR structures before they exit these structures (relatively long waiting times between scattering events) is captured statistically by the specification of a Lévy distribution for waiting times in the fractional diffusion-advection equation that puts additional weight on long waiting times (le Roux & Zank 2021). In our simple model of the disturbed-particle guiding center trajectory following a helical SMFR field line through the SMFR structure, the particles were estimated, on average, to experience a modest time delay due to the increased path length associated with the helical field line. This suggests the linking of parallel superdiffusion generated by subdiffusive pitch-angle transport to weak particle trapping in SMFRs. Our estimates based on reasonable parameters at 1 au suggest that the majority of particles in the pitch-angle distribution, ∼80%, experience weak trapping while the rest might undergo stronger trapping in the form of magnetic mirroring between the strongest magnetic fields inside SMFR structures (Guidoni et al. 2016).
The key results of our study of time-dependent superdiffusive shock acceleration with our fractional Parker transport theory at a parallel shock, as found from analytical solutions, are as follows: (1) The timescale for superdiffusive shock acceleration of energetic particles to reach a certain momentum is fractional with an index 0 < β < 1 and depends on the value of the parallel superdiffusion coefficient. This is different from superdiffusive shock acceleration where particle transport near the shock is modeled with a normal time derivative and fractional spatial derivatives (Lévy flights) where the timescale for superdiffusive shock acceleration is not fractional (Zimbardo et al. 2017).
(2) The superdiffusive-shock-acceleration time was estimated to be considerably longer compared to the normal DSA time for β = 0.5 (〈|Δz| 2 〉 ∝ Δ 2−β=3/2 ) at parallel shocks. This result is in disagreement with earlier work that suggests superdiffusive shock acceleration to be more efficient (Zimbardo & Perri 2013;Perri & Zimbardo 2015). These authors based there conclusion on taking into account that the particle step-size probability distribution is a power law (Lévy distribution) only beyond a certain minimum step size and that the minimum step size is small at heliospheric shocks. This issue needs to be further investigated in future work. (3) At late times, the accelerated energetic-particle spectrum at the shock converged to a steady state that agrees with the result from classical steady-state DSA theory. We think that this outcome followed from the assumption that subdiffusive pitch-angle scattering allowed for enough particle scattering on large spatial scales to maintain a nearly isotropic particle distribution, as discussed above. Our result is in contrast to the earlier work that suggested that superdiffusive-shock-acceleration results in a Figure 3. Solutions for time-dependent superdiffusive shock acceleration at a parallel shock with a compression ratio of 2 after ∼9 days of acceleration. Energeticparticle escape from the shock by subdiffusive perpendicular diffusion is included. Plotted is the direction-averaged particle distribution f 0 (p) as a function of particle momentum p normalized to the source particle momentum p 0 . The curves for f 0 are normalized to a value of 1 at p/p 0 = 1. For all solutions shown, parallel superdiffusion is specified as 〈|Δz| 2 〉 ∝ Δt 1.5 . Solutions are presented for perpendicular subdiffusive shock escape where 〈|Δx| 2 〉 ∝ Δt 0.9 (sold black curve), 〈|Δx| 2 〉 ∝ Δt 0.7 (red curve), and 〈|Δx| 2 〉 ∝ Δt 0.5 (green curve). The dashed curve is a reference solution without particle escape. particle spectrum that is harder than predicted by classical DSA (e.g., Zimbardo & Perri 2013). (4) At late times, the accelerated particle distribution converges to a steady-state uniform distribution ahead of the parallel shock in contrast to the prediction of an exponential spatial decay in the upstream distribution according to the standard steady-state DSA theory. This unexpected result is a direct consequence of the spatial advection of energetic particles toward the shock becoming increasingly inefficient compared to parallel superdiffusion away from the shock. As a result, the particles at late times are freely undergoing parallel superdiffusion away from the shock without significant modulation from the incoming solar wind flow. It is unlikely that such a uniform upstream particle distribution will ever form though, because there are indications that in complex systems anomalous diffusion processes do cross over into normal diffusion states for system evolution times longer than a certain correlation time (e.g., Ito & Miyazaki 2003). In the solar wind, a crossover from anomalous diffusion to normal diffusion is also expected to occur when energetic particles temporarily trapped in SMFR structures escape and propagate along diffusing magnetic field lines (Ruffolo et al. 2003), or when solar energetic particles, escaping nearly scatter free in front of traveling shocks, scatter on magnetic turbulence at later times (e.g., Effenberger & Litvinenko 2014). (5) The shockaccelerated particle distribution behind the shock converges to a uniform spatial distribution at late times as predicted by the standard steady-state DSA theory. This evolution toward a plateau is controlled not only by the increase in the advection distance of energetic particles with time but also by the effective, fractional parallel-diffusion-length scale downstream. The formation of the downstream plateau for particle distribution is predicted to develop more rapidly compared to upstream.
Additionally, more flexible semianalytical solutions revealed for a given acceleration time are as follows: (1) The more superdiffusive the parallel transport across the shock is, the less efficient the shock acceleration becomes. This was manifested in accelerated spectra at the shock that becomes steeper with a spectral rollover at smaller particle momenta, and increased characteristic spatial scales of decay of the distribution function ahead of the shock as a function of an increasing level of superdiffusiveness. Thus, the increasingly inefficient superdiffusive shock acceleration appears to be linked to the increasing characteristic distances of parallel superdiffusion, back and forth across the shock. (2) Increasing the superdiffusiveness of parallel transport results in the accelerated particle distribution upstream, exhibiting to a higher degree the characteristics of a wave solution and less the features of a diffusive solution. This can be seen as the energetic-particle distribution upstream deviates increasingly from a purely exponential spatial decay to a spatial decay with a prominent rollover. This happens as the leading particle pulse ahead of the shock increasingly exhibits a finite . Solutions for time-dependent superdiffusive shock acceleration of energetic particles with p/p 0 = 3 at a parallel shock with a compression ratio of 2 after ∼9 days of acceleration. Energetic-particle escape from the shock by subdiffusive perpendicular diffusion is included. Plotted is the direction-averaged particle distribution f 0 (z) as a function of distance z in astronomical units (au) relative to the shock. The curves normalized to a value of 1 at the shock located at z = 0. z > 0 indicates distance upstream of the shock whereas z < 0 represents distance downstream. For all solutions shown, parallel superdiffusion is specified to be 〈|Δz| 2 〉 ∝ Δt 1.5 . Just like in Figure 3, solutions are displayed for perpendicular subdiffusive shock escape where 〈|Δx| 2 〉 ∝ Δt 0.9 (sold black curve), 〈|Δx| 2 〉 ∝ Δt 0.7 (red curve), and 〈|Δx| 2 〉 ∝ Δt 0.5 (green curve). The dashed curve is a reference solution without particle escape.