Stochastic Turbulent Acceleration in a Fractal Environment

Magnetic reconnection, weak turbulence, and shock waves surrounded by passive scattering centers upstream and downstream were for years the prominent acceleration mechanisms in most astrophysical and laboratory plasmas (Melrose 1994, 2009). Recent magnetohydrodynamic (MHD) and kinetic simulations, as well as analytical work, have shown that magnetic reconnection can lead to self-generated turbulence (Matthaeus & Lamkin 1986; Drake et al. 2006; Onofri et al. 2006; Daughton et al. 2011; Oishi et al. 2015; Isliker et al. 2019); driven strong turbulence can also host reconnecting and non-reconnecting current sheets (Biskamp & Welter 1989; Lazarian & Vishniac 1999; Biskamp & Müller 2000; Arzner et al. 2006; Servidio et al. 2011; Isliker et al. 2017), and in shock waves turbulent reconnection will be present mainly downstream (Matsumoto et al. 2015; le Roux et al. 2016; Garrel et al. 2018). Similarly, coherent structures including reconnecting current sheets are now established to be key components of turbulence in magnetized plasmas (Matthaeus & Velli 2011; Cargill et al. 2012; Karimabadi et al. 2013, 2014; Karimabadi & Lazarian 2013; Vlahos & Isliker 2019). In most explosive space, astrophysical, or laboratory plasmas, e.g., flares, unstable astrophysical flows (solar wind and astrophysical jets), or large-scale shocks (bow shock, Heliospheric termination shock, coronal mass ejections, supernova remnants), the heating and acceleration of particles is due to the synergy of large-amplitude magnetic disturbances (stochastic energization) and magnetic reconnection and/or shocks (systematic energization; Comisso & Sironi 2018, 2019; Pisokas et al. 2018).

Acceleration of particles inside fully developed MHD turbulence is a very complex problem and depends on many important factors. (1) The nature of the interaction of particles with the "scattering centers" can be stochastic, systematic, or the synergy of both. The scattering centers inside fully developed turbulence are either large-amplitude magnetic disturbances or coherent structures (current sheets or shocks). (2) The scaling properties of the scattering centers control the energy and space transport and play a crucial role in the acceleration time and the escape time inside the finite acceleration volume.

The processes put forward by Fermi at the beginning of the 1950s to describe particle acceleration inside fully developed turbulence are very broad in nature and include the well-known (i) stochastic (second-order Fermi; Fermi 1949) and/or (ii) the systematic (first-order Fermi; Fermi 1954) process. One can explore these processes by using the concept of a random walk inside a network of scattering centers (Manolakou et al. 1999; Arzner & Vlahos 2004; Vlahos et al. 2004, 2016; Turkmani et al. 2005; Onofri et al. 2006; Isliker et al. 2017; Pisokas et al. 2017, 2018; Garrel et al. 2018; Sioulas et al. 2020).

Fermi (1949) used several simplified assumptions in his analysis of the stochastic interaction of cosmic rays with large-amplitude MHD fluctuations. (1) The interaction of a particle with the large-amplitude magnet fluctuations ("magnetic clouds") is stochastic and the energy gain (δW) is given by the relation (Longair 2011)

$$\begin{eqnarray}&&\displaystyle \frac{\delta W}{W}\sim \displaystyle \frac{2}{{c}^{2}}({V}^{2}-{\boldsymbol{V}}\cdot {\boldsymbol{u}}),\end{eqnarray} \tag{ 1 }$$

where ${\boldsymbol{V}}$ is the characteristic velocity of the magnetic disturbance, ${\boldsymbol{u}}$ is the velocity of the charged particle, and c the speed of light. If ${\boldsymbol{V}}\cdot {\boldsymbol{u}}\lt 0$ the particles gain energy, and if ${\boldsymbol{V}}\cdot {\boldsymbol{u}}\gt 0$ the particles lose energy. (2) The scattering centers are uniformly distributed in space, and the interaction of the particles with the scattering centers is expected to follow Gaussian statistics. The particles execute a random walk with a characteristic mean free path λ_sc between the scattering centers. The acceleration time t_acc was estimated by the relation ${t}_{\mathrm{acc}}\sim [(3c)/(4{V}^{2})]{\lambda }_{\mathrm{sc}}.$ (3) The Fokker–Plank transport equation was used as the basic tool for the study of the statistical evolution of the particles. The transport coefficients were grossly simplified. The acceleration time t_acc is a measure of the energy transport of the particles inside the acceleration volume, and the escape time t_esc a measure of the transport properties in space, t_esc ∼ L²/D, where L is the characteristic length of the acceleration volume and D ∼ λ_scc is the spatial diffusion coefficient if the particles are following a random walk between the scattering centers. Here ${t}_{\mathrm{acc}}/{t}_{\mathrm{esc}}\sim [{c}^{2}/{V}^{2}{L}^{2}]{\lambda }_{\mathrm{sc}}^{2}$; therefore, the power-law index of the high-energy tail in the distribution function is strongly dependent on the mean free path λ_sc. The acceleration time and the escape time are functions of λ_sc and are estimated using the assumptions listed above.

Independently from Fermi's treatment and assumptions, it can be shown that the steady-state solution of the energy continuity equation inside a finite acceleration volume (leaky box approximation) for stochastic Fermi acceleration (Equation (1)) is

$$\begin{eqnarray}&&f(W)\sim {W}^{-(1+{t}_{\mathrm{acc}}/{t}_{\mathrm{esc}})}\end{eqnarray} \tag{ 2 }$$

(Longair 2011).

The results and the simplifications listed above for the stochastic Fermi acceleration have been questioned recently (Pisokas et al. 2017; Sioulas et al. 2020). Both the transport properties in space and energy are not normal and the interactions of the particles with the scatterers follow non-Gaussian statistics.

It is well known that large-scale magnetic disturbances and coherent structures in fully developed MHD turbulence follow monofractal or multifractal scalings, both in space and laboratory plasma (Tu & Marsch 1995; Marsch & Tu 1997; Shivamoggi 1997; Biskamp 2003; Dimitropoulou et al. 2013; Leonardis et al. 2013; Schaffner & Brown 2015; Isliker et al. 2019). Dimitropoulou et al. (2009) examined the relationship between the fractal properties of the photospheric magnetic patterns and those of the coronal magnetic field discontinuities (current sheets) in solar active regions. Isliker et al. (2019) analyzed the current fragmentation of a large-scale current sheet formed during magnetic flux emergence on the Sun and showed that the fragments have a fractal structure, with a fractal dimension D_F = 1.7–1.8.

After all, in fully developed turbulence the coherent structures and the large-amplitude magnetic fluctuations are located on a fractal set with dimension D_F, and the mean free path of the particles with the scattering centers (λ_sc) is not a simple constant (Isliker & Vlahos 2003).

In this Letter, we explore for the first time stochastic Fermi acceleration when the large-amplitude MHD magnetic fluctuations have a fractal structure in space, and the particles are executing a random walk in this environment. In Section 2, we briefly outline the essential characteristics of the random walk in a fractal environment. In Section 3, we present our Monte Carlo simulation model, and in Section 4, we analyze our results. In the final section, we discuss the implication of our results for turbulent stochastic Fermi acceleration.

Isliker & Vlahos (2003) have analyzed the random walk in the environment of a natural fractal, where the fractal is embedded in 3D space and the particles move freely in the empty space not occupied by the fractal until they occasionally collide with parts of the fractal set, where they undergo some kind of scattering. The particles thus move across the fractal, not along it. The fractal is natural in the sense that it is made up of small elementary and finite volumes (and not of points, line segments, etc., as in the case of mathematical fractals), and it also is of finite, though usually large size, such that a clear fractal scaling holds from the fractal's size down to the size of its elementary volumes. The nature of the random walk is illustrated in Figure 1.

Zoom In Zoom Out Reset image size

Isliker & Vlahos (2003) derived the probability density function (pdf) p_F(dr) of the distances dr that a particle travels in between subsequent encounters with the fractal, assuming that initially, a particle resides on a part of the fractal and then moves freely into a random direction until it hits another part of the fractal. For fractals with fractal dimension D_F less than 2 (the case of interest here), this pdf turns out to be of power-law form in good approximation,

$$\begin{eqnarray}&&{p}_{F}(r)=A\,{{dr}}^{{D}_{F}-3},\end{eqnarray} \tag{ 3 }$$

with A a normalization constant, which is a function of the size of the natural fractal and the size of the elementary volumes it is constituted of. With D_F < 2, it follows that the power-law index of p_F(dr) lies in the range −3 < D_F − 3 < − 1, which means that p_F(dr) has the same asymptotic (large dr) functional form as the stable Levy distributions. Particles thus occasionally perform large spatial jumps or "Levy flights," and spatial transport must be expected to be anomalous (Vlahos et al. 2008). A peculiarity of the pdf p_F(dr) is that it is defective, i.e., it is normalized to a value less than one, which implies a finite probability for direct escape in one step, without any secondary encounter with the fractal.

We construct a 3D box of linear size L = 10¹⁰ cm. We initiate the simulation by uniformly placing 10⁶ particles in the interior of the acceleration volume. At time t = 0, the energy distribution of the particles is a Maxwellian with temperature T. We then allow each particle to perform a free flight of length dr_i^(j), before it meets a scatterer (i.e., it undergoes an energization event), where it gains or loses energy stochastically according to Equation (1). The scatterers in our model are assumed to form a fractal set of dimensions D_F = 1.8 (see Sections 1 and 2). From Equation (3), the probability density $P({dr})\sim {{dr}}^{-\gamma }$, with γ = 1.2, yields the length of the spatial step dr ${}_{i}^{(j)}$ each particle performs. We assume that spatial steps range from ${{\lambda }_{\mathrm{sc}}}_{\min }={10}^{2}\,$ cm to ${{\lambda }_{\mathrm{sc}}}_{\max }={10}^{10}\,$ cm. The turbulent volume is a multi-scale environment. The range of the steps used in this study covers the entire range from the kinetic to the MHD scale, the lower limit is of the order of several ion gyroradii, and the upper limit basically equals the size of the acceleration box. Our results are not sensitive to the exact values used for the lower and upper step limits, as long as it holds that ${{\lambda }_{\mathrm{sc}}}_{\min }\ll {{\lambda }_{\mathrm{sc}}}_{\max }$. As a result, there are "long flights," where particles are carried in one step over large distances, in some cases almost through the entire system, before they encounter a scatterer.

To completely specify the coordinates of a particle each time it encounters a scatterer, we also generate a random number for the azimuthal angle ϕ, 0 < ϕ < 2π , and one for $\cos (\theta )$, $-1\lt \cos (\theta )\lt 1$, with θ the polar angle. We then can determine the coordinates of each particle according to

$$\begin{eqnarray*}\begin{array}{rcl}{x}_{i}^{(j)} & = & {x}_{i}^{(j-1)}+{{dr}}_{i}^{(j-1)}\cos \phi \sin \theta \\ {y}_{i}^{(j)} & = & {y}_{i}^{(j-1)}+{{dr}}_{i}^{(j-1)}\sin \phi \sin \theta \\ {z}_{i}^{(j)} & = & {z}_{i}^{(j-1)}+{{dr}}_{i}^{(j-1)}\cos \theta ,\end{array}\end{eqnarray*}$$

where $i=1,2,\ldots ,{10}^{6}$ is the particle index, and $j=1,2,\ldots ,{N}_{i}$ is the number of encounters a particle undergoes, with N_i being the total number of encounters each particle is subjected to before it reaches the final simulation time or escapes from the acceleration volume.

During the free motion, the velocity of a particle remains constant, and since we know the length ${{dr}}_{i}^{(j)}$ and the energy of the particles after an acceleration event, we can keep track of the time elapsed during the free flight as ${{dt}}_{i}^{(j)}={{dr}}_{i}^{(j)}/| {v}_{i}^{(j)}| .$ Therefore, after a total number of j encounters, the time elapsed for each particle is ${\tau }_{i}^{(j)}={\sum }_{j=1}^{j}{{dt}}_{i}^{(j)}.$ We continue to keep track of the particles' energy and transport properties until they reach the final simulation time or cross the boundaries of the box and, therefore, escape from the acceleration volume at time $t={t}_{\mathrm{esc},i}$, which is, of course, different for each particle. In Figure 2, typical orbits in space are presented for a number of selected particles. Obviously, a standard orbit of the particles consists of a combination of long "flights" and efficient "trapping" in localized spatial regions.

Zoom In Zoom Out Reset image size

The conditions we simulate in this article are close to those found in the lower solar corona. We use as strength of the magnetic field B = 100 G, as density of the plasma ${n}_{0}={10}^{9}\,$ cm³, and as ambient temperature T = 100 eV. The Alfvén speed is V_A ∼ 7 × 10⁸ cm s⁻¹, a value close to the thermal speed of the electrons. With these parameters, the energy increments are close to ($\tfrac{\delta W}{W}$) $\sim {\left(\tfrac{{V}_{{\rm{A}}}}{c}\right)}^{2}\sim {10}^{-4}$ (see Equation (1)).

4.1. Spatial Diffusion in the Turbulent Volume

In order to estimate the mean square displacement of the particles, we monitor their positions at prescribed and equi-spaced monitoring times tⁿ (n = 1, ..., N). At time tⁿ a particle's displacement from its initial position is ${\rm{\Delta }}{{\boldsymbol{r}}}_{i}^{n}={{\boldsymbol{r}}}_{i}^{n}-{{\boldsymbol{r}}}_{0i}$ and the mean square displacement for the ensemble of particles is

$$\begin{eqnarray}&&\langle {({\rm{\Delta }}{r}^{n})}^{2}\rangle =\displaystyle \frac{1}{{N}_{p}}\sum _{i=1}^{{N}_{p}}{({\rm{\Delta }}{r}_{i}^{n})}^{2}.\end{eqnarray} \tag{ 4 }$$

We first assume that a particle's encounter with a scatterer solely influences its direction of motion, leaving its energy unchanged. In Figure 3(a), we show the mean square displacement as a function of time. The diffusion for the particles interacting with the passive scatterers is ballistic, the scaling with time has a power-law index close to 2. This result agrees with the results obtained by Isliker & Vlahos (2003, see Figure 10 and Figure 11 therein), where the particles also perform a random walk in an environment where a fractal with dimension D_F < 2 resides.

Zoom In Zoom Out Reset image size

We now turn to the case where the particles gain or lose energy stochastically through their interaction with active scatterers (see Equation (1)). The mean square displacement of the electrons is shown in Figure 3(b), it exhibits a superdiffusive scaling, $\lt {({\rm{\Delta }}r)}^{2}\gt \sim {t}^{2.2}$, with the power-law index decreasing to 1.88 after t ∼ 0.2 s.

In Figure 4(a) we show the distribution of the total number of times the particles encounter a scatterer. The number of encounters strongly varies, ranging from 2 to 2000, with a mean of ∼140 energization events per particle. From Figure 4(b) it is obvious that the particles trapped inside the acceleration volume are those accelerated most efficiently, yet only a fraction of the particles are subjected to a number of energization events that is high enough to be accelerated to superthermal energies.

Zoom In Zoom Out Reset image size

The time spent by the electrons inside the acceleration volume is very important for our study, profoundly affecting the power-law index of the kinetic energy distribution (see Section 1). As Figure 4(c) shows, most of the accelerated electrons escape from the volume quite early, while, for larger escape times, their distribution forms a power-law with an index close to 2.9. The mean value of the escape times yields t_esc ∼ 1.9 s. Comparing this result to Pisokas et al. (2017), where the acceleration process is taking place in an environment where the scatterers are uniformly distributed inside the acceleration volume, we observe a significant decrease in the escape time of the particles.

4.2. Diffusion of Electrons in Energy Space

Equally important for our study are the transport properties of the kinetic energy of the energized particles. In an encounter with a scattering center, a particle (with index i) departs from the scatterer with renewed energy,

$$\begin{eqnarray*}&&{W}_{i}^{j+1}={W}_{i}^{j}+\delta {W}_{i}^{j},\end{eqnarray*}$$

where $\delta {W}_{i}^{j}$ is given by Equation (1), and j counts the number of energization events for the particle. In Figure 5(a), the energization process is presented for several typical particles, revealing its stochastic nature, but also exhibiting a slight predilection for encounters leading to energy gain.

Zoom In Zoom Out Reset image size

Using the set of predefined monitoring times tⁿ outlined in Section 4.1, we keep track of the particles' energies Wⁿ_i at these times. If we denote by W_i⁰ the particles' initial energy, we can define the energy displacement as ${\rm{\Delta }}{W}_{i}^{n}=({W}_{i}^{n}-{W}_{i}^{0})$, and calculate the mean displacement in energy through the relation

$$\begin{eqnarray} &&\langle {\rm{\Delta }}W \rangle ({t}^{n})\equiv \langle {\rm{\Delta }}{W}^{n} \rangle =\displaystyle \frac{1}{{N}_{p}}\sum _{i=1}^{{N}_{P}}{\rm{\Delta }}{W}_{i}^{n},\end{eqnarray} \tag{ 5 }$$

while the mean square displacement in energy is given by

$$\begin{eqnarray} &&\langle {({\rm{\Delta }}W)}^{2} \rangle {({t}^{n})\equiv \langle ({\rm{\Delta }}{W}^{n})}^{2} \rangle =\displaystyle \frac{1}{{N}_{p}}\sum _{i=1}^{{N}_{P}}{({\rm{\Delta }}{W}_{i}^{n})}^{2}.\end{eqnarray} \tag{ 6 }$$

In general, we can assume that the mean energy displacement has a power-law form, $\langle {\rm{\Delta }}W\rangle (t)={F}_{W}{t}^{{a}_{W}}$, and the index a_W can be estimated through a power-law fit. Figure 5(b) shows $\langle {\rm{\Delta }}W\rangle (t)$; there is indeed a power-law scaling with a slope a_W ∼ 0.33 for times up to 0.1 s, and a_W ∼ 0.43 for larger times. The insert figure shows the evolution of the kinetic energy for the electrons remaining inside the acceleration box as a function of time. From the exponential fit we can estimate the acceleration time as ${t}_{\mathrm{acc}}\sim 1/0.38=2.6\,{\rm{s}}$ (Longair 2011). Figure 5(c) presents a_W as a function of W_esc, from which it follows that there is no systematic acceleration for electrons with escape energy smaller than 10⁴ eV. For the high-energy particles, the scaling index gradually increases with energy, reaching a value close to a_W ∼ 1.5.

Similarly, in the case of the mean square displacement in energy we expect a power-law form $\langle {({\rm{\Delta }}W)}^{2}\rangle (t)={D}_{{W}^{2}}{t}^{{a}_{{W}^{2}}}$. In Figure 6 the mean square displacement in energy is presented. For times up to t = 0.1 s, the scaling is slightly subdiffusive, following a power-law with index ${a}_{{w}^{2}}=0.76$. For larger times, the power-law index is ${a}_{{w}^{2}}=1.21$, indicating a superdiffusive behavior. In Figure 6(b) we show a_W² as a function of the energy with which the electrons escape from the acceleration volume. As in the case of convective transport, electrons with energies smaller than 10 keV have on average a scaling index a_w² close to zero. For the superthermal particles, we observe a substantial increase of the scaling index with increasing escape energy, moving from subdiffusive to superdiffusive, even attaining values close to ${a}_{W}^{2}=4$ for the highest-energy particles.

Zoom In Zoom Out Reset image size

In Figure 7(a), we show the histogram of the kinetic energies for the particles that remain in the simulation box, normalized to unity, for the injected distribution and the one at time t = 2.7 s, along with a Maxwellian fit at low energies that yields a temperature T = 10 keV. In the first few milliseconds of the simulation, the low-energy particles are actually already heated, and the high-energy particles are already accelerated and form a power-law tail with index k ∼ 3.8. Figure 7(b) presents the evolution of the power-law index of the tail. After 2.7 s (which is equivalent to the acceleration time t_acc), the initially appearing power-law index k ∼ 3.8 has decreased to an asymptotic value of about k ∼ 2.3, the case shown in Figure 7(a).

Zoom In Zoom Out Reset image size

The power-law index of the tail of the kinetic energy distribution can also be estimated through Fermi's expression $k=1+{t}_{\mathrm{acc}}/{t}_{\mathrm{esc}}\sim 2.37$ (see Equation (2)), which is close to the direct result from the power-law fit in Figure 7(a).

When reducing the size of the acceleration box, the particles do not have time to reach a steady-state distribution before escaping from the box, and the slope of the distribution at high energies becomes steeper and the maximum energy reached smaller. For example, for L = 10⁹ cm the energy distribution remains the same in shape as the one for L = 10¹⁰ cm in Figure 7, yet at an earlier time than in Figure 7(a). Thus, the acceleration time becomes much shorter, the slope of the high-energy tail gets steeper, k ∼ 3.4 at t = 0.2 s (in complete accordance with Figure 7(b)), and the maximum energy reached is 1 MeV. Also, the heated Maxwellian distribution at low energies remains unaffected by a reasonable reduction of the acceleration volume (considering again earlier times than in Figure 7(a)). These results agree very well with the current observations from solar flares and space plasmas (Oka et al. 2018). When increasing the size of the simulation box above 10¹⁰ cm, the energy distribution remains unaffected when comparing at equal times, since the energized particles are able to reach a steady state.

According to Oka et al. (2018), the observed index of the slope of the energetic particles is between 3 and 5 for most solar flares, which, based on our results, suggests that the acceleration box size is about 10⁸–10⁹ cm.

Stochastic turbulent acceleration and transport in space and astrophysical plasma has been analyzed so far with the use of the Fokker–Planck equation and the quasilinear approximation. Both approaches are appropriate for weak turbulence when the wave–particle interaction is a correct representation of the scattering of particles by the normal modes of an unstable plasma. Obviously, in strong and fully developed turbulence these approximations break down since the dominant acceleration mechanisms are large-amplitude magnetic disturbances and coherent structures (current sheets and shocks). Following the initial suggestion by Fermi (1949), we have explored the idea of particle acceleration and heating in the form of a random walk inside a network of scatterers. Fermi assumed that the scatterers (magnetic clouds) are uniformly distributed in space and the mean free path λ_sc is constant. The mean free path plays a key role in the estimates of the acceleration and escape time and controls the power-law index of the high-energy tail. As we outlined in the 1, numerous numerical studies suggest that the spatial scaling of large-amplitude magnetic disturbances and coherent structures inside fully developed turbulence are located on a well-defined fractal topology. We have explored here the role of the fractal scaling in stochastic Fermi acceleration.

The main results in this study are:

• 1.  
  The stochastic interaction of particles with fractal large-amplitude magnetic fluctuations results in the heating and acceleration of particles.
• 2.  
  The high-energy particles are accelerated by a combination of intense trapping within small-scale structures and delayed escape from the acceleration volume, undergoing up to thousands of energization events.
• 3.  
  The combined effects of trapping particles on small scales and of long "flights" dramatically affect the acceleration and escape time of stochastic acceleration. In particular, the acceleration time is strongly reduced when compared to acceleration in non-fractal, uniform environments.
• 4.  
  The spatial and energy transport of the high-energy particles is superdiffusive. The Fokker–Planck equation for the study of the spatial and energy transport of high-energy particles is inappropriate, although is valid for thermal particles.
• 5.  
  The small-scale interactions enhance the acceleration of particles from the thermal pool.
• 6.  
  We simulate in our study explosive phenomena (flares) in the low solar corona, using a simulation box with characteristic length L = 10¹⁰ cm. We have injected a very large number of electrons with a Maxwellian energy distribution with a temperature of 100 eV. In about 2 s, the energy distribution reaches an asymptotic shape, with a superhot plasma with temperature 10 keV, and a power-law tail above 100 keV with a power-law index −2.3, and reaching 10 MeV.
• 7.  
  When reducing the size of the box, e.g., to 10⁸–10⁹ cm, the particles do not have time to reach a steady-state distribution before escaping from the system, and the power-law slope of the high-energy tail becomes steeper, in agreement with the current observations from solar flares and space plasmas (Oka et al. 2018). Increasing the size of the acceleration box to L > 10¹⁰ cm does not affect the energy distribution, since the particles in any case can reach a steady-state distribution.

We confined our study to the stochastic Fermi acceleration of particles in a fractal turbulent environment, which turned out to be a very efficient and important mechanism for many turbulent astrophysical sources, beyond the case of solar flares studied here. Our next step is to incorporate coherent structures (reconnecting current sheets), as they are present in fully developed plasma turbulence.

We thank Theophilos Pisokas for his help in the initial phase of this project.