2. Non-diffusive transport in laboratory plasmas

Due to the importance of controlling transport in magnetically confined plasmas, many numerical and theoretical efforts have been done to understand the properties of plasma transport under different configurations and experimental setups (e.g. Dendy et al. Reference Dendy, Chapman and Paczuski2007; Sanchez et al. Reference Sanchez, Newman, Leboeuf, Decyk and Carreras2008; Dewhurst et al. Reference Dewhurst, Hnat and Dendy2010). In recent years, suprathermal ion transport has been extensively studied on the Toroidal Plasma Experiment, TORPEX, an experimental device particularly suited to single out the influence of various parameters in determining the nature of transport.

TORPEX (Fasoli et al. Reference Fasoli, Avino, Bovet, Furno, Gustafson, Jolliet, Loizu, Malinverni, Ricci and Riva2013) is a simple magnetized torus (SMT) with a major radius $R=1~\text{m}$ and a minor radius $a=0.2~\text{m}$ , in which helical, open magnetic field lines are created by combining a dominant toroidal magnetic field $B_{t}\simeq 75~\text{mT}$ with a much weaker vertical magnetic field $B_{v}\simeq 2~\text{mT}$ . The SMT configurations are characterized by the number $N$ of toroidal turns performed by a field line before hitting the vessel. Plasmas of different gases can be produced and sustained by injecting microwaves at 2.45 GHz in the electron cyclotron (EC) frequency range. Typical electron temperature and density are $T_{e}\simeq 1{-}6~\text{eV}$ and $n_{e}\simeq 10^{15}{-}10^{16}~\text{m}^{-3}$ , while the ions are cold ( $T_{i}<1~\text{eV}$ ). TORPEX allows for easy diagnostic access and well-characterized plasma scenarios, thus overcoming the difficulties in direct measurements of suprathermal ion transport. These are limited by the harsh plasma environment in fusion-grade devices and by the difficulty of performing measurements in distant astrophysical plasmas. In addition, detailed knowledge of the turbulence characteristics and of the background plasmas, necessary to realistically model the transport of suprathermal ions, has been achieved in recent years in TORPEX (Ricci et al. Reference Ricci, Theiler, Fasoli, Furno, Gustafson, Iraji and Loizu2011; Gustafson et al. Reference Gustafson, Ricci, Bovet, Furno and Fasoli2012a ,Reference Gustafson, Ricci, Furno and Fasoli b ).

For the present experiments, plasmas are created by a small level of microwave power ( ${\approx}400~\text{W}$ ) on the high-field side of the TORPEX device and with $N\approx 2$ . This scenario is characterized by the presence of ideal interchange modes (with a perpendicular wave number $k_{{\it\Delta}}\simeq 35~\text{rad m}^{-1}$ and a wave number parallel to the magnetic field, $k_{\Vert }\simeq 0$ ). These ideal interchange modes are driven by the magnetic curvature and the pressure gradient and intermittently generate field-aligned plasma structures termed ‘blobs’, which propagate radially outward (Theiler et al. Reference Theiler, Furno, Ricci, Fasoli, Labit, Muller and Plyushchev2009; Furno et al. Reference Furno, Spolaore, Theiler, Vianello, Cavazzana and Fasoli2011).

Figure 1.

Poloidal suprathermal ion current profiles (units in $\text{A m}^{-2}$ ) at different toroidal distances for two ion energies (a,b) $E=70~\text{eV}$ , (c,d) $E=30~\text{eV}$ . The vertical drift due to the curvature and gradient of the magnetic field is shown. The spreading due to the interaction with the turbulent plasmas is more important for 30 eV ions. (e) Radial width of suprathermal ion current profiles as a function of toroidal distance. Red squares and blue circles represent experimental measurements for ions emitted at 30 eV and 70 eV, respectively. Continuous bands are obtained from a synthetic diagnostic using numerical simulations for 30 eV (red) and 70 eV (blue) ions. The width of the bands is obtained by varying the simulation input parameters within experimental uncertainties. The different beam spreadings as a function of the toroidal distance indicate different transport regimes.

In these plasmas, suprathermal lithium $6^{+}$ ions are injected in the blob region using a miniaturized ion source and detected using gridded energy analysers (GEAs), which were specifically developed for these experiments (Bovet et al. Reference Bovet, Furno, Fasoli, Gustafson and Ricci2012, Reference Bovet, Furno, Fasoli, Gustafson and Ricci2013, Reference Bovet, Fasoli and Furno2014a ,Reference Bovet, Gamarino, Furno, Ricci, Fasoli, Gustafson, Newman and Sánchez b ). The suprathermal ion source is mounted on a moving system, which can continuously position the source at different toroidal locations. Two GEA detectors moving across almost the entire cross-section are installed at different toroidal distances along the suprathermal ion beam. This setup allows measurements of the three-dimensional properties of the suprathermal ion beam as it interacts with the plasma turbulence. For these measurements, two detection schemes can be used: (1) time-resolved measurements allowing for the evaluation of the statistical moments of the suprathermal ion current fluctuations and (2) synchronous detection to improve the signal to noise ratio, which allows for the measurement of the time-averaged profiles.

In figure 1, we show the spreading of the time-averaged suprathermal ion current profiles in the presence of plasma, measured at different toroidal locations for two different energies, $E=30~\text{eV}$ (figure 1 a,b) and $E=70~\text{eV}$ (figure 1 c,d), using the synchronous detection technique. The spreading of the suprathermal ion beam is quantified by the radial variance of the two-dimensional poloidal profile as a function of the toroidal distance, as shown in figure 1(e). In the present simple magnetized torus regime in the absence of a plasma, the motion of suprathermal ions corresponds to the gyromotion along the magnetic field lines at the cyclotron frequency, $f_{i}\simeq 188~\text{kHz}$ , which is superposed to an upward drift due to the curvature and gradient of the magnetic field. Figure 1 shows that on top of the unperturbed motion, a broadening due to the interaction with the plasma turbulence is present, which is different for the two suprathermal ion energies, hinting at different transport regimes. We note that, although the three-dimensional suprathermal ion beam profile is directly measured in TORPEX, in order to identify the transport regime the radial transport exponent, ${\it\alpha}_{R}$ , should be computed from the temporal dependence of the variance of the displacement of the suprathermal ions. The experimental data of figure 1 are time averaged and the radial spreading of the ions is only accessible as a function of the toroidal distance.

To gain insight into the nature of transport, the experimental measurements are compared and analysed using results from numerical simulations with the help of a synthetic diagnostic reproducing the suprathermal ion current profile measured with the detectors. A large number of tracers ( $1.6\times 10^{5}$ ) are injected in global fluid turbulent simulations that were previously validated against experimental data. The tracer trajectories are computed using the Newton equation of motion. Two simulations reproducing the experimental results are shown in figure 1(e) and are used to draw the bands in this figure, which quantify the simulation uncertainties.

Using the simulation results reproducing the experimental conditions, the value of the radial transport exponent, ${\it\alpha}_{R}$ , can be readily computed from the evolution of the variance of the ion radial displacements as a function of time, ${\it\sigma}_{R}^{2}(t)$ . At first, the transport of the ions is ballistic ( ${\it\sigma}_{R}^{2}(t)\propto t^{2}$ ), during a phase in which the ions have not yet interacted with the plasma and are not yet magnetized. The ions enter a second spreading phase as they start to interact with the plasma turbulence. In this phase, different transport regimes are observed, according to the energy of the ions and the character of the turbulence. A numerical study showed that in this interaction phase the transport can vary from a subdiffusive to a superdiffusive regime depending on two parameters that determine the relative sizes of the ion orbits and the turbulent structures: the injection energy normalized to the electron temperature, $E/T_{e}$ , and the normalized fluctuations amplitude, $e\tilde{{\it\phi}}/T_{e}$ (Gustafson et al. Reference Gustafson, Ricci, Furno and Fasoli2012b ). Here $\tilde{{\it\phi}}$ is the fluctuating electric potential due to the plasma turbulence. Fitting the temporal evolution of ${\it\sigma}_{R}^{2}$ to power laws provides the values of the transport exponents in the different phases.

In the interaction phase, an exponent ${\it\alpha}_{R}=0.51\pm 0.01$ is found for ions of 70 eV ( $E/T_{e}\simeq 54$ ) and ${\it\alpha}_{R}=1.20\pm 0.04$ for ions of 30 eV ( $E/T_{e}\simeq 23$ ), indicating that the transport varies from subdiffusive to superdiffusive as the energy of the ions is decreased. For ions of 30 eV, after the superdiffusive phase, a phase where the transport is close to diffusive ( ${\it\alpha}_{R}=0.92\pm 0.04$ ) is visible in figure 1(e) after approximately 1 m. This phase appears when the size of the beam becomes sufficiently large that ions sample regions of the plasma with a different fluctuation amplitude, originating an average transport close to diffusive (Bovet et al. Reference Bovet, Fasoli and Furno2014a ,Reference Bovet, Gamarino, Furno, Ricci, Fasoli, Gustafson, Newman and Sánchez b , Reference Bovet, Fasoli, Ricci, Furno and Gustafson2015).

In order to relate the above results to theoretical studies, we can notice that there is a formal analogy between the transport of particles in the presence of electric fluctuations, if they obey the drift approximation, and the transport of magnetic field lines perpendicular to the magnetic field, in the case of two-dimensional fluctuations. Actually, in both cases a Hamiltonian structure for the field line equations or the particle motion equations can be found. In the former case, the electric drift velocity can be written as $\boldsymbol{v}_{E}=c\boldsymbol{E}\times \boldsymbol{B}/B^{2}$ . Considering $xy$ as the plane perpendicular to $\boldsymbol{B}$ , we have

(2.1a,b ) $$\begin{eqnarray}\frac{\text{d}x}{\text{d}t}=-\frac{c}{B}\frac{\partial \tilde{{\it\phi}}}{\partial y},\quad \frac{\text{d}y}{\text{d}t}=\frac{c}{B}\frac{\partial \tilde{{\it\phi}}}{\partial x},\end{eqnarray}$$

where $x$ corresponds to the coordinate, $y$ to the impulse, and the Hamiltonian is given by $H=-c\tilde{{\it\phi}}/B$ . Given that $\tilde{{\it\phi}}$ and $B$ depend only weakly on the coordinate along $\boldsymbol{B}$ , in the limit of small gyroradius the particles follow the contours of $\tilde{{\it\phi}}$ . In a similar way, for magnetic field line transport we have $\text{d}x/\text{d}z=B_{x}/B_{z}$ and $\text{d}y/\text{d}z=B_{y}/B_{z}$ , which, in terms of the vector potential $\boldsymbol{A}(x,y)$ , assuming two-dimensional magnetic fluctuations, yields

(2.2a,b ) $$\begin{eqnarray}\frac{\text{d}x}{\text{d}z}=\frac{1}{B}\frac{\partial A_{z}}{\partial y},\quad \frac{\text{d}y}{\text{d}z}=-\frac{1}{B}\frac{\partial A_{z}}{\partial x},\end{eqnarray}$$

with $z$ playing the role of time and the Hamiltonian given by $H=A_{z}/B$ . Hence, magnetic field lines follow the contours of the vector potential.

This allows for a possible interpretation of TORPEX experimental results in terms of previous numerical studies of field line random walk, with the transport in the radial direction corresponding to transport in the plane perpendicular to the magnetic field. Several numerical studies show that anomalous transport (both subdiffusion and superdiffusion) is found for a small Kubo number  $R$ . The Kubo number is given, for magnetic field lines, by the product of the turbulence level ${\it\delta}B/B_{0}$ times the ratio of parallel and perpendicular correlation lengths ${\it\lambda}_{\Vert }/{\it\lambda}_{\bot }$ (e.g. Isichenko Reference Isichenko1992). For the $\boldsymbol{E}\times \boldsymbol{B}$ drift, the Kubo number is given by the drift speed $c{\it\delta}E/B_{0}$ times the ratio of the correlation time of ${\it\delta}E$ (with ${\it\delta}E$ the fluctuating electric field amplitude, which is related to the fluctuating electric potential by ${\it\delta}\boldsymbol{E}=-\boldsymbol{{\rm\nabla}}\tilde{{\it\phi}}$ ) over the perpendicular correlation length ${\it\tau}_{corr}/{\it\lambda}_{\bot }$ . What is found is that for $R<1$ we have: (i) anomalous transport, both for field lines and for particles (Zimbardo et al. Reference Zimbardo, Veltri and Pommois2000b ; Pommois et al. Reference Pommois, Zimbardo and Veltri2007); (ii) the structure of magnetic flux tubes becomes distorted but does not look Gaussian, rather it looks convoluted and fractal. When this structure is crossed by spacecraft (because of advection by the solar wind) the so-called solar energetic particle dropouts are observed: these are due to the spacecraft seeing both regions filled with SEPs, magnetically connected with the energetic particle source, and empty regions, which are not magnetically connected to the source; some events also show a fine structure of filled/empty regions (Giacalone et al. Reference Giacalone, Jokipii and Mazur2000; Mazur et al. Reference Mazur, Mason, Dwyer, Giacalone, Jokipii and Stone2000; Zimbardo et al. Reference Zimbardo, Pommois and Veltri2004; Trenchi et al. Reference Trenchi, Bruno, Telloni, D’Amicis, Marcucci, Zurbuchen and Weberg2013). On the contrary, for $R\gtrsim 1$ we have mostly normal diffusion, a nearly Gaussian structure for the magnetic flux tubes (sometimes very dense), and the trajectory of a random walker has fractal dimension equal to 2 (for time going to infinity) (Isichenko Reference Isichenko1992).

In addition, when superdiffusion is found, one can argue that the spatial structure of the random walker trajectory is fractal, and the fractal dimension can be related to the superdiffusion exponent (Bouchaud & Georges Reference Bouchaud and Georges1990; Isichenko Reference Isichenko1992). If superdiffusion is due to a Lévy process described by the free path probability in (1.2), it is found that the fractal dimension of the trajectory in a plane is given by $D_{F}={\it\mu}$ . If one takes a one-dimensional cut of the trajectory in the plane perpendicular to the magnetic field, its fractal dimension is $D_{F}^{cut}=D_{F}-1={\it\mu}-1$ . We can try to relate the fractal dimension of the distribution of superdiffusive ions to the observed superdiffusion exponent ${\it\alpha}_{R}\simeq 1.24$ (Bovet et al. Reference Bovet, Fasoli and Furno2014a ), so that $D_{F}={\it\mu}=2/{\it\alpha}\simeq 1.61$ . Crossing such a structure, as happens in TORPEX because of the transverse displacements of the plasma column due to turbulence, can give an intermittent signal as the one found for lithium ions in the case of energy equal to 30 eV (Bovet et al. Reference Bovet, Gamarino, Furno, Ricci, Fasoli, Gustafson, Newman and Sánchez2014b ), with a predicted fractal dimension in time of order $D_{F}^{cut}={\it\mu}-1\simeq 0.6$ .

In summary, the intermittency in the ion current which is found in TORPEX in the case of superdiffusion could be interpreted on the basis of the fact that for Kubo numbers of order one or less, anomalous transport perpendicular to the magnetic field is possible, and the corresponding structure of the magnetic flux tubes which are populated by lithium ions becomes convoluted and fractal, corresponding to the properties (i) and (ii) above. Indeed, the Kubo number for TORPEX experiments is estimated to be of order one. Therefore, this suggests that the anomalous behaviour is due to the structures induced by electric field turbulence, in agreement with the fact that superdiffusion and the intermittency of the current are not found for higher energy ions, which, having a larger gyroradius, are less influenced by the spatial structure of  ${\it\delta}\boldsymbol{E}$ .

3. Superdiffusive transport at heliospheric shocks

In numerical simulations of transport as well as in the particular case of lithium ions in the TORPEX device, one can trace the particle trajectories and hence directly obtain the mean square displacement and the corresponding transport properties. However, this is hardly possible in space plasmas, and certainly not in astrophysical plasmas. Therefore, for many systems the transport properties have to be obtained indirectly, considering the effects that transport has on some observable quantities. Space plasmas are in an intermediate situation, since a number of plasma properties can be measured in situ by spacecraft. In this context, Perri & Zimbardo (Reference Perri and Zimbardo2007, Reference Perri and Zimbardo2008) and Perri et al. (Reference Perri, Zimbardo, Effenberger and Fichtner2015) have developed some diagnostic tools to extract the transport properties of energetic particles accelerated at heliospheric shocks by the measured intensity profiles of the energetic particles.

The method is based on the fact that in the case of superdiffusive transport described by a Lévy walk, the probability density function of positions, i.e. the propagator $P(x,t)$ , is not a Gaussian but rather a probability distribution with power law tails, similar to a Lévy distribution. This can be obtained by solving the Montroll–Weiss equation in Fourier–Laplace space (e.g. Zumofen & Klafter Reference Zumofen and Klafter1993): for small values of the scaling variable ${\it\xi}=(x/\ell _{0})/(t/t_{0})^{1/{\it\mu}}\ll 1$ one has a modified Gaussian,

(3.1) $$\begin{eqnarray}P(x,t)\simeq \frac{{\it\Gamma}(({\it\mu}+1)/{\it\mu})}{{\rm\pi}(Ct)^{1/{\it\mu}}}\exp \left[-\frac{{\it\Gamma}(3/{\it\mu})}{2{\it\Gamma}(1/{\it\mu})}\left[\frac{x}{(Ct)^{1/{\it\mu}}}\right]^{2}\right],\end{eqnarray}$$

where ${\it\Gamma}$ is the Euler gamma function and the constant $C$ is given by

(3.2) $$\begin{eqnarray}C=2\frac{{\it\mu}-1}{{\it\mu}+1}\left|\cos \left(\frac{{\rm\pi}}{2}{\it\mu}\right)\right|{\it\Gamma}(-{\it\mu})\frac{\ell _{0}^{{\it\mu}}}{t_{0}}\end{eqnarray}$$

(Zumofen & Klafter Reference Zumofen and Klafter1993; Zimbardo & Perri Reference Zimbardo and Perri2013; Perri et al. Reference Perri, Zimbardo, Effenberger and Fichtner2015). In the definition of ${\it\xi}$ , $\ell _{0}$ is the shortest length for which the free paths have a power law probability distribution, see (1.2), and $t_{0}=\ell _{0}/v$ is the corresponding time. Conversely, for large distances, ${\it\xi}\gg 1$ but $|x|<vt$ , one obtains a power law propagator,

(3.3) $$\begin{eqnarray}P(x,t)\simeq \frac{{\it\Gamma}({\it\mu}+1)}{{\rm\pi}}\sin \left[\frac{{\rm\pi}}{2}{\it\mu}\right]|C|\frac{t}{|x|^{{\it\mu}+1}},\end{eqnarray}$$

with the propagator going to zero for $x>vt$ (Blumen et al. Reference Blumen, Klafter and Zumofen1990; Zumofen & Klafter Reference Zumofen and Klafter1993). As can be seen, the long-distance propagator has a power law form, basically different from the normal Gaussian propagator. Other forms of non-Gaussian propagators in the astrophysical context have been discussed by Ragot & Kirk (Reference Ragot and Kirk1997) and by Webb et al. (Reference Webb, Zank, Kaghashvili and le Roux2006).

The propagator allows us to compute the density of particles of a given energy at any position in space and time as (e.g. Ragot & Kirk Reference Ragot and Kirk1997)

(3.4) $$\begin{eqnarray}n(x,t)=\int _{-\infty }^{\infty }\,\text{d}x^{\prime }\int _{-\infty }^{t}\,\text{d}t^{\prime }P(x-x^{\prime },t-t^{\prime })Q_{sh}(x^{\prime },t^{\prime }),\end{eqnarray}$$

where $Q_{sh}={\it\Phi}_{0}{\it\delta}(x^{\prime }-V_{sh}t^{\prime })$ indicates the particle source, located at the shock which is moving with speed $V_{sh}$ . Here ${\it\Phi}_{0}$ is the flux of energetic particles, assumed to be injected at the shock. When using the propagator for Lévy walks, (3.4) allows us to obtain some of the most important modifications of the standard scenario due to superdiffusive transport. First, we can compute the particle density far upstream of the shock, where the form (3.3) of the propagator applies, which yields

(3.5) $$\begin{eqnarray}n(x,t)\simeq \frac{{\it\Gamma}({\it\mu}+1)}{{\rm\pi}}\frac{\sin ({\rm\pi}{\it\mu}/2)|C|{\it\Phi}_{0}}{{\it\mu}({\it\mu}-1)}\frac{(x-V_{sh}t)^{1-{\it\mu}}}{V_{sh}^{2}}\end{eqnarray}$$

(Perri & Zimbardo Reference Perri and Zimbardo2007, Reference Perri and Zimbardo2008). We can set the observer time at $t=0$ (like in a snapshot), so that far off the shock the spatial profile of the accelerated particles is a power law decay with slope $a={\it\mu}-1$ . This is markedly different from the predictions of normal diffusion, which yields an exponential profile upstream of the shock (e.g. Lee & Fisk Reference Lee and Fisk1982). On the other hand, considering the varying level of Alfvén waves generated by cosmic rays streaming away from the shock, Bell (Reference Bell1978) obtained an energetic particle profile which is a power law with slope $a=1$ . However, superdiffusive transport yields a power law profile even in the case of constant turbulence level, constant turbulence anisotropy (Perri et al. Reference Perri, Yordanova, Carbone, Veltri, Sorriso-Valvo, Bruno and Andre2009; Perri & Balogh Reference Perri and Balogh2010), and constant transport properties (Perri & Zimbardo Reference Perri and Zimbardo2012b ): accordingly, an energetic particle power law profile with $0<a<1$ implies superdiffusive transport with ${\it\alpha}=3-{\it\mu}=2-a$ . In figure 2, we show an example of energetic particle fluxes measured by the Voyager 2 spacecraft at the termination shock, the first boundary of the heliosphere. The obtained values of the slope are $a=0.69{-}0.71$ , so that superdiffusion with ${\it\alpha}\simeq 1.3$ is deduced for ions with energies from 200 keV to 3.5 MeV (Perri & Zimbardo Reference Perri and Zimbardo2009b ).

Figure 2.

Energetic ion fluxes for several energy channels (as indicated) measured upstream of the termination shock, as a function of ${\rm\Delta}t=|t-t_{shock}|$ . The fits to the power law profiles are shown by dashed lines. Adapted from Perri & Zimbardo (Reference Perri and Zimbardo2012b ).

It is interesting to note that similar results for the power law density profile upstream of the shock can be obtained using a fractional diffusion–advection equation of the form (e.g. Stern et al. Reference Stern, Effenberger, Fichtner and Schäfer2014)

(3.6) $$\begin{eqnarray}\frac{\partial n}{\partial t}={\it\kappa}\frac{\partial ^{{\it\mu}}n}{\partial |x|^{{\it\mu}}}+u\frac{\partial n}{\partial x}+{\it\delta}(x),\end{eqnarray}$$

where ${\it\kappa}$ is the superdiffusion coefficient, with dimensions $\text{length}^{{\it\mu}}$ /time, $1<{\it\mu}<2$ is the fractional exponent for superdiffusion, and $u$ is the background advection speed. The fractional spatial derivative is given by the Riesz derivative, defined as

(3.7) $$\begin{eqnarray}\displaystyle \frac{\partial ^{{\it\mu}}n(x,t)}{\partial |x|^{{\it\mu}}}=\frac{1}{{\rm\pi}}\sin \left(\frac{{\rm\pi}}{2}{\it\mu}\right){\it\Gamma}(1+{\it\mu})\int _{0}^{\infty }\frac{n(x+{\it\xi})-2n(x)+n(x-{\it\xi})}{{\it\xi}^{1+{\it\mu}}}\,\text{d}{\it\xi} & & \displaystyle\end{eqnarray}$$

(Samko et al. Reference Samko, Kilbas and Marichev1993; Saichev & Zaslavsky Reference Saichev and Zaslavsky1997). This can be regarded as a fractional generalization of the usual Laplace operator, so that (3.6) is the fractional analogue of the standard diffusion–advection equation. The latter is recovered in the limit ${\it\mu}=2$ , and gives the above-mentioned well-known exponential upstream density profile for shock-injected particles. Asymptotic expressions for solutions to the fractional diffusion–advection equation are discussed in Litvinenko & Effenberger (Reference Litvinenko and Effenberger2014). We only give a short summary of these results here.

An important observation is that the asymptotic form of the propagator obtained from (3.6), without advection and a delta source in space and time, is identical to the leading order power law form of the Lévy walk propagator, i.e. (3.3) upon identifying the constant $|C|$ with  ${\it\kappa}$ . The solution to (3.6) is given by

(3.8) $$\begin{eqnarray}n(x,t)=\int _{0}^{t}P(x+ut^{\prime },t^{\prime })\,\text{d}t^{\prime }.\end{eqnarray}$$

From this relation, two limiting cases for $n(x,t)$ with $x\gg ut$ and $0<x\ll ut$ can be derived. By changing the reference frame, the power law expression for the upstream particle density, as given by (3.5), can be recovered from these results in the limit of a very distant initial shock position. The actual approach of the solution to such a steady-state situation is an interesting field of study that is accessible now with a formulation in terms of a fractional transport equation. Litvinenko & Effenberger (Reference Litvinenko and Effenberger2014) verified their results with complementary solution methods using a weak diffusion approximation and a Fourier-series solution on a finite domain (Stern et al. Reference Stern, Effenberger, Fichtner and Schäfer2014). Furthermore, the relation ${\it\alpha}=3-{\it\mu}$ is recovered for the time dependence of the variance of particle displacement when a finite particle speed is considered (see also the discussion in Perri et al. Reference Perri, Zimbardo, Effenberger and Fichtner2015).

An additional, complementary solution method for fractional diffusion equations is based on a generalization of the equivalence between stochastic differential equations (SDEs) and the Fokker–Planck equation (Gardiner Reference Gardiner2009). The Wiener process representing the stochastic driver in the usual SDE formulation is generalized to Lévy motions $L_{{\it\mu}}(t)$ obeying the Fourier transform characteristic

(3.9) $$\begin{eqnarray}{\mathcal{F}}\{\text{e}^{\text{i}kL_{{\it\mu}}(t)}\}=\text{e}^{-t|k|^{{\it\mu}}}.\end{eqnarray}$$

Identifying again ${\it\kappa}$ and $|C|$ , the Fourier inversion of this expression yields the limiting forms of the propagator given in (3.1) and (3.3). The reader is referred to Magdziarz & Weron (Reference Magdziarz and Weron2007), Bovet et al. (Reference Bovet, Fasoli and Furno2014a ), and Effenberger (Reference Effenberger, Pogorelov, Audit and Zank2014) for more details on solution methods and applications of this approach.

4. Superdiffusive shock acceleration and hard electron spectra at SNR shocks

An important implication of anomalous diffusion is that particle acceleration mechanisms based on diffusive, stochastic motion are modified because the underlying statistical process is not Gaussian. In particular, the cosmic ray energization mechanism is thought to be the diffusive shock acceleration (DSA), a first-order Fermi process which works both for ions and electrons (e.g. Bell Reference Bell1978; Drury Reference Drury1983; Perri et al. Reference Perri, Zimbardo and Greco2011; Balogh et al. Reference Balogh, Bykov, Lin, Raymond and Scholer2013; Giacalone Reference Giacalone2013; Amato Reference Amato2014). It was shown by Duffy et al. (Reference Duffy, Kirk, Gallant and Dendy1995) and by Kirk et al. (Reference Kirk, Duffy and Gallant1996) that in the case of subdiffusion the acceleration time and the energy spectral index are modified, and an application to the diffuse radio emission observed from the Coma cluster of galaxies was proposed (Ragot & Kirk Reference Ragot and Kirk1997). The extension to the superdiffusive case was carried out by Perri & Zimbardo (Reference Perri and Zimbardo2012a ) and by Zimbardo & Perri (Reference Zimbardo and Perri2013), and similar modifications to DSA in the case of Richardson diffusion were recently considered by Lazarian & Yan (Reference Lazarian and Yan2014) (see below). Here we propose for the first time how application of superdiffusion to electrons accelerated at supernova remnants (SNRs) can help to explain the observed radio synchrotron radiation spectra.

Supernova remnant shocks are the most popular candidate for the acceleration of galactic cosmic rays (CRs) (e.g. Helder et al. Reference Helder, Vink, Bykov, Ohira, Raymond and Terrier2012; Blasi Reference Blasi2013; Amato Reference Amato2014). This association has gained important observational support in recent years thanks to the data collected by X-ray telescopes such as Chandra and XMM-Newton and gamma-ray telescopes such as Fermi and AGILE, with the former instruments showing accelerated electrons with energies up to tens of TeV (Reynolds et al. Reference Reynolds, Gaensler and Bocchino2012; Vink Reference Vink2012) and the latter highlighting for the first time the presence of relativistic protons in SNRs (Abdo et al. Reference Abdo, Ackermann, Ajello, Baldini, Ballet, Barbiellini, Baring, Bastieri, Baughman and Bechtol2010a ,Reference Abdo, Ackermann, Ajello, Baldini, Ballet, Barbiellini, Bastieri, Baughman, Bechtol and Bellazzini b ; Tavani et al. Reference Tavani, Giuliani, Chen, Argan, Barbiellini, Bulgarelli, Caraveo, Cattaneo, Cocco and Contessi2010; Giuliani et al. Reference Giuliani, Cardillo, Tavani, Fukui, Yoshiike, Torii, Dubner, Castelletti, Barbiellini and Bulgarelli2011). However, a number of issues on the actual acceleration processes remain open.

Figure 3.

Distribution of synchrotron spectral indices ${\it\alpha}_{{\it\nu}}$ for shell-type SNRs. Adapted from Reynolds et al. (Reference Reynolds, Gaensler and Bocchino2012).

Information on the electron acceleration can be obtained from synchrotron emission at shell-type SNRs (Green Reference Green2009; Morlino et al. Reference Morlino, Amato, Blasi and Caprioli2010; Reynolds et al. Reference Reynolds, Gaensler and Bocchino2012; Vink Reference Vink2012). For these SNRs one can confidently assume that electron acceleration happens at the blast wave of the supernova ejecta. The synchrotron radiation of several SNRs is observed both in the X-rays and in the radio bands. The latter radiation is due to electrons with energies in the range between a few hundred MeV and a few GeV, and it is characterized by a spectral density $S_{{\it\nu}}$ at frequency ${\it\nu}$ :

(4.1) $$\begin{eqnarray}S_{{\it\nu}}\propto {\it\nu}^{-{\it\alpha}_{{\it\nu}}}\end{eqnarray}$$

(Reynolds et al. Reference Reynolds, Gaensler and Bocchino2012), where ${\it\alpha}_{{\it\nu}}$ is the power law spectral index in the radio band. This is related to the spectral index ${\it\gamma}$ of the electron differential energy distribution by (e.g. Longair Reference Longair1994)

(4.2) $$\begin{eqnarray}{\it\gamma}=2{\it\alpha}_{{\it\nu}}+1,\end{eqnarray}$$

so that observing the radio synchrotron emission gives information on the electron energy distribution. On the other hand, the well-known relation between the energy spectral index for ultrarelativistic particles and the shock compression ratio $r$ given by DSA is

(4.3) $$\begin{eqnarray}{\it\gamma}=\frac{r+2}{r-1},\end{eqnarray}$$

so that the compression ratio might be deduced from the observation of ${\it\alpha}_{{\it\nu}}$ , if we restrict our reasoning to the test particle approach of DSA.

A puzzling property of the observed radio synchrotron spectral indices (Onić Reference Onić2013), as obtained for instance from Green’s catalogue (Green Reference Green2009), is that the distribution of ${\it\alpha}_{{\it\nu}}$ is peaked at 0.5, implying ${\it\gamma}=2$ and $r=4$ , but values of ${\it\alpha}_{{\it\nu}}$ substantially smaller than 0.5 are also found, see figure 3. These imply ${\it\gamma}<2$ , which, according to (4.3), would require a compression ratio larger than 4. However, for non-relativistic shock speeds, and if no distinction is made between the gas compression ratio and that due to the motion of the scattering centres (e.g. Vainio & Schlickeiser Reference Vainio and Schlickeiser1999), $r=4$ is the largest compression ratio that can be obtained within the test-particle regime. We note that for the blast wave of SNRs, the shock speeds are in the range $10^{3}{-}10^{4}~\text{km s}^{-1}$ , and the sound speed in the interstellar medium is of the order of a few tens of $\text{km s}^{-1}$ , so that the appropriate thermodynamical description is that of a strong non-relativistic shock, namely the Mach number $Ma\gg 1$ and the compression ratio $r\approx 4$ , unless efficient particle acceleration occurs.

In this latter case the total effective compression ratio can be increased by the fact that accelerated particles take away energy from the system, making the shock effectively radiative. Indeed, nonlinear shock acceleration theory, a framework which takes into account the backreaction of accelerated particles on the shock properties, predicts the formation of an extended cosmic ray precursor ahead of the fluid discontinuity, which now takes the name of subshock (e.g. Drury Reference Drury1983; Amato & Blasi Reference Amato and Blasi2005; Bykov et al. Reference Bykov, Ellison and Renaud2012; Reynolds et al. Reference Reynolds, Gaensler and Bocchino2012): the overall compression ratio from far upstream to downstream becomes larger than what the Rankine–Hugoniot conditions would lead us to infer, while the compression at the fluid subshock becomes less than the expected value. Such a precursor is actually observed for the solar wind termination shock (TS) crossed by the Voyager 2 spacecraft (Decker et al. Reference Decker, Krimigis, Roelof, Hill, Armstrong, Gloeckler, Hamilton and Lanzerotti2008; Richardson et al. Reference Richardson, Kasper, Wang, Belcher and Lazarus2008). Indeed, for the TS the energetic particle precursor allows us to go from the compression ratio of $r\simeq 2$ observed just at the shock (Decker et al. Reference Decker, Krimigis, Roelof, Hill, Armstrong, Gloeckler, Hamilton and Lanzerotti2008; Richardson et al. Reference Richardson, Kasper, Wang, Belcher and Lazarus2008) to a global compression ratio of $r\simeq 2.4{-}2.8$ (Florinski et al. Reference Florinski, Decker, le Roux and Zank2009; Arthur & le Roux Reference Arthur and le Roux2013). For SNRs, the extended precursor with increased compression ratio, perhaps as large as $r=7{-}8$ , leads to a spectral hardening for higher energy particles which sample a larger region around the shock (Reynolds et al. Reference Reynolds, Gaensler and Bocchino2012; Amato Reference Amato2014). However, this is not the case for the GeV electrons emitting synchrotron radiation at radio wavelengths. Indeed, because of their limited energy and Larmor radius, they are expected not to diffuse far in the precursor, but rather stay close to the ‘fluid subshock’, and only experience the compression ratio there. Various models give a compression ratio at the subshock in the range $r=3{-}3.8$ (Reynolds Reference Reynolds2008; Reynolds et al. Reference Reynolds, Gaensler and Bocchino2012), corresponding to ${\it\gamma}=2.07{-}2.5$ and to ${\it\alpha}_{{\it\nu}}=0.53{-}0.74$ . Therefore, the observed values of ${\it\alpha}_{{\it\nu}}>0.5$ can be explained by DSA with a compression ratio less than 4, but values of ${\it\alpha}_{{\it\nu}}<0.5$ (which are quite a few; see figure 3) cannot be interpreted within the same scheme. A somewhat similar problem is found for heliospheric shocks (Lee et al. Reference Lee, Mewaldt and Giacalone2012).

Here we propose that superdiffusive shock acceleration (SSA) (Perri & Zimbardo Reference Perri and Zimbardo2012a ) can explain those hard spectral indices. Indeed, a lower particle density far downstream leads to a modified probability of escape from the acceleration region in DSA (Kirk et al. Reference Kirk, Duffy and Gallant1996). A simple computation based on (3.4) shows that, instead of (4.3), SSA yields a new differential energy spectral index for shock accelerated particles,

(4.4) $$\begin{eqnarray}{\it\gamma}=\frac{6}{r-1}\frac{2-{\it\alpha}}{3-{\it\alpha}}+1,\end{eqnarray}$$

where $r$ is the shock compression ratio, ${\it\alpha}$ is the superdiffusion exponent, and where a corresponding expression can be obtained for non-relativistic particles (Perri & Zimbardo Reference Perri and Zimbardo2012a ; Zimbardo & Perri Reference Zimbardo and Perri2013). We can say that the spectral index quantifies the statistical distribution of cosmic ray energies: since superdiffusion is due to a process, the Lévy walk, characterized by a statistics different from the Gaussian one, see (1.2), (3.1), and (3.3), the distribution of energies is also different, even if it is still a power law. Note that ${\it\gamma}$ depends both on the shock compression ratio and on the exponent of superdiffusion ${\it\alpha}$ , so we can obtain a wide range of values for ${\it\gamma}$ . In particular, for ${\it\alpha}>1$ , harder spectral indices than those predicted by DSA are obtained. The above expression of ${\it\gamma}$ has been applied to the interpretation of energetic ion spectra at the termination shock by Perri & Zimbardo (Reference Perri and Zimbardo2012a ) and to electron spectra at interplanetary shocks by Perri et al. (Reference Perri, Zimbardo, Effenberger and Fichtner2015).

Figure 4.

Superdiffusion exponent ${\it\alpha}$ versus the inferred electron spectral index ${\it\gamma}$ , obtained from (4.5) for $r=4$ . The exact values of ${\it\alpha}$ are indicated close to each data point.

For SNRs, assuming that we can obtain ${\it\gamma}$ from the observed spectral index ${\it\alpha}_{{\it\nu}}$ for radio synchrotron emission, (4.2), and the compression ratio from modelling or other physical arguments, we can infer the superdiffusion exponent as

(4.5) $$\begin{eqnarray}{\it\alpha}=3\frac{({\it\gamma}-1)(r-1)-4}{({\it\gamma}-1)(r-1)-6}.\end{eqnarray}$$

Assuming, as the most unfavourable case to superdiffusion, $r=4$ , we obtain the values of the superdiffusion exponent plotted in figure 4. We can see that hard electron spectra imply very strong superdiffusion with ${\it\alpha}$ up to 1.88. It is interesting to notice that these values of ${\it\alpha}$ are consistent with those found for electrons in the heliosphere (Zimbardo et al. Reference Zimbardo, Perri, Pommois and Veltri2012). As a further example, let us consider young SNRs, for which ${\it\alpha}_{{\it\nu}}\simeq 0.6$ is typical. This can be explained by DSA with $r\leqslant 4$ . However, many models predict a compression ratio $r\sim 3$ at the thermal subshock (Helder et al. Reference Helder, Vink, Bykov, Ohira, Raymond and Terrier2012; Reynolds et al. Reference Reynolds, Gaensler and Bocchino2012); this would give ${\it\gamma}\sim 2.5$ and ${\it\alpha}_{{\it\nu}}\sim 0.75$ , which is larger than the observed value. Therefore, superdiffusion can allow us to explain the observations in this case too, SSA being able to yield harder spectral indices for a given compression ratio.

The fact that the electron energy spectral index can be inferred from the synchrotron spectral index in the radio band, see (4.2), allows us to infer that superdiffusion can explain those observations. On the other hand, superdiffusion can be a very common process in astrophysics: recently Lazarian & Yan (Reference Lazarian and Yan2014) emphasized the importance of Richardson superdiffusion in MHD turbulence for galactic cosmic rays. Richardson superdiffusion predicts that the spatial separation of magnetic field lines in the direction perpendicular to the mean magnetic field grows as $z^{3}$ , with $z$ the coordinate along the average magnetic field. This superdiffusive regime holds up to field line separations comparable to the injection scale of turbulence, which in the galactic plane is of the order of 100 parsec, and hence much larger than the typical SNR size. This implies that both ions and electrons can undergo perpendicular superdiffusion; even in the case of parallel normal diffusion, $z\propto \sqrt{D_{\Vert }t}$ , the perpendicular separation would grow as $\langle {\rm\Delta}x^{2}\rangle \propto t^{1.5}$ , following the same mechanism of compound diffusion (Webb et al. Reference Webb, Zank, Kaghashvili and le Roux2006). Clearly, the exponent of perpendicular superdiffusion can be larger in the case of parallel superdiffusion. Lazarian & Yan (Reference Lazarian and Yan2014) also found that superdiffusion yields harder energy spectral indices than DSA (although they considered explicitly only the case ${\it\alpha}=1.5$ ), so that this can be a quite general property of cosmic rays accelerated in supernova environments dominated by MHD turbulence.