Recent Developments in Particle Acceleration at Shocks: Theory and Observations

Energetic particles are ubiquitous in space and astrophysical environments. In the interplanetary space, energetic particles populating the suprathermail tail of the low-energy cosmic ray spectrum range from solar energetic particles (SEPs) (Lin 1974, 2005) to pick-up ions, to corotating interaction region particles, energetic storm particles, impulsive-solar-flare-related particles, to anomalous cosmic rays observed close to the termination shock (TS) of the solar wind (Giacalone et al. 2012) (see Lee et al. (2012), Mewaldt (2000) for a general overview), and to galactic cosmic rays (Potgieter 1998). Understanding the different acceleration mechanisms represents an open challenge for modern astrophysics. It has recently been discovered that the presence of coherent structures, such as current sheets, flux tubes, or multiple reconnecting islands, represents an efficient stochastic acceleration mechanism for charged particles accounting for the generation of power law tails of suprathermal particles (Drake et al. 2006; Hoshino 2012; Vlahos et al. 2004, 2016). Actually, the presence of vortical structures and reconnecting magnetic islands can also play a fundamental role during particle shock acceleration, as pointed out by Zank et al. (2015a,b) and, e.g., le Roux and Zank (2021). Indeed, considering that the plasma downstream of a shock wave is highly compressed and turbulent, it is natural to consider that a sea of fluctuations, flux ropes, and magnetic islands emerge causing an efficient pitch-angle scattering of particles and a local acceleration in response to island contraction and electric fields from merging structures. Including such effects in the particle transport equation (Zank et al. 2014; le Roux et al. 2015), Zank et al. (2015b) found that the peak of the energetic particle intensity is located downstream behind the shock front; this was observed to be in qualitative agreement with Voyager 2 observations of energetic particle fluxes at the TS of the solar wind (Decker et al. 2008).

Major sources of suprathermal particles in the heliosphere are interplanetary shocks, i.e., driven by coronal-mass-ejections while expanding into the interplanetary medium (Vourlidas et al. 2003; Cliver and Ling 2009; Grechnev et al. 2015), corotating interaction region shocks (Crooker et al. 1999; Allen et al. 2021), and the heliospheric TS (Decker et al. 2005, 2008). However, these mechanisms are also invoked in astrophysical systems to account for galactic (Aharonian et al. 2004; Helder et al. 2012; Reynolds et al. 2012) and extra galactic (Van Weeren et al. 2016; Kang et al. 2012) cosmic rays. The acceleration process is based on the first order Fermi mechanism, where turbulent magnetic field fluctuations in the upstream and in the downstream regions scatter particles back and forth across the shock front several times, thus increasing the efficiency of the acceleration. Particles can also gain energy by drifting along the shock front, due to the effect of the gradient in the magnetic field magnitude, in the same direction as the motional (\(\mathbf{V}\times \mathbf{B}\)) electric field. In the limit of no pitch-angle scattering, this process is known as shock drift acceleration. Both shock-drift acceleration and first-order Fermi acceleration are subsets of the general process known as diffusive shock acceleration (DSA) (Bell 1978; Lee and Fisk 1982; Jokipii 1987; Liu and Jokipii 2021). It leads to a power law energy distribution \(dN/dE\sim E^{-\gamma}\) where

and \(r=V_{u}/V_{d}\) is the compression ratio of the shock. \(V_{u}\) is the unshocked plasma and \(V_{d}\) is the shocked plasma speed measured in the shock frame. In the strong shock regime \(r=4\), the spectral index \(\gamma =2\). Furthermore, starting from a test particle approach, assuming a planar shock, stationary conditions, and a spatially independent diffusion coefficient, the solution of the transport equation is given by (Drury 1983)

where \(x<0\) indicates the upstream region and \(x\ge 0\) the downstream region. \(g_{u}(p)=f_{d}(p)-f_{u}(p)\) depends on the particle momentum and ensures that the particle distribution is continuous across the shock, \(D_{x}\equiv D_{x}(p)\) is the spatial diffusion coefficient component along the shock normal, and \(f_{d}(p)\) is the finite value reached downstream for the particle distribution function. Thus, the “standard” DSA prediction implies that the energetic particle flux falls off exponentially with increasing distance from the shock front far upstream and remains constant in the downstream side. Another important aspect that needs to be considered is the influence of the turbulent properties of the magnetic field fluctuations on their interaction with energetic particles; for example, it has been shown that at perpendicular shocks a strong correlation exists between the maximum energy reached by accelerated particles and the slope of the power spectral density of the magnetic field, i.e., a steeper power spectrum leads to higher maximum energies (Dosch and Shalchi 2010). It is also well known that magnetic turbulence has a great impact on the properties of the energetic particle transport, since transport regimes can be different in the directions parallel and perpendicular to the average magnetic field (Jokipii 1966; Matthaeus et al. 2003; Webb et al. 2006; Shalchi and Kourakis 2007; Dosch and Shalchi 2010; Shalchi 2011; Zimbardo et al. 2012; Ferrand et al. 2014). At the same time the cosmic ray driven instabilities are able to amplify the background plasma magnetic fluctuations (see e.g. Bell 2004; Bykov et al. 2013a; Schure et al. 2012; Marcowith et al. 2016) to produce strong magnetic turbulence. The effect is widely discussed to operate in DSA of cosmic rays in supernova remnants. Recent direct, multi-spacecraft THEMIS/ARTEMIS observations revealed the effect of a pool of energetic non-thermal ions on the very strong super-adiabatic magnetic field amplification during a passage of a rotational discontinuity through the Earth bow shock in accordance with the results of the hybrid-code simulations (Kropotina et al. 2021).

Numerical simulations of energetic particle transport in synthetic magnetic turbulence have shown that in quasi-slab turbulence the transport along the direction parallel to the mean field can be faster than normal, i.e., superdiffusive, while along the perpendicular direction it can be slower than normal, namely subdiffusive (Zimbardo et al. 2006; Pommois et al. 2007; Shalchi and Kourakis 2007). Further, using a combined slab/2D turbulence model, Hussein and Shalchi (2016) have observed normal diffusion for energetic particles, with mean free paths in excellent agreement with the Palmer consensus (Palmer 1982). The key role of magnetohydrodynamic (MHD) turbulence on the cosmic ray propagation has been investigated in depth by Lazarian and Yan (2014), where perpendicular superdiffusion at scales smaller than the turbulence injection scale has been found. Such a behavior has been ascribed to the turbulent wandering of the magnetic field lines and has effects on the acceleration efficiency at shocks. Thus, the interaction of suprathermal particles with magnetic field turbulence and the related transport properties can affect the shock acceleration process (Perri and Zimbardo 2012a; Giacalone 2013; Zimbardo and Perri 2013; Perri and Zimbardo 2015b). Since the discovery of the DSA mechanism (Axford et al. 1977; Krymskii 1977; Bell 1978; Blandford and Ostriker 1978) many interesting results were obtained in the application of DSA to astrophysical objects of different scale sizes (see for a review e.g. Blandford and Eichler 1987; Berezhko and Krymskiĭ 1988; Jones and Ellison 1991; Malkov and Drury 2001; Ptuskin et al. 2013; Sironi et al. 2015; Marcowith et al. 2016; Blasi 2019). In this review we will primarily focus on the shock acceleration mechanism, giving emphasis to the broad number of observations of shock waves in the interplanetary space, as a place to infer acceleration properties for shocks in the outer heliosphere with potential applications for astrophysics. We review current theoretical interpretations of the observations and highlight the remaining unanswered questions.

The manuscript is organized as follows: in Sect. 2 we review observations of supra-thermal particle fluxes during spacecraft shock crossings and we detail recent interpretations in terms of anomalous transport; we focus on energy spectra, acceleration times, and magnetic field amplification. In Sect. 3 the role of magnetic turbulence in the shock acceleration process is discussed. Section 4 is dedicated to recent observations of flat energy spectra of energetic particles at interplanetary shock waves that challenge standard theory predictions. In Sect. 5 we review a mechanism for electron acceleration at the TS, addressing the role of plasma instabilities that eventually give rise to preferential electron heating. In Sect. 6 we discuss a multiple shock acceleration process as a mechanism connected to the acceleration of cosmic rays at multiple corotating interaction regions in the heliosphere. In Sect. 7 we summarize the content and address unanswered questions.

Although DSA can be considered one of the fundamental processes that can account for cosmic ray acceleration, there are some observational evidences that question its validity. Indeed, Voyager 2 LECP observations of energetic ions accelerated at the solar wind TS indicate that the measured spectral index is consistent with a compression ratio of \(r=3\) according to Eq. (1), while the observed compression ratio is \(r=2\) (Decker et al. 2008). Furthermore, the PAMELA experiment has found evidence of different spectral indices for hydrogen and helium ions (within the 30–1000 GV rigidity range) (Adriani et al. 2011), which is not consistent with the DSA prediction of a spectral index that depends on the compression ratio only. Caprioli et al. (2017), by means of 2D hybrid simulations, have shown that the time-dependent acceleration process of ions at non-relativistic shocks has a dependence on the ratio between the ion atomic mass (\(A\)) and the charge (\(Z\)). At fixed \(Z\) the thermal peak in the energy spectra is shifted towards high energies as \(A\) increases (mass dependence), while the non-thermal tail in the spectra rolls of at the same \(E_{max}/Z\), since DSA is a rigidity dependent process. Also the efficiency of injection into DSA depends on \(A/Z\). At late times all the species tend to show the same universal DSA spectrum. However, cosmic ray spectral differences still remain an open question (Aguilar et al. 2021). Observations of suprathermal particle fluxes from \(\sim 100\) keV to MeV energies accelerated at interplanetary shock waves do not always show an exponential time decrease from upstream of the shock and a constant profile downstream, as predicted by Eq. (2) (Perri and Zimbardo 2008; Sugiyama and Shiota 2011; Giacalone 2012; Perrone et al. 2013; Tessein et al. 2015). In this Section we review recent works that have included anomalous transport of particles in the shock acceleration process to explain some observations in the heliosphere.

2.1 Superdiffusive Transport as Inferred from Energetic Particle Fluxes Far Upstream of Shocks

In a series of papers Perri and Zimbardo (2007, 2008, 2009) have found evidence of energetic ion and electron fluxes, accelerated at interplanetary shocks, decaying as power laws far upstream. In the events analyzed, the magnetic field fluctuations in the upstream region, computed at the time scale in resonance with energetic particles, are independent of the shock distance, so that the observed power law time profile cannot be ascribed to a spatial dependent parallel diffusion coefficient (Perri and Zimbardo 2012b). One example of the event investigated is displayed in Fig. 1, which reports a corotating interaction region (CIR) (indicated in the figure by the vertical dashed line) crossed by the Ulysses spacecraft on 1992 October 11. From top to bottom: the hourly resolution radial solar wind bulk speed and the proton temperature from SWOOPS (Bame et al. 1992), the electron and the proton fluxes in some energy channels (see legend) as measured by LEFS 60/HI-SCALE, both plotted in log-lin axis, and the magnetic field normalized variance computed as \(\sigma ^{2}_{|B|}/|B|^{2}=\sum _{i=1}^{3} \sigma _{i}^{2}/|B|^{2}\) over a time scale \(T=10\) s (where \(\sigma _{i}^{2}=\langle (B_{i}-\langle B_{i} \rangle )^{2} \rangle _{T}\) is calculated over the \(i\)-th component of the magnetic field). Such a time scale corresponds to the Larmor radius \(\rho _{e}\) of electrons of \(\sim 200\) keV propagating in the upstream magnetic field of the event reported in Fig. 1, since \(f=1/T=V_{SW}/(2\pi \rho _{e})\). A closer look at the electron fluxes as a function of the distance from the front in log-log axis is reported in Fig. 2 (Perri et al. 2009), where the far upstream decay can be fitted by a power law \(J(t)\propto t^{-\beta}\). The slopes of the power law found in Perri et al. (2009) for the electron fluxes are \(\beta \in [0.30-0.56]\). Notice that the normalized magnetic field variance in Fig. 1 from about 10 to 100 hr (delimited by the thin red vertical lines) from the shock crossing, i.e., the time range over which the best fits of the electron fluxes have been performed, does not show a significant trend with time in the upstream. Energetic particle power law time profiles have also been detected during the TS crossing by the Voyager 2 spacecraft on August 2007 (Perri and Zimbardo 2009). The spacecraft crossed a supercritical quasi-perpendicular shock (Richardson et al. 2008) and the fluxes of energetic ions from about 540 keV to 3500 keV were observed to decay far upstream of the shock with slope \(\beta \in [0.68-0.71]\), as shown in Fig. 3. All those evidences of power law energetic particle time profiles upstream of shock waves have been ascribed to the possibility of parallel anomalous transport for particles during the shock acceleration process. In particular, a transport faster than normal diffusive, i.e., superdiffusive, can explain those observations. It is worth recalling that superdiffusion can be described by a mean square displacement \(\langle \Delta r^{2} \rangle \) growing with time \(t\) as

where \({\mathcal{D}}_{\alpha}\) is the anomalous diffusion coefficient and \({\alpha}\) is the anomalous diffusion exponent, which is \(>1\) for superdiffusion. The anomalous diffusion coefficient \({\mathcal{D}}_{\alpha}\) has dimensions \([L^{2}/t^{\alpha}]\) different from the standard (i.e., normal) diffusion coefficient. Recently, superdiffusion has been found in several physical systems (Sokolov and Klafter 2005; Zimbardo et al. 2015; Sioulas et al. 2020). For example, numerical simulations have shown that particle transport in presence of magnetic turbulence can be superdiffusive along the mean magnetic field direction (Zimbardo et al. 2006; Pommois et al. 2007; Shalchi and Kourakis 2007). However, when the level of fluctuations is very low, superdiffusion can also characterize the particle transport along the perpendicular direction (Pommois et al. 2007). Further, it has been highlighted how superdiffusive transport has a great influence on shock acceleration (Ragot and Kirk 1997; Perri and Zimbardo 2012a; Zimbardo and Perri 2013; Bykov et al. 2017), leading to a reformulation of the energetic particle spectral index and of the acceleration time (Perri and Zimbardo 2012a; Zimbardo and Perri 2013; Perri and Zimbardo 2015a,b), as it will be discussed in Sect. 2.4.

Overview of a CIR event detected by the Ulysses spacecraft on 1992 October 11. From top to bottom, the panels show the hourly resolution radial solar wind bulk speed and the proton temperature from the SWOOPS experiment (Bame et al. 1992), the electron fluxes in the energy channels 42–65 keV, 65–112 keV, 112–178 keV, and 178–290 keV, the proton fluxes in the energy channels 546–761 keV, 761–1223 keV, 1223–4974 keV, and 178–290 keV as measured by LEFS 60/HI-SCALE, the 10 sec resolution total magnetic field normalized magnetic field variance. The shock jump is indicated by the vertical dashed line. The thin vertical red lines delimit the time range where the electron fluxes decay as a power-law in time (see Fig. 2)

Full size image

Electron fluxes upstream of the reverse shock of 1992 October 11 in log–log scale. Dashed lines represent the best power-law fits. Energies are indicated in the figure legend (Adapted from Perri et al. (2009))

Full size image

Ion fluxes upstream of the termination shock of the solar wind, as measured by the Voyager 2 spacecraft, in log–log scale. Dashed lines represent the best power-law fits. Energies are indicated in the figure legend (Adapted from Zimbardo and Perri (2010))

Full size image

There are several observations that indicate anomalous transport in astrophysical systems: the transport of non-relativistic electrons accelerated in solar events (Lin 1974), the analysis of transport properties from X-ray synchrotron emission from relativistic electrons accelerated at supernova remnants (Perri et al. 2016; Perri 2018), as well as the analysis of diffuse emission from electrons in the Coma cluster, which claims for nearly ballistic motion of energetic electrons (Ragot and Kirk 1997).

Therefore, both observations and numerical studies have pointed out how particle transport largely influences the efficiency of shock acceleration. Indeed, energetic particles can be scattered by the magnetic field fluctuations in resonance. The efficiency of such a scattering strongly depends on the fluctuation amplitude at the resonant scales. In case of very low \((\delta B/B)^{2}\) at the resonant scales, the particle scattering can be reduced, leading to an almost scatter-free propagation, i.e., non diffusive transport (Bieber et al. 1994; Zimbardo et al. 2012).

The analysis method adopted by Perri and Zimbardo (2007, 2008) and Zimbardo and Perri (2013) for deriving the energetic particle fluxes under the assumption of superdiffusive transport is based on the propagator formalism (Zumofen and Klafter 1993). A propagator \(P(x-x', t-t')\) gives the probability of observing a particle at \((x, t)\) if injected at \((x', t')\). Its mathematical form does depend on the particle transport regime: normal diffusion is properly described by a Gaussian propagator with a width proportional to the diffusion coefficient and to \(t-t'\) (Metzler and Klafter 2000); superdiffusion is modeled by a power-law propagator (Zumofen and Klafter 1993):

with \(2 <\mu < 3\). A “macroscopic” description of the anomalous transport by means of fractional Fokker-Planck equations in space and time is described in Engelbrecht et al. (2022) in this volume.

Under the hypothesis of a planar 1D shock moving at constant speed \(V_{\mathrm{sh}}\) with respect to the observer into a statistically homogeneous medium characterized by an almost constant level of magnetic field turbulence (Lee 2005), and injecting the energetic particles at the shock, the particle density with energy \(E\) can be calculated via

where \(S(x',E,t')=\varPhi _{0}(E) \delta (x'-V_{\mathrm{sh}}t')\) is the source term that depends on the flux of particles of energy \(E\) injected at the shock front. Further, assuming a shock that comes from \(x<0\) for \(t< 0\), the profile in equation (5) can be easily calculated by inserting the form of the propagator that describes the process. In case of normal diffusion, i.e., for a Gaussian propagator, we recover equation (2). So that for an observer at \(x=0\) and a shock coming from far upstream (i.e., \(t<0\)) we end up with \(n(0,E,t)\sim \exp{(V_{\mathrm{sh}}^{2}t/D})\). In case of superdiffusion, i.e., for the power-law propagator, a power-law decay is derived as (Perri and Zimbardo 2008)

where the exponent \(\beta = \mu -2\) ranges between \(0 <\beta < 1\), being \(2 <\mu < 3\). Notice that the exponent \(\beta \) is directly related to the exponent of the mean square displacements in superdiffusive transport (see equation (3)), namely \(\alpha = 2-\beta \) (Perri and Zimbardo 2007, 2008).

2.2 Kappa Distributions and Non-thermal Particle Populations

In the context of the acceleration of suprathermal particle populations an analytical description of their velocity or momentum or energy spectra is of interest. Particular useful tools are so-called kappa distributions, which have been employed for solar energetic particles (Oka et al. 2013; Effenberger et al. 2017), for the suprathermal solar wind components (Lazar 2017; Yoon et al. 2018), for pick-up ions (Fahr et al. 2016; Heerikhuisen et al. 2019), for cosmic rays (Treumann and Baumjohann 2018) as well as for studies of various wave modes in the corresponding plasma milieus (e.g., Husidic et al. 2020, and references therein).

In all of these studies kappa distributions of the form

or variants thereof have been employed. Therein \(v\) is the particle speed, \(\theta _{M}\) is a reference speed corresponding to the thermal speed of the associated Maxwellian, and \(A_{\kappa}\) is a normalisation constant.

The above representation can be referred to as an (isotropic, non-relativistic) Olbertian (see the left panel of Fig. 4) since it dates back to Olbert (1968, see also Vasyliunas (1968)) for a detailed overview see Lazar and Fichtner (2021) and for relativistic kappa distributions see Han-Thanh et al. (2022). In the following, we discuss a few topics that may, so far, not be as well known as they deserve to be. This will comprise a brief review of different versions of kappa distribution, their application-related limitations, and their connection to nonlinear and anomalous diffusion.

Left panel: The (isotropic) Olbertian Kappa (red solid line). Right panel: The modified kappa distribution (blue solid line). Both standard kappa distributions (SKDs) are shown for the same kappa value along with their associated Maxwellians (dotted and dash-dotted, respectively). Note that for the Olbertian the suprathermal tails form on expense of its core population resulting in a higher temperature than that of the associated Maxwellian, while in case of the modified kappa distribution both have the same temperature. Adapted from Lazar and Fichtner (2021)

Full size image

2.2.1 Different Versions of Kappa Distributions

While it is straightforward to generalize any given isotropic kappa distribution to its anisotropic counterpart (Scherer et al. 2019; Lazar and Fichtner 2021), for simplicity we consider only the former here. Closer inspection of the above Olbertian reveals that it describes suprathermal distributions characterized (i) by the same constant reference speed \(\theta = \theta _{M}\) and (ii) by temperatures that are increasing with decreasing kappa parameter according to

with \(m\) being the particle mass, \(k_{B}\) the Boltzmann constant, and \(T_{M}\) the temperature of the associated Maxwellian. This kappa dependence of the temperature is a feature that was desired by Olbert (1968) to quantify the deviations of observed energetic (in his case magnetospheric) electrons from thermal equilibrium and has, more recently, been shown to hold in the solar wind as well (Lazar et al. 2017). Note that, of course, the velocity dependence can straightforwardly translated in to one on momentum or on kinetic energy, so that in the following we continue to use the velocity formulation.

The first variant of the Olbertian was originally introduced by Matsumoto (1972) and consists in an interchange of the kappa-dependence of the reference speed and the temperature, so that the latter becomes independent of kappa:

This modified kappa distribution (see the right panel in Fig. 4), which has been favoured in recent years by Livadiotis (2015), is, however, only applicable to specific physical systems as discussed in Lazar et al. (2016). One of the features of this particular distribution to be explained is the enhancement of particles in the core region compared to the associated Maxwellian, which is, thus, not approximating the core, as in the case of the Olbertian (see Fig. 4).

Beyond these principal differences between these two so-called standard kappa distributions (SKDs), which have significant implications, e.g. for the properties of waves maintained in a suprathermal plasma (Lazar et al. 2015), both alternatives face common problems. First, the velocity moments, which are necessary for the derivation of a macroscopic fluid description from kinetic theory, are diverging for \(l \geq 2\kappa -1\) with \(l\) denoting the order of the moment. Second, all moments are diverging for \(\kappa \leq 3/2\). On the one hand, this prevents the formulation of a consistent hydrodynamical treatment of a given physical system and, on the other hand, power laws with \(\kappa \lesssim 2\) cannot justifiably represented with a SKD. Third, it has been demonstrated that a SKD implies a non-extensive entropy (Tsallis 1988), which is a questionable property (Nauenberg 2003).

In order to remove these conceptual difficulties, Scherer et al. (2017) introduced the regularized kappa distribution (RKD, see Fig. 5), which equips the SKDs with an exponential, Maxwellian-like cut-off:

Again, \(A_{RKD}\) is a normalisation factor. The cut-off is controlled with the parameter \(\alpha \). Such a cut-off can be motivated theoretically by the finiteness of scales: in a given physical system an acceleration process takes place within a maximum spatial scale and within a finite period of time. These constraints limit any acceleration process correspondingly, so that a resulting spectrum of accelerated particles cannot extend to arbitrarily high velocity (momentum, energy) but rather must be cut off. Given that a consistency of kinetic and fluid theory requires the existence of all velocity moments, such cut-off must be exponential.

The (isotropic) Olbertian (red solid line, for \(\kappa = 2\)) compared to the corresponding RKD (blue dashed line) for \(\alpha = 0.2\), and to the Mxwellian limit. Taken from Lazar and Fichtner (2021)

Full size image

The RKD has not only the advantage of having finite velocity moments at arbitrary order, but it also avoids all divergences for \(\kappa \leq 3/2\) so that it is defined for all positive kappa values. The latter property makes it possible to treat significantly harder power laws than a SKD would allow. Furthermore, almost surprisingly but very satisfactorily, it is consistent with an extensive entropy (Fichtner et al. 2018) and, thus, with the fundamental laws of thermodynamics. Consequently, the entropy of a RKD can then be compared to that of a Maxwellian and, as to be expected, it is always lower for a given system.

While the RKD can, as outlined above, be constructed on purely theoretical grounds, it can also be motivated by observations. Scherer et al. (2017) pointed out that both pick-up ion and anomalous cosmic ray spectra have previously been described by RKD-like distributions, see, e.g., Fisk and Gloeckler (2012) and Steenberg and Moraal (1999), respectively. Further information and various recent applications of the concept of kappa distributions to suprathermal particle populations can be found in Lazar and Fichtner (2021). We note that the acceleration mechanism of anomalous cosmic rays was the subject of the review by Giacalone et al. (2012). It is briefly been summarized again in the review of Giacalone et al. (2022) (this volume), which includes some recent insights.

2.2.2 Kappa Distributions and Generalized Diffusive Transport

It is also interesting to note that kappa distributions appear to be related to both nonlinear and anomalous diffusion. Tsallis and Bukman (1996) have considered a nonlinear Fokker-Planck (i.e. diffusion) equation and, by using the SKD-related concept of non-extensive entropy, have demonstrated that this way anomalous transport, i.e. sub- and superdiffusion, can be described. They refer to this type of anomalous diffusion, that was also found by Litvinenko et al. (2017), as ‘driven’ or ‘correlated’ and distinguish it from ‘Levy-like’ anomalous diffusion described by linear but fractional Fokker-Planck/diffusion equations, see, e.g., Stern et al. (2014), Perri et al. (2015), Zimbardo et al. (2015, 2017). The significance of RKDs in this context still needs to be explored.

Combining the observations of power-law energetic particle time profiles described above with the findings about nonlinear Fokker-Planck equations, one can indeed, in a first approach, attempt to fit the time profiles on the basis of a simple one-dimensional nonlinear diffusion equation of the form

Doing so for shock-accelerated electrons observed by the Ulysses spacecraft gives the results shown in Fig. 6. While such studies are still in their infancy, it is, nonetheless, obvious that solutions of a nonlinear diffusion equation allow to fit the data reasonably well. Thus, the exploration of such equation in application to energetic particle populations looks promising and is presently ongoing (e.g., Litvinenko et al. 2019; Walter et al. 2020).

Differential intensity of energetic electrons accelerated at corotating interaction regions as observed by the Ulysses spacecraft in the upstream region. For the three events the mean electron energy is given along with the parameter \(\mu \) obtained from fits (blue line) of the solution of Eq. (11) to the data (red symbols). For the latter see, e.g., Perri et al. (2015) or Perri and Zimbardo (2012b)

Full size image

2.3 The Possibility of a Non-markovian Pitch-Angle Scattering

Now, it is important to address at least one of the physical arguments that are at the basis of the anomalous diffusion. A relation between pitch-angle scattering and parallel spatial diffusion in non-diffusive regimes has been discussed in Shalchi et al. (2012). Here we start from some evidence that emerged in both numerical studies and natural environments. Recently, simulating the transport of energetic particles in static 3D isotropic turbulence, Pucci et al. (2016) have highlighted how particles that diffuse normally in the direction parallel to the mean magnetic field, exhibit a power-law distribution of the pitch-angle scattering times with no characteristic values close to the quasi-linear prediction (Kennel and Petschek 1966; Amato 2014). Such a distribution tend to broaden over a larger range of values of the scattering times for high level of magnetic turbulence intermittency (i.e., when the degree of spatial non-homogeneity of the turbulence increases) and for \(\delta B/B_{0}\ll 1\) (Pucci et al. 2016; Perri et al. 2019). A power-law distribution of the particle scattering times have also been found for energetic particles upstream of interplanetary shock waves (Perri and Zimbardo 2012b; Perri et al. 2021), where the scattering times have been evaluated via the power of the magnetic field fluctuations at the particle resonant scale (Crooker et al. 1999). A theoretical explanation has been formally given in Zimbardo and Perri (2020), where the power-law distribution of the particle scattering times has been shown to be the basis of anomalous pitch-angle scattering. This eventually gives rise to spatial superdiffusion in the direction parallel to the mean magnetic field.

In another numerical study that simulates the collisionless plasma close to a shock region, Riordan and Pe’er (2019) found a distribution of the particle scattering time that extends over several decades, depending on the magnetic turbulence level \(\delta B/B_{0}\) (spanning from weak to strong turbulence).

Thus, pitch-angle scattering of energetic particles in presence of magnetic turbulence can be modeled as a non-Markovian process with a power-law distribution of the particle scattering times. Indeed, the distribution of the particle free path lengths can be described in the framework of the Lévy random walks, i.e.,

with \(C\) being a normalization constant, \(\ell \) is the particle jump length, \(v_{||}\) is the particle’s parallel speed, and \(\tau \) represents the scattering time. Equation (12) is valid for large \(\ell \) (Metzler and Klafter 2004). Integrating over the path lengths, it is possible to get \(\psi (\tau )= C |v_{||}|^{-\mu} \tau ^{-\mu}\). For \(2<\mu <3\) superdiffusive transport is obtained. Thus, Zimbardo and Perri (2020) has shown that assuming a distribution of the particle scattering times \(\psi (\tau ) \propto \tau ^{-(1+\beta )}\), with \(0<\beta <1\), it is possible to formally derive a non-Markovian Fokker-Planck equation, where the time derivative is replaced by a fractional time derivative (Caputo 1967). By solving the fractional Fokker-Planck equation for two different forms of the pitch-angle scattering coefficient, the parallel velocity autocorrelation function has been found to decay as a power-law in time (long range correlations) and the spatial transport along the mean field direction to be superdiffusive.

2.4 The Influence of Superdiffusive Transport on the Shock Acceleration Process

Starting from a first order Fermi process, where test particles of velocity \(v\) can cross the shock front several times, gaining a relative energy at each crossing \(\Delta E/E=(4/3) (V_{u}-V_{d})/v\), where \(V_{u}\) and \(V_{d}\) are the fluid speeds in the upstream and downstream regions, respectively, it is possible to include the effect of superdiffusive motion in the shock acceleration process. Perri and Zimbardo (2012a), Zimbardo and Perri (2013) have generalized the DSA mechanism including anomalous, faster, transport. Assuming a shock coming from \((x,t)=(-\infty ,-\infty )\) that is moving to the observer according to \(x'=V_{\mathrm{sh}}t'\), where the shock speed \(V_{\mathrm{sh}}\) is kept fixed to a constant value in the downstream frame, Perri and Zimbardo (2012a) have derived the probability of escape for particles from the acceleration region. Such a probability is directly related to the spectral index and is given by the ratio of the fluxes upstream and downstream, namely \(P_{\mathrm{esc}}=\varPhi _{d}/\varPhi _{u}=n_{d}V_{d}/(n_{0}v/4)\). Notice that the flux upstream is derived under the usual isotropic conditions (Drury 1983; Gaisser 1990), and \(n_{0}\) is the density of particles coming from upstream with speed \(v\) at the shock front. Further, following Kirk et al. (1996) and assuming stationary conditions, far downstream the particle density is proportional to the flux of particles at injection, i.e., \(n_{d}=\varPhi _{0}(E)/V_{\mathrm{sh}}\). The crucial point is to derive the particle density at the shock \(n_{0}\) assuming a superdiffusive motion. This has been computed by using equation (5) with the correct form of the propagator (i.e., see equation (4)). Setting \(\int _{-\infty}^{\infty}P(x,t)dx=1\), Perri and Zimbardo (2012a), Zimbardo and Perri (2013) end up with

where the relation that links \(\mu \) to the exponent of superdiffusion \(\alpha \) in equation (3), namely \(\mu =4-\alpha \), has been used. Equation (13) basically relates the density at the shock to the density far downstream. It is worth addressing that Equation (13) tends to its diffusive limit \(n_{d}/n_{0}=1\) for \(\alpha \rightarrow 1\). This directly leads to an escape probability \(P_{\mathrm{esc}}=(n_{d} V_{d})/(n_{0} v/4)=8 (\mu -2)/(\mu -1)V_{d}/v=8 (2- \alpha )/(3-\alpha )V_{d}/v\).

2.4.1 Particle Energy Spectra

From equation (13), it is straightforward to derive the exponent of the integral energy spectrum (i.e., \(N(\ge E)\propto (E/E_{0})^{-\tilde{\gamma}}\)) of relativistic particles,

Again, for \(\alpha \rightarrow 1\) (diffusive transport) equation (14) tends towards \(\tilde{\gamma}=3/(r-1)\) (Fisk and Lee 1980), and then the exponent of the differential energy spectrum towards \(\gamma =\tilde{\gamma}+1=(2+r)/(r-1)\). Perri and Zimbardo (2012a), Zimbardo and Perri (2013) also derived the spectral indices for non relativistic particles. Considering a relative gain in momentum during the acceleration process \(\Delta p/p=(4/3) (V_{u}-V_{d})/v\), they get

The exponent of the differential energy spectrum is found to be \(\gamma =3/(r-1)(2-\alpha )/(3-\alpha )+1\) (Zimbardo and Perri 2013). Basically, the superdiffusive correction in equations (14) and (15) is due to the ratio between the density far downstream and the density at the shock as in equation (13). Examples of the spectral indices as a function of the compression ratio of the shock are reported in Fig. 7.

Energy spectral index as a function of the compression ratio of the shock. The spectral index predicted by assuming a superdiffusive transport for energetic particles are displayed for four different values of \(\alpha \). The DSA prediction is also shown with a black solid line both for relativistic (left panel) and for nonrelativistic (right panel) particles. (Adapted from Zimbardo and Perri (2013))

Full size image

Superdiffusion leads to smaller spectral indices and thus harder energy spectra than DSA. This has been interpreted as due to the fact that downstream of the shock two physical effects are at play, namely advection (\(\Delta x\sim V_{d} t\)) and random superdiffusive motion \(\langle \Delta x^{2}\rangle \sim t^{\alpha}\), with \(\alpha >1\). In the limit of long times even in the presence of advective motion, energetic particles can propagate upwind more easily than in the DSA framework, so that they have a higher probability to return to the shock. This will also lead to shorter acceleration times, as it will be discussed in Sect. 2.4.2. Notice that we have not considered here the effect of adiabatic cooling, which can have influence on the acceleration process when the particle mean free-path is large upstream (Ruffolo 1995).

An example of particle differential energy spectrum calculated averaging the ion fluxes within a time window of about 80 min downstream is displayed in Fig. 8 for an interplanetary shock crossing of the Advanced Composition Explorer (ACE) spacecraft at 1 AU that occurred on 2010 August 3. Both the DSA prediction (black line) and the spectra derived within the superdiffusive framework are shown (red lines). In these events upstream fluxes have been found to follow a power-law similar to the ones displayed in Fig. 2. Notice that the two red lines are calculated by considering the lowest and the highest \(\alpha \) values obtained by fitting the upstream energetic particle fluxes in different energy channels (see Sect. 2.1). A better agreement is found with the energy spectrum predicted assuming superdiffusive motion of energetic particles.

Differential energy spectrum of ions in five different energy channels for a shock event detected by ACE on 2010 August 3rd. The solid line is the spectrum predicted by DSA, while the red lines are the ones predicted assuming superdiffusive motion for energetic particles; they have been computed using the highest and the lowest \(\alpha \) values obtained from the best power-law fits of the ion fluxes in the upstream region. Error bars represents the size of each energy bin of the detector

Full size image

2.4.2 Acceleration Times

The occurrence of superdiffusive transport for energetic particles during the shock acceleration process affects also the acceleration time. In the framework of DSA the acceleration time is given by

where \(D_{ux}\) (\(D_{dx}\)) is the upstream (downstream) diffusion coefficient perpendicular to the shock front (e.g., Drury 1983). They are related to both the parallel and the perpendicular diffusion coefficient, although in case of a dominant parallel transport it is possible to get \(D_{x}\sim D_{||}\cos ^{2}{\theta _{Bn}}\), where \(\theta _{Bn}\) is the angle between the normal to the shock front and the upstream magnetic field. Perri and Zimbardo (2012a) have derived a similar relation for the acceleration time in presence of superdiffusion, by equating the macroscopic plasma advective motion, \(\Delta x=V \Delta t\), with the microscopic particle superdiffusive motion \(\langle \Delta x^{2}\rangle \sim {\mathcal{D}}_{\alpha x}\Delta t^{\alpha}\) along the direction of propagation of the shock. The residence time within the acceleration region can be considered equal to the ratio between the number of particle present in a region of length \(\Delta x\) and the incident energetic particle flux coming from upstream (e.g., Drury 1983), i.e., \(t_{res}=n\Delta x/(nv/4)\). Substituting \(\Delta x\) obtained from the balance of advective and superdiffusive motions and resolving for the acceleration time \(t_{\mathrm{acc}}^{S}=t_{res}/(\Delta p/p)\) it is easy to obtain (Perri and Zimbardo 2012a)

The key point in such an expression is the derivation of the anomalous diffusion coefficient \(D_{\alpha}\). Indeed, as shown in Zimbardo and Perri (2013), it is related to the exponent of superdiffusion \(\alpha \), to the particle speed, and also to a length \(\ell _{0}\) that comes from the probability \(\psi (\ell , t)\) of making for a particle a jump of length \(\ell \) in a time \(t\) (Klafter et al. 1987; Zumofen and Klafter 1993) shown in equation (12), namely

which is the parallel anomalous diffusion coefficient (\({\mathcal{D}}_{ \alpha}={\mathcal{D}}_{\alpha x}/\cos ^{2}{\theta _{Bn}}\)). It is worth noting that for \(\ell \rightarrow 0\) \(\psi (\ell ,t)\) is bell-shaped, so that we have no divergence in equation (12). Thus, \(\ell _{0}\) represents the minimum length at which the power-law probability for the jump lengths holds in the particle random walk. As shown in Perri and Zimbardo (2015b), it is possible to derive \(\ell _{0}\) directly from spacecraft in-situ observations. Moreover, coming back to the propagator formalism, it has been shown that in case of superdiffusion the propagator in equation (4) holds far from the source of energetic particles, while close to it the form of the propagator changes towards a modified Gaussian distribution (Zumofen and Klafter 1993). Zumofen and Klafter (1993), Zimbardo and Perri (2013), Perri and Zimbardo (2015b) have defined a scaling variable \(\xi =(x/\ell _{0})/(t/t_{0})^{1/\mu}\) (where \(t_{0}=\ell _{0}/v\)) such that for \(\xi \ll 1\) particles are close to the shock and for \(\xi \gg 1\) particles are far upstream from the source. Perri and Zimbardo (2015b) have suggested that the change in the functional form of the propagator corresponds in terms of spacecraft observations to a change in the energetic particle fluxes in the entire upstream region (from close upstream towards far upstream). This happens when the spacecraft detects a power-law decay in the particle flux at a certain distance/time from the shock, i.e., \(t_{\mathrm{break}}^{\mathrm{s/c}}=x_{\mathrm{break}}/V_{\mathrm{sh}}^{\mathrm{s/c}}\), where \(V_{\mathrm{sh}}^{\mathrm{s/c}}\) is the relative speed between the shock front and the satellite.

In Fig. 9 it is reported an example of identification of \(t_{\mathrm{break}}\) directly from observations. Figure 9 is relevant to a CIR crossing by Ulysses on 1993 May 10. Thus, if at \(\xi =1\) the two forms of the propagator match, Perri and Zimbardo (2015b) have derived

Using equation (19) an estimation of equation (18) from satellite data can be obtained. By applying such an analysis to a series of interplanetary shock waves and to the TS of the solar wind, Perri and Zimbardo (2015b) have found that the acceleration times for energetic protons in the presence of superdiffusive transport are largely reduced with respect to the ones derived from equation (16). They are also shorter than the typical shock lifetimes. This has been ascribed to the scale-free nature of the particle jump lengths (see equation (12)), since the probability \(\psi \) of very long displacements is non-negligible with respect to a Gaussian distribution. Further, also very short jumps are highly probable increasing both the probability to return to the shock for particles travelling far upstream, and the probability for particles close to the front to cross the shock several times. For interplanetary shocks at 1 AU, superdiffusion reduces the acceleration time to a fraction of day, while for corotating interaction regions in the outer heliosphere and for the TS the acceleration time has been estimated to be ∼ few days.

Electron flux at about 100 keV upstream of a CIR detected by the Ulysses spacecraft on 1993 May 10 as a function of the time distance from the front. The flux decay as a power-law from 3 to 100 hrs upstream from the shock. Closer to the front the electron flux tends to become flatter. This change in shape is marked as \(t_{\mathrm{break}}^{\mathrm{s/c}}\) by a red dashed line. (Adapted from Perri et al. (2015))

Full size image

2.4.3 Magnetic Field Amplification in Superdiffusive Shock Acceleration in Supernova Remnants

Microscopic plasma modeling of cosmic ray (CR) acceleration by collisionless shocks (both PiC and hybrid) demonstrated an efficient formation of non-thermal particle population simultaneously with the accompanied spectra of magnetic fluctuations (for a review see Marcowith et al. 2016). To model cosmic ray acceleration in supernova shocks computer simulations with a wide dynamical range are needed and instead of computationally expensive PiC simulations one can use simplified kinetic equations or Monte Carlo models. The models use some more or less justified recipes to connect the CR transport to magnetic fluctuation spectra. A nonlinear Monte Carlo model of collisionless shock flow structure with super-diffusive propagation of high-energy particles coupled to particle acceleration and magnetic field amplification was presented by Bykov et al. (2017). The cosmic ray anisotropy provided strong non-adiabatic amplification of magnetic fluctuations due to the resonant and non-resonant instabilities in the model that allows the Lévy-walk type propagation of CRs in the shock flow. The model uses a semi-phenomenological equation to calculate the local spectral energy densities of magnetic fluctuations assuming a balance between the magnetic field growth and the spectral cascading. Contrary to the widely used advection-diffusion equations, the Monte Carlo technique includes the full angular anisotropy of the particle distribution at all spatial positions. A distinctive feature of the non-linear Monte Carlo model of super-diffusive shock acceleration is the relatively strong quadruple anisotropy in the precursor CR distribution, which produces a nonresonant mirror-type instability that amplifies compressible wave modes with wavelengths longer than the gyroradii of the high-energy protons produced by the shock. In Fig. 10 taken from Bykov et al. (2017) we illustrated the magnetic field amplification by CR-driven instabilities in shocks with super-diffusive transport. The strong super-adiabatic amplification of the magnetic field is apparent both in the case of the turbulence cascading to the short scales and with suppression of the cascade along the mean large scale magnetic field.

Spectral energy densities of the CR-driven magnetic fluctuations derived in the non-linear model of super-diffusive shock acceleration by Bykov et al. (2017). The spectra are measured at three positions at a quasi-parallel shock: (a) the free escape boundary (FEB), (b) \(z = -0.01z\) FEB, and (c) in the downstream region. The upper panel shows the spectra simulated in the model with quenched turbulent cascade as it is the case in the models of anisotropic turbulence, while the bottom panel illustrates the magnetic field amplification in the models with the Kolmogorov type cascade

Full size image

A pivotal role in the evolution of the shock waves is played by their interaction with magnetic turbulence of the ambient medium. Indeed, it is well known that the magnetic field turbulent cascade spontaneously generate coherent structures and wave modes over a broad number of spatial scales (Bale et al. 2005; Greco et al. 2010; Gary et al. 2012; Karimabadi et al. 2013). Shock waves that propagate through such a complex medium will be highly distorted with amplification of the upstream turbulent fluctuations at the shock (Zank et al. 2002), a spatial (and time) change in the compression ratio of the shock (Giacalone and Neugebauer 2008), and in the geometry along the shock surface (Giacalone 2017; Kajdič et al. 2019). Such effects are worth being investigated since they can have a large affect on the process of particle acceleration at the shocks. Using multi-spacecraft data techniques, distorted fronts were observed in the near-Earth environment as non-planar rippled surface with a radius of curvature consistent with the interplanetary magnetic field correlation length and characterized by different energetic particle fluxes at different spacecraft locations (Neugebauer and Giacalone 2005).

Recently, Kajdič et al. (2019), taking advantage of the multispacecraft Cluster mission, reported four observations of irregular shock surfaces in the near-Earth environment. During the four encounters, the inter-spacecraft distance were varying, ranging from few ion inertial lengths to about \(100 d_{i}\) (\(d_{i}=c/\omega _{pi}\), where \(c\) is the speed of light and \(\omega _{pi}\) is the ion plasma frequency). In terms of distances they correspond to tens km up to \(10^{4}\) km. This was a unique opportunity to explore the spatial variations of the shocks at different spatial scales. Further, the shocks analyzed exhibit a geometry that varies from quasi-parallel to quasi-perpendicular. Kajdič et al. (2019) computed the normal to the shock front using a minimum variance analysis technique (Sonnerup and Cahill 1968) and consequently estimated the \(\theta _{bn}\) angle, with an uncertainty of about \(5^{\circ}\). In addition, thanks to the multispacecraft measurements, they were able to compute the angle between two shock normals for pairs of satellite: such an angle was found to be small for small interspacecraft separation and up to \(30^{\circ}\) for largely separated satellites. Thus, the angle \(\theta _{bn}\) tends to vary from probe to probe, especially for large spacecraft separations. These evidences imply that the shock geometry changes locally at the different satellite positions, probably caused by an irregular shock surface. Data were also supported by a 2D hybrid simulation of a shock with an Alfvénic Mach number \(M_{A}\sim [4,6.6]\). Pre-existing compressive structures in the upstream, once convected toward the shock front, tend to rise in amplitude and merge with the front causing a local shock reformation. Such a process makes the shock surface highly corrugated and irregular, with a local change in the shock geometry. Thus, the presence of upstream structures can cause strong deformation of the shock front when interacting with them. Giacalone and Jokipii (2007) investigated a similar effect induced by turbulent density fluctuations upstream of a hydrodynamic shock, the interaction between the shock wave and these turbulent fluctuations also caused the formation of a turbulent postshock region and magnetic field amplification.

A more realistic situation has been explored by Trotta et al. (2021) using a 2D hybrid code, where a turbulent jet generated by a 2.5D compressible MHD simulation can interact with an oblique supercritical shock wave. This aims at reproducing the shock dynamics when expanding in a highly turbulent medium. The level of turbulence is controlled by the values of \(\delta B/B_{0}\) and the upstream condition is passed to the hybrid code when the MHD turbulence is already fully developed and coherent structures are formed. For typical heliospheric conditions, i.e., \(\delta B/B_{0}\sim 0.8\), the interaction between upstream turbulent structures/waves and the shock causes a high distortion of the propagating front. In such a scenario, energy spectra of particles are characterized by a high energy field aligned beam that moves towards higher energy and spread as the \(\delta B/B_{0}\) increases and by a formation of a very low energy population of particles that is probably emerging due to trapping mechanisms between particles and coherent turbulent structures. The effect of turbulence on the particle spatial and velocity diffusion has deeply been explored by using a coarse-graining of the Vlasov equation, in order to get an equation for the reduced particle density. The terms in the equation corresponding to spatial and velocity diffusion were computed within the simulation domain finding a large enhancement in both space transport and velocity transport upstream as \(\delta B/B_{0}\) increases. This implies that the acceleration process is not only limited to the shock front but can involve the turbulent fluctuations in the box. The velocity diffusion was also found to be correlated with the electric field parallel to the mean field direction, which is a source of particle energization upstream. Interestingly, it was discussed that the fast particle reflection at the shock can lead to superdiffusion upstream (as already pointed out in a 2D hybrid simulation by Trotta et al. (2020), where superdiffusive transport was found for energetic particles in the upstream region), while the emergence of coherent structures can act as particle traps and accelerators.

Exploring years of ACE measurements, Tessein et al. (2013, 2015) found a good correlation between the presence of energetic particle enhancements and magnetic field discontinuities, as computed from the Partial Variance of Increments (PVI) technique (Greco et al. 2008, 2018). PVI basically measures the local magnetic energy over a given time scale with respect to the average magnetic energy of the data-set. Beside particle flux enhancements close to shock crossings, discontinuities characterized by high values of PVI (strong variations in the magnetic field) have associated high energetic particle fluxes. Such coherent structures are also detected in the vicinity of shocks, thus probably contributing to the particle acceleration process. Tessein et al. (2015) also pointed out how energetic particle flux decays can be observed near magnetic discontinuities that act as a magnetic barrier. In either case it is an important open question how structures emerging from a turbulent background can influence the processes of particle acceleration and energization.

These studies point out the fundamental role of magnetic turbulence in the shock acceleration process from the thermal pool towards supra-thermal energies (Giacalone 2005a,b).

Energetic particle observations that question standard theory results were investigated by Lario et al. (2018). They found several large CME-shock related SEP events which revealed a common, but previously unreported characteristic that the SEP energy spectra upstream of the shock were independent of energy over a fairly broad range of energy. An example is shown in Fig. 11. The SEP event in this figure was included in the analysis of Lario et al. (2018), but we present a re-formatted version for this paper in order to highlight a few key points. The figure shows data from the MAG and EPAM instruments onboard ACE. Panel (a) shows the fluxes of energetic ions with the energies shown in the legend. This data comes from the LEMS120 detector on EPAM. Panels (b) and (c) show the magnetic field magnitude and components, respectively. The vertical dashed lines and numbers are used to indicate some key time periods. The shock arrives at the end of period 2, and period 3 is associated with the CME sheath region, between the shock and CME flux rope, which crosses the spacecraft at the dashed line between periods 3 and 4. Period 4 is the CME flux rope. Consider time period 1 in Fig. 11. This is a ∼ 17-hr period, over which the flux of 68-keV - 1.9-MeV ions have very nearly the same flux at all energies. Thus, in terms of flux, this is a flat energy spectrum. During period 2, the fluxes increase rapidly to where they peak at the shock, and do so in a manner which is dependent on energy since just downstream of the shock, the fluxes clearly depend on energy, despite that they are nearly constant in time.

ACE MAG and EPAM observations of a fast interplanetary shock driven by a coronal mass ejection. (a) Differential flux of ions from the LEMS120 detector as a function of time, with the energies indicated in the legend; (b) magnitude of the magnetic field; and (c) the three components of the magnetic field in the usual RTN coordinate system. The vertical dashed lines indicate regions of interest. The CME-driven shock crossed the spacecraft at the end of period 2. Period 3 is region between the CME flux rope and shock, while period 4 is the CME-flux rope itself. This figure was generated for this paper, but shows data that was also presented in Lario et al. (2018)

Full size image

Since the differential flux, \(J\) (flux per energy) is related to the distribution function, \(f\) according to \(f=p^{2}J\), where \(p\) is the momentum, we note that an energy-independent differential flux spectrum corresponds to \(f\propto p^{-2}\). Note also that the fluxes all rise at the same rate from the background level near the start of the time period, to the end of the time period, increasing by a factor of about ten. Thus, it seems that not only is the energy spectrum independent of energy, the diffusion coefficient associated with the transport of these particles must also be independent of energy. This is an important constraint on particle-transport theory. For example, the prediction from quasi-linear theory (e.g. Jokipii 1966; Hasselmann and Wibberenz 1970; Giacalone and Jokipii 1999) for the parallel diffusion coefficient of charged particles moving in a turbulent magnetic field having a Kolmogorov-like power spectrum gives \(\kappa _{\parallel}\propto p^{4/3}\) (for non-relativistic energies). This is clearly inconsistent with the observations. It is also important to note that the start of this period is so far in advance of the shock that it is reasonable to expect that the magnetic turbulence during this period is likely quite similar to the ambient turbulence in the solar wind. Thus, the observations of flat-energy spectra upstream of shocks provide important clues into the nature of energetic particle transport in the ambient solar wind. It is noteworthy that there is a significant discrepancy between the estimates of the diffusion coefficient (or mean-free path of scattering) from quasi-linear theory and observations (e.g Palmer 1982; Bieber et al. 1994).

Previous numerical modeling of SEP events reveal some results which show a flat energy spectrum upstream of the shock (c.f. Fig. 18 of Prinsloo et al. 2019). These authors suggest that the cause of this in their simulations is the competing effects of shock acceleration and adiabatic cooling of the particles, with cooling being the dominant factor for shocks that are farther from the Sun (close to 1 AU). The model of Ng et al. (2003) also reveals spectra that have some similarities to the observations; however, there is no long-lasting period where the spectrum is flat in these results (see Figs. 1 and 8 of Ng et al. 2003). In this model, energetic particles excite magnetic fluctuations upstream of the shock, leading to a spatial variation in the turbulent magnetic power spectrum upstream. The source at the shock is assumed, based on standard shock acceleration theory. In each of these models, there are notable differences between the behavior of the time intensity profile in these simulations and the observations shown in Fig. 11. Thus, the cause of the flat-spectra period upstream of the shock remains a poorly understood phenomenon.

The transition from period 1 to period 2 in Fig. 11 also provides an important constraint on acceleration and transport theory. The fact that the intensities rise dramatically over a short (∼ 10-30 minute) period before the shock crossing is not unusual (c.f. Giacalone 2012, and references therein), and is thought to be due to reduced transport in increased magnetic fluctuations near the shock, possibly due to the energetic particles themselves (e.g. Lee 1983; Kennel et al. 1986). However, this particular observation, and several authors as reported by Lario et al. (2018), suggest that not only is there a reduction in the transport coefficient near the shock, but there is also a change in its dependence on energy. These observations, therefore, present important challenges to theory, and to our general understanding of how charged-particles scatter in interplanetary magnetic fields.

For many years now the solar wind TS has been treated as a multifluid phenomenon. The fluids, not separately but jointly, have to fulfil the Rankine-Hugoniot conservation laws (e.g., Fahr and Siewert 2013, 2015, and references therein). It turned out more recently that pick-up ions play a fairly important role for the entropy mediation at the TS, since they absorb, in form of their downstream temperature, a substantial part of the upstream ion kinetic energy. Hence these ions need to be considered as an independent entropy-transporting fluid (e.g., Decker et al. 2008; Fahr and Siewert 2015). Kinetic features of this multifluid ion shock have been studied by, e.g., Zank et al. (2010) and Fahr and Siewert (2013): they found quantitative relations between upstream and downstream ion distribution functions that are substantially different for solar wind protons and pick-up protons.

Surprisingly these refinements did not lead to a completely satisfying representation of the most relevant properties of the downstream plasma as observed by Voyager 2 (Richardson et al. 2008). As demonstrated by Chalov and Fahr (2013) these data can be explained when allowing also the solar wind electrons to thermodynamically behave as a separate, independent fluid. In a simple parameter study, allowing electrons to be heated at the passage over the shock more effectively than protons by a factor \(\xi \gg 1\), they obtained best-fitting results to the relevant Voyager 2 data for \(\xi \geq 10\). Such a specific heating of electrons at fast mode MHD shocks had already earlier been pointed out as a rather likely process by Leroy and Mangeney (1984), Tokar et al. (1986), and Schwartz et al. (1988). Furthermore, a preferential shock-heating of electrons has become manifest in plasma shock simulations in which electrons are treated kinetically (Lembège et al. 2003, 2004). Due to strongly acting two-stream and viscous interactions, they become de-magnetized and do not get Lorentz-wound by the frozen-in, convected magnetic fields, but instead attain downstream temperatures increased by factors of \(\xi \geq 50\) compared to protons.

To achieve a similar level of consistency Fahr et al. (2012) have described the conditions of the upstream and downstream plasma in the bulk frame systems with frozen-in magnetic fields. Describing via the Liouville-Vlasov theorem the shock transition as an instantaneous conversion of the velocity distribution functions, they could determine all relevant downstream plasma quantities. The excessive heating of electrons at the solar wind TS as described by these authors (see also Fahr and Siewert 2013) was obtained with a semi-kinetic theory of the multi-fluid shock with electrons and protons reacting in mass- and charge-specific forms to the electric shock ramp. Again these studies show that electrons do enter the downstream side as a strongly energized, mass-less plasma fluid that even dominates the downstream plasma pressure (for latest results see Fahr and Heyl 2020). These results are also supported by Zieger et al. (2015) through a multi-fluid approach to the solar wind TS in which the plasma electrons are admitted as a separate third fluid. This study gives full support for claims raised by Chalov and Fahr (2013), Fahr and Siewert (2015), and Fahr et al. (2015), namely that electrons are preferentially heated when converting their overshoot velocities into thermal energies, and, w.r.t. their pressure, then play an important role for the downstream plasma flow. Fahr and Verscharen (2016) found that the self-consistent transition structure at the shock in a stationary flow can be described by three consecutive regions, region I with strongly upflashing charge-induced electric fields, region-II with Lorentz reactions of the electrons to the co-convected frozen-in magnetic fields, and region III with a relaxation towards isotropic thermal distributions characterizing the joined downstream plasma flow (see Fig. 12). In the following we explain more carefully, why electrons not only can, but must be expected and respected as a preferentially heated fluid at the passage over the shock describing them by their adequate distribution functions.

Illustration of the concept of a three-region shock. Region I: electric forces dominate over Lorentz forces and arrange the deceleration of the ion plasma flow while enforcing the acceleration of the electrons; Region II: magnetic Lorentz forces dominate, the Buneman instability rapidly converts the electron overshoot kinetic energy into electron thermal energy; Region III: quasi-equilibrium relaxation to the downstream plasma flow (Adapted from Fahr and Verscharen (2016))

Full size image

5.1 The Liouville Theorem and the Electron Distribution Function in Shock-Region I

For arbitrarily, but isotropically distributed particles of a medium, if only co-moving with the plasma bulk velocity \(U\), in absence of particle sources or sinks, the Liouville theorem requires the following form of a conservation of the phase space current (e.g., Goedbloed and Lifschitz 1997; Fahr and Lay 2000)

Here \(f_{1}\) and \(f_{2}\) denote the distribution functions of the species at coordinates \(z_{1},v\) and \(z_{2},v\). Considering now infinitesimally neighbored places on the \(z\)-axis, i.e. with \(z_{2}=z_{1}+dz\), the above relation expresses itself in the more evolved differential form

simply requiring that

The latter requirement can be expressed in a more evolved form by

Now we assume that - while doing a space step \(dz\) - particles in this shock region I exclusively experience the action of the mono-directional shock-electric field \(E=E_{z}\) along the flow direction and find

If applied to electrons co-moving with their genuine bulk velocity \(U_{e}>U_{p}\), one obtains

and analogously for each individual electron with the velocity \(v_{z}\) in the co-moving electron frame. Hence deriving in \(z\) the given velocity space volume \(d^{3}v_{1}=dv_{x}dv_{y}dv_{z}\) one can express

If we now make use of the fact that all electrons co-move with \(U_{e}\), then \(dz/dt=dz^{\prime }/dt=U_{e}\), one can conclude that the velocity space volume \(d^{3}v_{1}\), transported with the bulk flow, appears as an invariant along the \(z\)-axis. The differential equation further simplifies to

Now we use the result of Fahr and Verscharen (2016) that the electric field for the deceleration of the protons (or vice-versa acceleration of the electrons) is given by

we find the following solution of the differential equation in the co-moving electron frame

Adopting upstream conditions with identical bulk velocities of electrons and protons, i.e. \(U_{p1}= U_{e1}=U_{1}\), integrating by parts, and since the deceleration of protons and acceleration of electrons are simply complimentary effects in energy, being related by

the above solution reduces to

Here the function \(f_{1}(z_{1},v)\) is the initial electron distribution in the co-moving upstream frame and may simply be expected to be a Maxwellian \(f_{1}(z_{1},v)=Max(n_{1},U_{1},T_{e1})\), with \(n_{1},U_{1},T_{e1}\) denoting the first three Maxwellian velocity moments, i.e. number density \(n_{e1}=n_{1}\), bulk velocity \(U_{e1}=U_{1}\) and temperature \(T_{e1}\). In case good physical reasons exist to assume upstream electrons as non-Maxwellian, one can simply replace \(f_{1}(z_{1},v)\) by the expression favoured.

The above expression reveals that the temperature of the electron population in region I of the shock due to \((4\pi /3)m\int f_{1}(z_{1},v)v^{2}dv=n_{e}K\) \(T_{e1}\) is conserved, while the number density of the population changes inversely proportional to \(U_{e}\), and the kinetic electron bulk energy changing with \(z\). The same is valid for the shocked protons along analogous arguments, but of course with their genuine bulk velocity, i.e. \(U=U_{p}\ll U_{e}\).

5.2 Electric Currents and Excited Instabilities

As mentioned, the shock transition region can be divided into three different regions I,II,III amongst which the region I is responsible for a basic proper oscillation. As soon as the shock-electric field has self-installed in this region, protons are slowed down and electrons are accelerated, meaning that the local bulk velocities of protons and electrons \(U_{p}\ll U_{e}\) are becoming very different in magnitudes whereby also the quasi-neutrality of the upstream plasma is locally violated, however leading in a time average to a stationary (w.r.t. shock standing) charge structure. Let us assume a linear extent \(\varLambda \) of this electric field structure, then upstream the electron density is \(n_{e1}\), but downstream it only is \(n_{e1}(U_{1}/U_{e})\). In a linear approach this resembles an electric capacitor with a surface charge of \(C=e(n_{e1}- n_{e1}(U_{1}/U_{e})) \simeq en_{e1}\) and a constant electric field \(E\) between the capacitor plates. The linear dimension \(\varLambda \) is related to the restoring polarisation forces with \(eE\varLambda =(1/2)m_{e}U_{e}^{2}\) and can be expressed through its 1D-Debye analogue:

where \(\lambda _{D1}\) as the upstream Debye length, \(M_{1}\simeq 10\) as the upstream sonic Mach number, and the shock compression ratio \(s\simeq 3\). Reminding that the upstream Debye length is calculated by \(\lambda _{D1}= \sqrt{\frac{KT_{1}}{4\pi e^{2}n_{e1}}}\simeq 3000\) cm, then one finds \(\varLambda = 12\) km.

Now one can calculate the typical time period \(\tau _{\varLambda }\) during which an electron is passing over this transition length to change the capacitor charge:

Connected with this period also is a proper oscillation frequency for the transition region which is given by

where \(\omega _{e1}\) is the upstream electron plasma frequency. This hence means that the jet of overshooting electrons is oscillating with typical frequencies of the order \(\omega _{\varLambda }\simeq 45\omega _{e1}\) and hence excites turbulence for shock region II.

5.3 The Electron Distribution Function in Shock Region II

Now we start with the distribution function of the electrons derived for shock region I. In region II this distribution function due to the upcoming dominance of Lorentz forces has to adapt to the magnetic field, frozen-in to the downstream bulk solar wind proton flow, moving there with the plasma bulk velocity \(U_{2}\). Hence the bulk of the overshooting electrons moves with a velocity \(\Delta U=U_{e}-U_{2}\) relative to the bulk of the protons and consequently has velocity components \(\Delta U_{\perp }= \sin \alpha _{2}\cdot \Delta U\) and \(\Delta U_{\parallel }= \cos \alpha _{2}\cdot \Delta U\) perpendicular and parallel to the downstream magnetic field \(\mathbf {B}_{2}\). Here \(\alpha _{2}\) denotes the angle between the downstream frozen-in field \(\mathbf {B}_{2}\) and the velocity \(\Delta \mathbf {U}\). Due to the acting Lorentz forces electrons are forced to first produce a torus distribution with an average velocity \(\Delta U_{\perp }\) and a thermal spread around with a velocity of \(\delta v=\sqrt{3KT_{1}/m_{e}}\). If this toroidal distribution is subject to a fast and effective pitch-angle scattering process ((e.g., Wykes et al. 2001), then an isotropic shell distribution evolves from the initially occurring toroidal distribution given by

This intermediate shell electron distribution function, which is again unstable as shown for the case of pick-up protons by Fahr and Fichtner (2011), will be subject to strong velocity diffusion due to whistler wave turbulence driven by the fluctuating jet of overshooting electrons. The original non-equilibrium shell distribution will asymptotically deform into a quasi-equilibrium distribution, which we assume to be a kappa shell distribution:

which may be considered as the solution of the following kinetic transport equation for the intermediate semi-stationary state arising as the first quasi-stationary case 〈 \(\partial f_{2e}(v)/\partial t\rangle _{t=\tau _{1}}=0\)

5.4 Instabilities Excited in Region II

After the passage through the layer of strong electric field (i.e. region I) the following physical situation is established. The electron bulk velocity, achieved at the end of region I, is higher than both the electron and the proton thermal velocity at its upstream border, i.e. \(U_{e}>c_{e,1} \gg c_{p,1}\) (see also Fig. 12). Hence at a first glance one is confronted with plasma conditions similar to a strong electron current flowing into a background magnetic field and a slow proton bulk flow. This, as is well known, can trigger the Buneman instability (e.g., Alexandrov et al. 2013; Chen 2016), possibly in competition with the modified two-stream instability, (e.g., Matsukiyo and Scholer 2003; Vasko et al. 2018; Egedal et al. 2021) both exciting whistler turbulence.

In view of the above statement, the assumption of a cold plasma is applicable, and the dispersion relation describing the Buneman instability in the reference frame of the proton bulk motion \(U_{p2}\ll U_{e}\) reads (Chashei and Fahr 2013):

Here \(k\) and \(\omega \) denote wave vector and frequency of the excited wave, and \(\theta \) is the inclination between the downstream magnetic field and the velocity \(\boldsymbol{\Delta} \mathbf {U}\). At the entrance to region II we furthermore can assume \(\mathbf{k}\cdot \boldsymbol{\Delta} \mathbf{U}=\mathbf{k}\cdot (\mathbf{U}_{e}-\mathbf{U}_{p})=k_{z}(U_{e}-U_{p})\). Looking now for maximum wave growth at \(\omega \ll k_{z}(U_{e}-U_{p}) \simeq \omega _{p_{e}}\) and keeping in mind that close to the TS one has the following relation between the electron plasma frequency and the gyro frequencies: \(\omega _{pe}\gg \omega _{Be}\gg \omega _{Bp}\), the above relation leads to the constraint

For \(\omega ^{2}\gg \omega _{Bp}^{2}\), one simply obtains

Now one can see that for wavelengths corresponding to \(\omega _{pe}^{2}\geq (k\Delta U\mathbf{)}^{2}\) there exists a wavepower growth (generation of whistler wave turbulence) according to a growth rate

This for \(\sin \theta \simeq 1\) clearly indicates that within typical growth periods \((1/\gamma (k))\) smaller than electron gyroperiods \(\tau _{ge}\)

the accumulated excess kinetic energy of the electrons, gained in form of overshoot velocities in region I, is then very rapidly redistributed, presumably equally branched into turbulent whistler wave energies and thermal energies of the electrons, respectively. This would, however, then mean that the action of the Buneman instability is so strong, that the Lorentz-induced winding of overshooting electrons practically does not occur, but instead the electrons are strongly decelerated by the instability down to the proton bulk velocity \(U_{p}\) (i.e. to \(\Delta U\simeq 0\)), and most probably the electron kinetic energy of the overshoot \(E_{over}=(1/2)m_{e}(U_{e}^{2}-U_{p}^{2})\) is redistributed in equal shares both to the whistler wave power and to the thermal energy of electrons, i.e. electron pressure.

What can be concluded from the above? The aforegoing investigation did show that the two-stream ion-proton distribution, which had developed in region I, as reaction to the strong electric fields of the shock, dissolves in region II, where Lorentz forces in connection with the frozen-in magnetic fields start controlling the ongoing physical motions. However, instead of a Lorentz winding of the motion, the Buneman instability within a period of only a few electron gyroperiods \(\tau _{ge}\) must be expected to remove the kinetic excess energy of the electrons, to take care of a bulk velocity adaptation of protons and electrons, i.e. \(U^{II}=U_{p}^{II}\simeq U_{e}^{II}\), and to arrange a thermalization or heating of the electrons to an energy equivalent of about one half of the original overshoot kinetic energy (i.e. a branching ratio of 0.5 between electron thermal and whistler turbulent energies, respectively). This nevertheless would mean that, after the Buneman instability has stopped operating, an electron temperature \(T_{e}^{II}\) of the order of \(T_{e}^{II}\simeq 0.5(1/4k)MU_{1}^{2}=2.5\cdot 10^{5}~\mbox{K}\) can be expected to prevail, meaning that (comparing this result with the proton temperatures of the order of \(T_{1}^{II}\simeq 1\cdot 10^{5}~\mbox{K}\) measured by Voyager-2), electrons, even after the occurrence of this relaxation process in shock region II, still reflect an effective thermal excess heating w.r.t. protons by at least a factor 2.5. That appears to clearly point out the fact that further downstream of the TS electrons have to be considered as the hotter, more pressurized plasma population compared to the shocked solar wind protons.

As we have shown in previous sections, particle acceleration at MHD shocks is directly observed in the heliosphere. Fast and efficient acceleration of ions at the Earth bow shock detected with Magnetospheric Multiscale mission (MMS) was reported recently in Turner et al. (2018), which is consistent with the models of first-order Fermi acceleration. Turner et al. (2018) identified an efficient form of Fermi acceleration that involves two large-scale converging shocks in some intermittent hot flow anomalies. Electron acceleration by non-resonant compressible fluctuations (magnetic pumping) in the Earth bow shock ripples is shown to be consistent with MMS observations (Lichko and Egedal 2020).

At larger scales in the heliospheric hierarchy the pile-up twin-shock collision scenario (Wang et al. 2019) is successful in modeling observed ground-level enhancements, which were associated with some of the fast coronal mass ejections produced by gradual solar events. Dynamic Monte Carlo simulations were used in Wang et al. (2019) to model the MeV regime ion enhancements in CME events. In the larger scale astrophysical events, like a collision of a supernova shock with stellar wind, the non-linear back-reaction effects may affect CR spectra evolution at colliding shock flows (Bykov et al. 2013b; Vieu et al. 2020). The colliding shock flow models are applied to studies of CR acceleration in the interaction between the strong stellar winds in the colliding-wind binaries (Grimaldo et al. 2019) and in the bow shock pulsar wind nebulae (Bykov et al. 2019a).

Re-acceleration of the fast solar flare electrons by shocks associated with CME events was suggested to explain the observed electron spectral hardening (Petrosian 2016). The combined effect of shocks and turbulent (resonant and non-resonant) acceleration (e.g. (Petrosian and Bykov 2008)) is widely used in solar events.

Particle acceleration at multiple shocks (both identical and with different strengths) and its effect on particle spectral index has been investigated in Pope and Melrose (1994). Recent particle-in-cell simulations of collisionless shock formation and non-thermal particle acceleration are very successful in modeling many astrophysical phenomena (Marcowith et al. 2016). They provided the microscopic description of particle acceleration in the magnetic field reconnection processes (Comisso and Sironi 2019). On the other hand, the limited computational resources do not allow at present to model non-thermal particle distributions with a wide dynamical range in complicated astrophysical multi-scale MHD flows with multiple shocks and the accompanied rarefaction plasma motions. In the extended astrophysical objects, like the clusters of galaxies, the energy containing scale of magnetic turbulence may exceed 10 kpc while the gyroradius of radio-emitting GeV regime electrons is below 10¹⁴ cm. In the case of supernova remnants and in the co-rotation interaction regions of the stellar winds the dynamical range is very wide as well to be covered by the exact PiC models.

Particle acceleration by a strong magnetic turbulence with coherent structures and plasma velocity gradients may need kinetic description outside the quasi-linear theory (e.g. Bykov and Toptygin 1993; Lemoine 2021). In many cases the intensive energy release in active regions is accompanied with formation of shocks of different strengths due to supersonic flow collisions or just non-linear breaking of Riemann type waves. The MHD flow where fast particles are interacting with large amplitude magnetic fluctuations in the presence of the flow discontinuities will produce the intermittent distributions of particles in phase space (Lemoine 2021).

The presence of intermittency requires special description even in statistically homogeneous strong compressible MHD turbulence with shocks comparing to that in the case of the weak turbulence. Apart from the pair correlation functions of the plasma bulk velocity some extra correlation functions describing the shock ensemble are needed to construct the kinetic equations (Bykov and Toptygin 1993). In the test particle case the intermittent distribution function of CRs have strong localized peaks in the vicinity of shock front surfaces. One can use local plain solutions in the vicinity of a shock of a given strength like that given by Eq. (2) to obtain the kinetic equation for \(F(\boldsymbol{r},p,t)\), the average phase space distribution function in between the shocks that is a smooth function in the phase space. The distribution function in Eq. (2) in the vicinity of each of \(i\)-th shock is parameterized by the local diffusion coefficient \(D_{x}^{i}\) and the index \(\gamma _{i}\) which is determined by the shock compression ratio \(r_{i}\). It should be noted here that \(F\cdot p^{2}\) corresponds to \(f_{u}(p)\) in Eq. (2).

The distribution function averaged over an ensemble of turbulent motions and shocks satisfies the kinetic equation

where the source term \(Q(p)\) is determined by injection of the nuclei at some momentum range. The operators \(\hat{L}\) and \(\hat{P}\) for the ensemble of strong shocks (providing the close values of the index \(\gamma \)) are given by

The kinetic coefficients of test particles \(D_{\alpha \beta}\), \(D^{pp}\), \(A,B\) and \(G\) in the case of the strong intermittent shock turbulence satisfy the renormalization equations (Bykov and Toptygin 1993; Bykov et al. 2019b, 2020). The spatial diffusion coefficient \(D_{\alpha \beta}\) describes the particle transport by strong long-wavelength MHD fluctuations in the presence of short scale resonant fluctuations. At low CR energies the transport is dominated by turbulent CR advection with \(D_{\alpha \beta}\) which does not depend on the CR energy. At high energies the CR energy dependent diffusion is determined by the spectrum of MHD turbulence (see e.g. Jokipii 1966; Bykov and Toptygin 1987; Bieber et al. 1994; Casse et al. 2001). The momentum diffusion coefficient \(D^{pp}\) is dominated here by the long-wavelength compressible fluctuations. The terms with \(\hat{L}\) are describing the CR spectra transformations at individual shocks, while the kinetic coefficients \(A,B\) described the cross-correlations of shocks with the MHD compressible turbulence.

Cosmic ray acceleration by large scale plasma motions with the ensemble of MHD shocks are expected to occur in a number of astrophysical objects including clusters of galaxies (Bykov et al. 2019b), clusters of massive stars with energetic stellar winds and supernovae shocks (Bykov et al. 2020), fast radiation-driven stellar winds (White and Chen 1994) and many others. The time dependent solutions of Eq. (43) for CR distribution started from a narrow range of particle momenta injected into the large scale compressible turbulence discussed in Bykov et al. (2019b) to model CRs in clusters of galaxies are approaching the asymptotic power law of the form \(F(p) \propto p^{-5}\). This is rather a soft CR spectrum which may explain the results of the recent analysis of the GeV regime gamma ray data from Fermi LAT on the Coma cluster (Adam et al. 2021). The direct measurements of the characteristics of the magnetic turbulence and cosmic rays are available in the heliosphere where detailed analysis of CR spectra in the quiet periods (i.e. in between of the events accompanied with shocks) revealed the spectra \(F(v) \propto v^{-5}\) (Fisk and Gloeckler 2006, 2014, 2017) which is consistent with that discussed above.

The time dependent particle acceleration model solutions for the averaged distribution function of CRs in between the shocks approached the time asymptotic behaviour \(F(p) \propto p^{-a}\). In Fig. 13 we illustrate the temporal evolution of the spectra of CRs accelerated by multiple shocks. The power-law exponent \(a\) depends on the effective adiabatic index of the system. The index \(a = 5\) for the systems with injection of non-relativistic particles as shown in Fig. 13. However, if the pressure of the accelerated relativistic CR component is comparable to the background plasma pressure as it is the case in the active regions of the ISM such as the superbubbles the power-law exponent is approaching \(a = 4\). This power-law distribution extends up to some characteristic momentum \(p_{c}\) which depends on the small scale MHD turbulence in the system (see Bykov et al. (2020) for details). Such a model can be applied to interpret the existing gamma-ray observations of both large scale superbubbles and compact clusters of young massive stars (Ackermann et al. 2011; Bykov 2014; Aharonian et al. 2019; Bykov et al. 2020).

The temporal evolution of the nonthermal protons accelerated by a multiple shock ensemble. Mono-energetic non-relativistic particles of the momentum \(p=p_{i}\) are injected at \(t=0\). The low energy spectra of the accelerated particles are shown at \(t= 2 \tau _{\mathrm{a}}\) (red curve), \(t = 4 \tau _{\mathrm{a}}\) (green curve) and the asymptotic spectra after \(t= 10 \tau _{\mathrm{a}}\) (blue curve). The acceleration time scale \(\tau _{\mathrm{a}} \sim 9L/U_{s}\) is determined by the energy containing scale \(L\) and the velocity amplitude \(U_{s}\) of the supersonic turbulence

Full size image

In this review we focus on recent developments in the shock acceleration mechanism, also in relation with particle transport. The results reported comes mainly from investigations of particle acceleration in the heliosphere but have strong implications also in distant astrophysical environments. We included evidences of anomalous transport as deduced from the analysis of in-situ data upstream of interplanetary shocks and in numerical simulations. In particular, superdiffusion has been found to give rise to power-law energetic particle fluxes far upstream of the shock when the level of magnetic field fluctuations, at the scales in resonance with the energetic particle Larmor radius, does not depend on the distance from the shock (which suggest a spatially independent diffusion coefficient) (Perri and Zimbardo 2012b). It also reduces the acceleration times and flattens the particle energy spectra with respect to the standard DSA prediction. Further, Monte Carlo models of collisionless shocks with superdiffusive energetic particle transport have addressed how superdiffusion plays an important role in magnetic field amplification in astrophysical shocks. Indeed, superdiffusion upstream of shocks gives rise to a strong particle anisotropy, which resultant pressure generates instabilities in the plasma that amplify the magnetic field. This is a point that is worth investigating using in-situ data in the heliospheric environment, as a natural laboratory for investigation, with the joint analysis of magnetic field, plasma, and energetic particles data. New space missions, such as MMS, Parker Solar Probe and Solar Orbiter will surely advance such a point thanks to the unprecedented high time cadence in plasma and energetic particle data. Also particle anisotropy is a key point that needs to be taken into account when including superdiffusion in shock acceleration.

Since superdiffusion is due to a stochastic process (Lévy walk) that is characterized by a power-law distribution of the free path lengths, the probability of occurrence for long displacements is higher than the one associated to a Gaussian process; very short displacements have an associated high probability, too. Such a scale-free nature of the walk allows particles to make both very long and very short jumps, reducing the time of the particle acceleration process (Perri and Zimbardo 2012a). Further, Zimbardo and Perri (2020) have shown that assuming a power-law distribution for the scattering times, it is possible to derive a fractional Fokker-Planck equation that yields to a superdiffusive spatial transport; anomalous transport can also emerge as solution of a non-linear Fokker-Planck equation in presence of an external force (Tsallis and Bukman 1996).

An efficient process for accelerating CR at high energies is the acceleration at pairs or multiple shocks. Such a mechanism is reported in this review and the governing kinetic equation is discussed, since it finds applicability within a broad range of astrophysical systems. Colliding supersonic flows occur in colliding wind binaries (e.g. Pittard et al. 2021) and can be realized in different astrophysical environments including relativistic pulsar winds or jets produced by black holes in gamma-ray binaries (Bykov et al. 2021). A very efficient first order Fermi-type acceleration can be expected in such of the flows which may form a hard spectral component at the highest energies of CRs accelerated in this sources (Bykov et al. 2013b; Grimaldo et al. 2019; Vieu et al. 2020; Bykov et al. 2021). This type of CR accelerators can produce the high energy CRs component beyond the knee of the galactic CRs spectrum.

One aspect that is fundamental in shock physics is the interaction of propagating shocks with upstream turbulence, since it can play a key role in the evolution of particle energy spectra. As we described in Sect. 3, numerical simulations have shown that the interaction between a propagating shock and turbulence upstream causes a high distortion of the shock front and high spatial and velocity particle diffusion when the amplitude of magnetic field fluctuations upstream is large. Fast particle reflection at the shock can probably be source of superdiffusion. However, further study is required especially from the observational side. Indeed, a possible relation between turbulence properties upstream and particle transport regimes is still poorly understood. Also how shock geometry influences energetic particle transport needs a deeper investigation.

Another observational evidence presented in this review that challenges the standard shock acceleration theory is the presence of flat energy spectra at some interplanetary shock waves. Whether the explanation has to be ascribed to changes in the connection with the satellite and the shock front, or to the interaction between energetic particles and pre-existing turbulent fluctuations that lead to a diffusion coefficient independent of energy remains an open question. Additional insights can of course come from the joint study of observations, analytical and numerical modellings.

Recent studies on particle acceleration at the TS address the importance of the energetic electron population as a major source of plasma heating (Chalov and Fahr 2013). Electrons are found to represent the preferentially heated fluid at MHD shocks, are hotter than shocked protons, and dominate the plasma pressure downstream of the TS. This has been taken into account when interpreting Voyager 2 observations during the TS crossings.