Relativistic Particle Transport and Acceleration in Structured Plasma Turbulence

We discuss the phenomenon of energization of relativistic charged particles in three-dimensional (3D) incompressible MHD turbulence and the diffusive properties of the motion of the same particles. We show that the random electric field induced by turbulent plasma motion leads test particles moving in a simulated box to be accelerated in a stochastic way, a second-order Fermi process. A small fraction of these particles happen to be trapped in large-scale structures, most likely formed due to the interaction of islands in the turbulence. Such particles get accelerated exponentially, provided their pitch angle satisfies some conditions. We discuss at length the characterization of the accelerating structure and the physical processes responsible for rapid acceleration. We also comment on the applicability of the results to realistic astrophysical turbulence.


INTRODUCTION
The fact that the cosmic-ray (CR) spectrum extends up to extremely high energies, as well as the difficulties encountered by standard acceleration mechanisms to energize particles up to such energies, has led to a wide search for ways to boost the maximum energy and/or alternative acceleration mechanisms. While we can safely say that CR acceleration at the shocks associated with supernova (SN) explosions is now confirmed experimentally, it is also true that convincing observational proof that SN remnants (SNRs) can accelerate protons to energies in excess of ∼ 100 TeV is, thus far, missing. Even from the theoretical point of view, SNRs associated with Type Ia and ordinary core-collapse SNe may be connected to CR acceleration up to ∼ 100 TeV, while higher energies require much more extreme conditions, perhaps to be found in rare and very energetic SN explosions (Cristofari et al. 2020).
While standard acceleration processes, such as the secondorder Fermi process (Fermi 1949) and diffusive shock acceleration (Axford et al. 1977;Krymskii 1977;Blandford and Ostriker 1978), have received a lot of attention throughout Corresponding author: Oreste Pezzi oreste.pezzi@istp.cnr.it the years, more recently, the acceleration of charged particles in realistic MHD turbulence has been attracting an increasing level of attention (e.g., Lazarian et al. 2020, and references therein), with special emphasis on magnetic reconnection. Early on, attention focused on contributions to energization by MHD activity in special regions within a dynamic plasma (Giovanelli 1947;Dungey 1953). Examination of diverse scenarios for energization has over time become increasingly varied and complex (Ambrosiano et al. 1988;Dmitruk et al. 2004;Drake et al. 2006Drake et al. , 2009Kowal et al. 2012), suggesting that a more general perspective may be led by a simpler view of the physical principles at work.
In 2D configurations, there exist simple conservation laws that define and constrain particle orbits (Sonnerup 1971;Matthaeus et al. 1984) so that when turbulent fluctuations are present, charged test particles can be confined and accelerated (Ambrosiano et al. 1988;Drake et al. 2006) in secondary magnetic flux structures or "islands". Initially applied to smaller-gyroradius particles such as electrons in the context of magnetic reconnection geometries, this idea was further developed by Oka et al. (2010), who noted that coalescence in a multiple-island context can efficiently accelerate electrons due to "anti-reconnection" electric fields associated with such mergers. Furthermore, when applied to "pickup" of protons in reconnection jets (Drake et al. 2009) that subsequently feed into a Fermi process as multiple magnetic islands contract and merge, a more complex proton energization process may be developed involving multiple islands (Drake et al. 2010). More recently, Trotta et al. (2020b) showed that transrelativistic electrons can be significantly accelerated while trapped in turbulent structures and experiencing curvature drift. Such a behavior, obtained in a 2D configuration where the background turbulent plasma is modeled with a hybrid particle-in-cell method while energetic electrons are treated as test particles, has also been verified at different intensities of the turbulent fluctuations. The basic physical elements of these models were formalized in transport theories (Zank 2014;le Roux et al. , 2018) that facilitate applications (see also the recent reviews by Khabarova et al. (2021); Pezzi et al. (2021b)).
Complementary to these developments, there has been a parallel line of studies that begin with the premise that charged particle energization might be treated as being due to plasma dynamics that from the onset is complex and, in fact, turbulent. Early efforts demonstrated the feasibility of efficient turbulent acceleration (Dmitruk et al. 2003) in three dimensions, although numerical limitations such as small system size and lack of resonant power at numerical grid scales (however see (Lehe et al. 2009)) cast doubt on conclusions concerning scaling laws.
Numerical experiments have improved ever since (Dmitruk et al. 2004;Kowal et al. 2012;Dalena et al. 2014), including turbulence effects associated with current sheets and reconnection, contracting islands, and proliferation of plasmoids and/or secondary islands (Kowal et al. 2011(Kowal et al. , 2012. Similar effects have been characterized in relativistic plasmas (Hoshino 2012;Zhdankin et al. 2017Zhdankin et al. , 2019Comisso and Sironi 2018;. Most important, it has been understood that the initial efforts to characterize acceleration in or around reconnection regions focused on the physical processes responsible for extracting particles from the thermal background and injecting them into some energization process at work on larger scales. Yet, the particle velocity, the thermal velocity, and the Alfvén speed were of the same order of magnitude, making it very difficult to identify the nature of the acceleration mechanism (first or second order) and the scaling laws that could be used for higher-energy particles.
Injecting test particles in a 2D or 3D simulation box of MHD turbulence has been an independent method to investigate both the transport of such particles and their energization, although in some of the previous articles (e.g. (Kowal et al. 2011(Kowal et al. , 2012) the Larmor radius of the particles at the beginning of the simulation was chosen to be much smaller than the grid spatial spacing, which makes transport unrealistic, due to the lack of resonance with the turbulence and thereby unreliable diffusive transport.
When self-consistent broadband turbulence is involved, particle energization becomes complex due to interactions with the internal structure of large-and small-scale flux tubes, as well as the current sheets, vortices, and reconnection sites that typify flux-tube boundaries and their mutual interactions. Common to a number of these treatments of energization including turbulence is the role of direct acceleration for particles of smaller Larmor radii, transitioning to the involvement of perpendicular acceleration at larger energy (Dmitruk et al. 2004;Dalena et al. 2014;Sironi 2018, 2021;Trotta et al. 2020b).
Although some works indicated that high energies can be also attained by second-order processes (e.g., Arzner et al. (2006); Sioulas et al. (2020)), the role of temporary trapping has been highlighted in various contexts as it can dramatically influence both transport and acceleration (Kowal et al. 2011(Kowal et al. , 2012Tooprakai et al. 2016;Dalena et al. 2014). Indeed, the ubiquity of turbulent coherent structures and islands/plasmoids/flux ropes produced by magnetic reconnection makes the potential of the trapping mechanism significant for different systems. In the solar and stellar coronae, magnetic loops may provide sufficient conditions to entrap particles (Vlahos and Isliker 2018). In the interplanetary medium, there are observations supporting the idea that particles are locally accelerated when trapped in merging or coalescing islands (Khabarova and Zank 2017;Malandraki et al. 2019). Particle trapping may also provide a source for particle reacceleration downstream of shocks Nakanotani et al. 2021). Other systems, e.g. the intracluster medium, could also encompass such a mechanism although observations are usually explained in terms of second-order processes (Vazza et al. 2016;Brunetti and Vazza 2020). The importance of trapping for accelerating particles has been also recently highlighted by Lemoine (2021) in the context of strong and intermittent turbulence by relating the energization process and the gradients of the bulk velocity and magnetic fields.
Intuitively, trapping can enhance acceleration when appropriate electric fields are encountered, but it can also inhibit stochastic acceleration when the required transport is thwarted. The full range of possibilities for such effects remains to be exhaustively explored. The study of Kowal et al. (2011, see also Kowal et al. (2012) is instructive on several salient features. In these numerical experiments, the authors observe two small magnetic islands merging with a central elongated island. A selected particle gains considerable energy after entering the system by following field lines through one small island, eventually becoming trapped within the large island, circuiting it numerous times, and gaining energy exponentially and mainly when passing near the region of merger with the small island.
These conclusions were illustrated in detail only for 2D turbulence in a configuration that was optimized to create reconnecting islands. However, the test particles injected in a snapshot of the simulation were initially subgrid, which means that their transport could not be described in a realistic way. Eventually, the energy of the particles, after an exponential increase, reaches the regime in which resonances could in principle be relevant for particle transport, but this phenomenon was not discussed. Despite these shortcomings, this approach demonstrates clearly the complex interplay between the transport effects that entrain the particle within the island and near the acceleration region, along with the special circumstances that support the electric field responsible for the energization itself.
In the present paper, we advance such a scenario for strong turbulence and particle energization in the highly relativistic particle regime, having in mind the implications that the process may have for the acceleration of very-high-energy CRs. In particular, we continue the examination of these complex interactions of charged particles with turbulence by performing relativistic test-particle simulations with a turbulent electromagnetic field produced by means of 3D MHD simulations, not specifically devised to produce reconnection regions. Our focus is on clarifying of the nature of the trapping or entrainment that leads particles to rapid acceleration to the highest energies.
We find that the bulk of the test particles injected in the simulation box went through a secular-second order acceleration due to the random plasma motions, which in turn induced random electric fields. The interaction of the particles with these plasma motions leads to stochastic energy change, which we characterize in terms of plasma properties and manage to associate with a diffusive motion in momentum space. The transport of the same particles in physical space is also found to be well described through a diffusive motion. The latter is in a range of scales where the anisotropic cascade of the turbulence with respect to the local magnetic field does not seem to have a visible effect as yet.
In addition to this acceleration process that is clearly at the second order in the quantity v A /c 1, we also identify a small fraction of particles that manage to get trapped in selected regions associated with the interaction between flux tubes. These particles all have pitch-angle cosine very close to zero, a necessary condition for trapping, and go through an exponential phase of energy increase. We provide a characterization of these regions in terms of physical observables that could be measured in the simulation. We also build a simple model of the region where this phenomenon occurs and manage to reproduce the main properties of the acceleration process and the time scales involved. The acceleration process is similar to a first-order mechanism in which the particle trajectory in the plane perpendicular to the local magnetic field encounters a gradient of plasma velocity (i.e. of the induced electric field).
The paper is structured as follows. In Sect. 2 we present the MHD simulations adopted in the present work, while in Sect. 3 we describe the test-particle code and the first results obtained in terms of physical space transport. Then, in Sect. 4 we focus on the main numerical outcomes of the work concerning particle energization. Sect. 5 discusses a simple model of the acceleration region and derives the main properties of the acceleration mechanism and the main time scales involved. Moreover, we discuss the implication of our findings for astrophysical systems. Finally, in Sect. 6 we conclude by summarizing our results and illustrating future developments.

MHD SIMULATION BACKGROUND
In order to study the transport and acceleration of charged test particles, we follow particle evolution in electromagnetic fields obtained through incompressible three-dimensional MHD simulations. These simulations solve the following set of equations: where u(r, t) is the magnetofluid speed composed only of its fluctuating part, and B(r, t) is the magnetic field that is decomposed into a uniform mean B 0 and a zero-mean fluctuation b, B(r, t) = B 0 + b(r, t) = B 0 e z + b(r, t). Furthermore, P is the thermal pressure, and ρ is the magnetofluid density. The current density is j = ∇ × B, while ν and η are the viscosity and resistivity, respectively. The flow is incompressible ∇ · u = 0, and the density is uniform ρ = const. Lengths, time, and velocities in Equations (1-3) are respectively normalized to a typical length L A , time t A and to the Alfvén speed v A = L A /t A =B/ 4πm pn , whereB andn are reference values for the magnetic field and for the background number density. We here adopt L A = 81.5 pc, corresponding to L box = 512 pc, whileB = 1 µG and n = 1 cm −3 . Unless specified, hereafter we assume normalized variables.
Equations (1-3) are solved in a 3D Cartesian periodic box of size L box = 2π, with spatial resolution N x = N y = N z = 1024 adopting a pseudo-spectral method in a Fourier basis. The time advancement is performed with a second-order Runge-Kutta scheme and the 2/3 rule for spatial dealiasing is chosen (Patterson and Orszag 1971). Small values of resistivity and viscosity η = ν = 2 × 10 −4 are introduced to define the well-resolved spectral domain. The dissipative Figure 1. Rendering of the current density j 2 (r) shows a plethora of intermittent coherent structures. Such structures form a template for the possibility of rare acceleration events.
wavenumber k diss (the reciprocal of the Kolmogorov length scale) for the considered runs is always smaller by a factor 2 than the maximum resolved wavenumber k max (for further details, see Bandyopadhyay et al. 2018).
Large-scale uncorrelated fluctuations of u and b are introduced at t = 0 and turbulence develops, producing smallscale fluctuations. We focused here on the case with u rms = b rms = 1 and B 0 = 0. The role of a finite background magnetic field and compressibility will be discussed in a separate forthcoming work. We then selected the time instant at which the turbulent activity is strongest (i.e. highest dissipation). The complex and highly structured pattern of the turbulence is displayed in Figure 1, showing the contour plot of j in the 3D domain. Vortices and magnetic islands, as well as intense current sheets where magnetic reconnection may be at work, naturally emerge as elementary structures of the turbulent flow. The omnidirectional spectrum of magnetic energy ( Figure 2) indicates that an inertial range, whose length is about a decade in wavenumber space, develops before dissipative effects steepen the spectrum at higher wavenumber k. In the inertial range, the slope is rather compatible with either the Kolmogorov or Kraichnan predictions; these are respectively displayed in green and orange dashed lines in Figure 2 (see also the inset in the same figure). By numerical evaluating the correlation length l c of the magnetic field, we find l c = 0.218 (l c = 17.7 pc in physical units), corresponding to protons with energy E ∼ 16 PeV in the typical fieldB.

METHODS: TEST-PARTICLE PROPAGATION DETAILS
We numerically integrate the motion equations of N p = 10 5 relativistic test particles of positive charge e and mass m p moving in the turbulent electromagnetic field obtained Omnidirectional spectrum of the magnetic energy. The green dashed (orange dotted-dashed) line shows the Kolmogorov (Kraichnan) prediction. The small inset displays the magnetic energy spectrum compensated by the Kolmogorov (green dashed) and the Kraichnan (orange dotted-dashed) slope. The gray dashed vertical lines indicate the wavenumber associated with the initial particle gyroradius.
by means of the incompressible MHD simulations described above. The normalized particle equations of motion are where x = (x, y, z), v, and p = γv are the particle position, velocity, and momentum, while E and B are the electric and magnetic fields. Equations (4-5) are scaled analogously to MHD simulations. In normalized units, the Lorentz factor reads γ = 1/ 1 − (β A v) 2 = 1 + (β A p) 2 , where β A = v A /c. The electric field in Equation (5) is derived through Ohm's law: E = −u × B + ηj. The parameter α = t A Ω 0 , where Ω 0 = eB/m p c is the proton cyclotron frequency, can be easily rewritten as α = L A /d p , with d p the proton skin depth of the background plasma. α is thus connected to the extension of the inertial range of the turbulence with respect to kinetic, dissipative scales (Dmitruk et al. 2004;González et al. 2016). In a β p ∼ 1 plasma (with β p the thermal to magnetic pressure ratio), the parameter α corresponds to the inverse of the normalized gyroradius of nonrelativistic particles moving with speed ∼ v A . Previous works considering the injection of thermal particles into the acceleration region were hence forced to reduce α to much smaller and computationally feasible values. Such a requirement provides particles with a gyroradius at least larger than the grid size, so that resonant scattering might be properly taken into account. On the other hand, relativistic particles moving at the speed of light have a much larger gyroradius because γ 1, thus removing the constraint on the value of α. For the parameters described above, α ∼ 10 12 and β A ∼ 10 −5 . To save computational resources, we only artificially increase β A = 5 × 10 −2 .
Because we are interested in the energization of relativistic particles moving in a nonrelativistic environment, we assume stationary electromagnetic fields, i.e. ∂B/∂t = ∂E/∂t = 0 (magnetostatic approximation), and we consider a static snapshot of these fields when turbulence is fully developed.
Eqs. (4-5) are integrated by adopting the relativistic Boris method (Ripperda et al. 2018;Dundovic et al. 2020). The electric and magnetic fields are interpolated at the particle position through a trilinear interpolation method (Birdsall and Langdon 2004). We verified that the results presented here are not affected by adopting a more accurate yet significantly slower 3D cubic spline method (not shown here). Particles are injected homogeneously throughout the computational box at a given energy and with isotropic velocity direction on the unit 3D sphere. The time step is set to 1/50 of the initial gyroperiod.
Most of the results here adopt the initial gyroradius to be r g,0 0.1l c = 0.02, corresponding to E 0 1.6 PeV. The resonant wavenumber k rg = 1/r g = 50, reported in Figure 2 with a vertical dashed gray line, resides in the inertial range of turbulence, and it is also quite far from the dissipative scales where the resistive electric field is expected to become important. This ensures that the acceleration process studied here is mainly driven by the inductive electric field, this being the most relevant term for analyzing the energization of relativistic particles whose gyroradius is much larger than the typical length where dissipative and resistive effects are expected to steepen the magnetic spectrum. To doublecheck, we also verified that our results are not affected by the resistive field, in that if we exclude the resistive component from the computation of the electric field, the energization process is basically unchanged. This shows that for the highenergy particles we are interested in, namely when the Larmor radius exceeds the thickness of the reconnection regions, the energization is not due to the resistive fields but rather to the induced electric fields due to the plasma motion.
Although here we are most interested in how particles react to electric fields in the simulation box, it is first worth studying how particles move in the magnetic field, especially to confirm that we find diffusive motion and to identify possible differences with respect to cases where turbulence is synthetic rather being the result of an MHD simulation.
In order to study particle transport in physical space, we performed a subset of test-particle simulations by excluding the electric field. A diffusive regime after a ballistic transient is always recovered. When reaching the diffusive plateau, the isotropic diffusion coefficient is computed as 10 −2 10 −1 10 0 10 1 rg/lc  Figure 3 shows the isotropic mean free path λ iso = 3D iso /c as a function of the particle gyroradius r g , normalized in the usual way to the correlation scale l c . The typical behavior of the path length as a function of energy is the same as that found in synthetic turbulence, with a lowenergy trend that reflects the shape expected from a given isotropic power spectrum (Subedi et al. 2017;Dundovic et al. 2020). In particular, the dotted-dashed red line in Figure 3 implements Equation (53) of Dundovic et al. (2020) for the Kolmogorov case, where l iso is the bend-over scale of the synthetic model in Dundovic et al. (2020), being l iso ∼ 2l c . At variance with synthetic models of the turbulent field, the slope of the isotropic power spectrum is not very well defined here because of the limited dynamical range (see Figure 2).
At high energies, when the gyroradius satisfies the condition r g > l c , the diffusion coefficient becomes weakly dependent upon the power spectrum, D iso ∼ r 2 g (blue dotteddashed line). For gyroradii r g /l c < 0.01, several numerical effects, especially dissipation, limit the validity of the approach by reducing the power available in the form of modes that the particle gyration can resonate with. As a result, the numerically computed path length departs from the dotteddashed red line at such low energies. The vertical dashed gray line identifies the energy of the particles used below for our investigation of energization.
This correspondence to the isotropic spatial diffusion theory is itself a result of some significance: on one hand, it validates the approach of Subedi et al. (2017) and Dundovic et al. (2020) with a "realistic" turbulent and intermittent magnetic field obtained from the numerical evolution of MHD equations. On the other hand, this may be considered somewhat surprising because the anisotropic cascade development with respect to the local magnetic field, expected in MHD (Sridhar and Goldreich 1994), seems to have little effect on the diffusion properties. It is likely that the effects of the anisotropic cascade are not fully developed as yet, due to the limited dynamical range imposed by the numerical constraints. These results seem to be in good agreement with those of Cohet and Marcowith (2016).

RESULTS ON PARTICLE ENERGIZATION
The turbulent motion of the plasma in the simulation box leads to the unavoidable creation of inductive electric fields, which are expected to have random orientations. As such, their presence is expected to lead to changes in the energy of the particles, due to the presence of such inductive electric fields in the Lorentz force. This effect is expected to cause the energy of a particle to increase or decrease depending on the relative orientation of the particle momentum and the local electric field, namely a typical second-order phenomenon. On the other hand, the complex structure of MHD turbulence is known to trigger additional phenomena that may lead to more rapid energization of the particles (see for instance Kowal et al. (2011)). Here we investigate all these phenomena in great detail, stressing that the simulation was not carried out to maximize the formation of reconnection regions or other peculiar structures. The phenomena we see are, in this sense, very generic.
We find that in, addition to an overall second-order stochastic acceleration mechanism, a first-order process is at work as well, due to the temporary trapping of particles in coherent structures. Figure 4 displays the running energy diffusion coefficient as a function of energy, assuming that in fact the motion of the particles can be described in terms of a random walk in momentum space:

Stochastic energization
It is evident that after a transient, the energy diffusion coefficient saturates at a roughly constant value. This implies the presence of diffusion in energy space, thus revealing the typical nature of a second-order process (Ostrowski and Siemieniec-Ozieblło 1997). The energy diffusion coefficient D EE 0.01v A E 2 0 /l c implies a characteristic time for the energy diffusion process τ diff,E = E 2 /D EE ∼ 10 2 l c /v A that is consistent with the large time scale for growth of the average energy that occurs at late times in our simulations (not shown here). The slight increase recovered in D EE for very large ∆t may be due to the fact that the average gyroradius starts to increase on this timescale due to other processes (see below), thus making D EE move away from the plateau.
Another signature of an active stochastic energization mechanism comes by looking at the probability density 10 −2 10 −1 10 0 10 1 ∆tv A /lc functions (PDFs) of the relative energy gain ∆E(∆t)/E, displayed in Figure 5 for ∆t = 0.05l c /v A (blue) and ∆t = 0.5l c /v A (black). The PDFs have been computed averaging on the initial time instants t up to t = t max 22l c /v A . Particles are likely to undergo both increases and decreases in energy. Although the distribution is peaked at small values, larger changes of energy up to β A (dashed gray lines in Figure 5) are allowed. The distribution functions are skewed toward the positive value of energy changes because the standardized skewness iss = 0.20 for ∆t = 0.05l c /v A ands = 0.12 for ∆t = 0.5l c /v A , wheres = s/σ 3 . Here, σ and s are, respectively, the standard deviation and the skewness (third-order moment) of the distribution function. The distribution function of the relative energy gains is manifestly non-Maxwellian for small ∆t and tends to recover the Maxwellian shape for larger ∆t. Indeed, the kurtosis κdefined as the fourth-order moment of the PDF normalized by σ 4 -is κ = 4.84 and κ = 3.69 for ∆t = 0.05l c /v A and ∆t = 0.5l c /v A , respectively.
The presence of positively skewed PDFs indicates that energy increase is more favorable than energy decrease. This provides the secular direction of the process, leading to a net energy gain.
Clearly, all acceleration processes are at work simultaneously, and it is not trivial to discriminate among them by looking at a collection of particles: What we can say is that basically all test particles launched in the simulation box suffer from the second-order process illustrated above. As we discuss below, a small fraction of particles happen to be trapped in selected structures and get energized through a first-order process. It is not clear to what extent these few particles can affect the shape of the high-energy-gain tail of the PDF shown in Figure 5. As a consequence of the fact that only a few particles experience trapping, it is in general rather difficult to evaluate statistical properties of the population of trapped particles, such as the transport coefficients. In future work, we plan to explore the potential of novel meth- ods, for instance, those commonly adopted in biophysics and based on single-particle trajectories (Golding and Cox 2004;Saxton 2012;Trotta et al. 2020a).

Particle trapping in coherent structures: First-order acceleration
In Figure 6 we show the spectra of particles in the simulation box, after a time t indicated in the figure. A few comments are in order: (1) Particles are injected at energy E 0 which, for the natural units adopted here, corresponds to E 0 = 1.6 PeV.
(2) As time evolves, the second-order process leads to a broadening of the distribution function, namely there are both particles losing energy and particles gaining energy. On average, however, the particle energy increases as one can see by noticing that the peak of the distribution moves toward higher energies. (3) Contemporaneously an approximate power law is created at high energies that eventually extends to particles with energies such that r g l c . This typically happens at times t 10l c /v A .
For the sake of comparison, we also report in Figure 6 a line indicating the spectrum ∝ E −5/2 , which was predicted and observed by different groups, although for very different systems and with different qualitative premises.
The pioneering numerical simulations of Ambrosiano et al. (1988) found evidence for accelerated particles with a powerlaw tail with a slope compatible with −5/2 using twodimensional simulations. As we discuss below, the dimensionality of the problem is very important in assessing the efficiency of trapping processes that are responsible for particle energization. Moreover, Ambrosiano et al. (1988) focused on the extraction of particles from the thermal bath, for which the particle velocity remains close to that of the background thermal particles. As we discuss below, the acceleration, the trapping, and escape from the acceleration region work in somewhat different ways for relativistic particles. Moreover, as it is well known, the slope in energy should not be the same in any case for relativistic and nonrelativistic particles: This is so even for particles accelerated at a strong shock, for which the spectrum of accelerated particles in momentum is f (p) ∝ p −4 , but when expressed in terms of energy, it is N (E) ∝ E −2 for relativistic particles and ∝ E −3/2 for nonrelativistic particles. The slope −5/2 was also predicted by de Gouveia dal Pino and Lazarian (2005) and del Valle et al. (2016), where a shock-like toy model was introduced to describe a reconnection region: The particles would be advected into the reconnection region with the inflowing plasma and would be expelled (the analog of escape to downstream in the case of a shock) at the speed of the reconnection exhaust. This determines a sort of universal spectrum for the accelerated particles. Such universality was later criticized by Drury (2012). In fact, for particles that are being energized, it is unlikely that the velocity of the exhaust plays any role in the shaping of the spectrum of accelerated particles because the particles' Larmor radius becomes quickly larger than the thickness of the current sheet. As we mentioned above, the spectrum E −5/2 is shown in Figure 6 only as a reference, while it appears to be asymptotically reached in our simulations only for exceedingly long times compared with the dynamical time of the MHD turbulence.
In the perspective of understanding the nature of the acceleration processes at work, in the top panel of Figure 7, we show the temporal evolution of the gyroradius averaged on the full particle ensemble (red dashed line), while the red shadowed area corresponds to the standard deviation of the averaged gyroradius. We clearly see that there is a secular increase in the particles' energy, which we attribute to the random interaction with the inductive electric fields in the simulation box. The black curve shows the temporal evolu-tion of the particle gyroradius that, after a time T ∼ 22l c /v A , turns out to be the most energetic particle in the simulation. It is worth noticing that for early times, namely on the left side of the first vertical dashed line (t ≈ 8l c /v A ), the fluctuations in the particle gyroradius (or its energy) are compatible with the fluctuations expected based on the bulk of the particles in the simulation (shaded area). Moreover, the particle energy is clearly carrying out a random walk, in that it increases and decreases, a typical feature of a second-order process.
At a time t ≈ 8l c /v A an exponential increase of the particle energy starts (the plot is in lin-log scale) and lasts for about ∼ 10l c /v A . The end of this period is marked by the rightmost vertical dashed line (t ≈ 17.5l c /v A ). This period of rapid energization suggests that a small number of particles experience some new phenomenon. This number must be small because the energy gained by such particles is visibly larger than the typical deviation from the mean (shadowed area).  Figure 7. Typical behavior of a trapped particle showing exponential growth of energy that occurs within the two vertical green dashed lines. The particle gyroradius for the trapped particles increases exponentially over a time-scale τ 10lc/vA. This growth is much faster than the growth of the averaged gyroradius, where the average is performed on the full ensemble of test particles (red dashed line). The dashed red area represents the standard deviation of the averaged gyroradius. The bottom panel shows the particle trajectory illustrating that the particle is trapped.
We notice that, superposed on the main regular energy growth, there remain visible smaller-scale oscillations of the particle energy (see inset in Figure 7). In particular, we visually identify a smaller-scale oscillation whose period is T 5 × 10 −2 l c /v A , which corresponds to the particle gyromotion τ g = 2π/Ω g = 2πr g /c. A larger-scale modulation with period T 5 × 10 −1 l c /v A is also observed, and this may be correlated with fluctuations of the magnetic field intensity on this timescale. This is suggestive of the simultaneous presence of additional processes, such as mixing of second-/first-order processes and the role of mirroring or drifts. The multiscale complexity of the overall energization dynamics is evidenced by the appearance of at least four timescales in Figure 7 -the exponential time scale, the gyromotion, the modulation seen in the inset, and the secondorder energy gain seen to the left and right sides of the exponential phase.
The peculiar behavior of the particles during the exponential phase is best illustrated in the bottom panel of Figure 7, where we show the particle's trajectory. One can see that, during the stage of exponential energy growth, the particle is trapped in a small region of the computational domain of size ∼ 0.5l c . In fact, the spatial excursion per unit l c /v A is about 10 times smaller between the dashed lines than outside that region. The particle escapes from the trapping region when its gyroradius becomes comparable with the island size l isl ∼ l c .
The phenomenology of this trapping can be also appreciated by looking at the particle trajectory in the 3D domain. Figure 8 shows the particle trajectory as dots colored with the particle energy, where the color scale goes from blue to red as particle energy increases. Magnetic field lines near the trapping region are also displayed, colored with the amplitude of the field itself (again going from blue to red as the magnetic field amplitude increases). When the particle is not trapped, it carries out an erratic motion in the whole computational domain, akin to an unconstrained random walk. The trapping is associated with a spherical-like motion constrained within a flux-tube-like structure. The energization occurs when the flux tube is perturbed by another large-scale structure, more easily appreciated in the right panel of Figure 8. This confirms the scenario that an intense acceleration can occur when magnetic islands and, more in general, largescale plasma structures, are interacting (collapsing, merging, etc) with other similar structures (Drake et al. 2006;Kowal et al. 2011), leading to a locally strong magnetic field gradient.
It is important to point out here that there is no evident association of the structure responsible for the exponential growth of the particle energy with the process of magnetic reconnection. Magnetic reconnection is a sufficient condition for generating large-scale islands where particles can be Figure 8. Particle trajectory in the 3D domain, with the points colored with the particle energy, where the color scale goes from blue to red as the particle energy increases. Magnetic field lines, colored with the magnitude of the magnetic field itself (again from blue to red as the magnetic field magnitude increases), indicate that the particle is trapped in a flux tube and that it is accelerated when the flux tube is feeling the gradients associated with the interaction with another large-scale structure. The right panel shows an inset of the left plot zoomed in the trapping region and limited in time to a few particle gyrations. The green line in each panel corresponds to the correlation length lc.
trapped. Indeed, it can be expected that when the magnetic field reconnects in a turbulent environment, the magnetic islands produced by reconnection interact, thus allowing an intense and fast energization process. However, it is apparently not necessary that reconnection be present during the energization process itself. Other configurations without the explicit invocation of magnetic reconnection, such as the interaction of two large-scale turbulent structures (e.g. flux ropes, as recently reported in recent Parker Solar Probe observations by Pecora et al. (2021)), may provide a similar behavior, provided that the magnetic geometry of the interaction region favors particle trapping. We remark that the direct acceleration due to the electric field at the reconnection site is negligible for the relativistic particles considered in the present work, given that such particles have a gyroradius much larger than the typical width of current sheets.

Characterization of trapping and concomitant energy gain
In order to characterize the coherent structure that entraps and gives a significant boost to the particle energy, we calculate the current density j = ∇ × B and the normalized magnetic helicity h m = a · B/(|a||B|), with B = ∇ × a, interpolated at the particle position. The current density is a direct proxy of the small-scale gradients of the magnetic field, and an intense current density is expected to highlight small-scale structures and current sheets where magnetic reconnection and, in general, dissipative processes may occur (see Pezzi et al. 2021a, and references therein). On the other hand, magnetic helicity measures the topology of the magnetic field: In particular, a nonnull magnetic helicity indicates twisted, helical magnetic structures.
The values of these variables at the particle position are displayed in Figure 9. The exponential phase is limited by the green dashed vertical lines. The structure responsible for the exponential growth of the particle energy is a relatively quiet region in which the current density is relatively smooth. In comparison, the current outside the structure easily reaches intense values j 4j rms , but such intense peaks are not evident within the structure. The magnetic fluctuations are also less intense within the structure, as the rms value of magnetic fluctuations is reduced there by a factor of 3 − 4 with respect to the global value.
The structure is furthermore characterized by a finite magnetic helicity, suggesting a flux-tube and/or plasmoid-like shape where magnetic field lines wrap helically on themselves. A finite magnetic helicity also suggests that the structure tends to be force free as a||b → j||b, i.e. it may be a large-scale quasi-equilibrium structure typical of intermittent plasma turbulence, where nonlinearities are depleted (Matthaeus et al. 2015). The properties of the particle trapped in the accelerating coherent structure are also remarkable. The top panel of Figure 10 illustrates the pitch-angle cosine of the particle, here defined as because the regular field is absent. Concurrently with the period when the particle is trapped, its pitch-angle cosine displays reduced oscillations around the mean zero value. In fact, µ loc oscillations, estimated as (∆µ loc ) rms , are weaker by a factor of 3 − 4 inside the structure with respect to outside.
The trapped particle has a peculiar motion configuration, with a pitch angle almost perpendicular to the local field. This indicates that the particle is trapped within an elongated 2D-like flux tube and, in particular, it moves in the plane perpendicular to the tube axis. As a consequence, within this period, the particle mainly experiences mainly a perpendicular energization. The evolution of the particle magnetic moment (not shown here) also supports this view because it shows a secular growth on the same timescale as the exponential energy growth of the particle, while the magnetic field is roughly constant within the same time window. Similar observations have been pointed out in the different contexts of the so-called "second stage" of acceleration of nonrelativistic particles (Dalena et al. 2014) and, more recently, of electron acceleration in nonrelativistic plasma turbulence (Trotta et al. 2020b) and particle acceleration in relativistic plasma turbulence including radiative losses (Comisso and Sironi 2021). To get insights into the nature of µ loc oscillations in the time period corresponding to trapping, the bottom panel of Figure 10 displays the Fourier time spectrum of µ loc , performed in such a window. The dashed yellow area corresponds to the frequencies associated with the particle gyroradius, which changes within the period due to the particle energization. High-frequency µ loc oscillations, associated with the particle gyromotion, are combined with lower-frequency fluctuations, possibly related to smaller-scale turbulent fluctuations or other effects, such as mirroring and drifts.
To demonstrate that the most intense energization is statistically associated with a small pitch-angle cosine, Figure  11 displays the PDFs of the pitch-angle cosine conditioned to the particle energy. In particular, we computed the PDFs by considering the full ensemble of particles up to the time T 22l c /v A and by setting the threshold E thr = 70%E max (blue) and E thr = 95%E max (red), where E max = 24.4 PeV (i.e. r g,max /l c = 2.2). Particles below the threshold show a distribution compatible with isotropy, with a mean pitchangle value of 1/2, reported in Figure 11 with a gray dashed line. On the other hand, the most energetic particles display a strongly anisotropic distribution, peaked at small µ loc , becoming more evident for larger thresholds. PDF E thr = 70%Emax E thr = 95%Emax Figure 11. PDFs of µ loc conditioned with the particle energy for the top 30% of the most energetic particles (blue) and for the top 5% (red), with the maximum energy being Emax = 24.4 PeV (rg,max/lc = 2.2). The horizontal dashed gray line displays the 1/2 value, corresponding to isotropic distribution. The exponential acceleration responsible for the high-energy tail is associated with a small pitch-angle cosine µ. Because we are not removing particles once they enter the exponential phase, they consequently undergo the standard pitch-angle process, thus producing the observed spreading in the µ loc distribution.
This confirms that the most-energetic particles preferentially have a small local pitch-angle cosine, i.e. that they move perpendicular to the local magnetic field. As we discuss in §5, this may be the very cause of the trapping: The velocity of the particles parallel to the local field is such that the particle stays in the structure long enough to be energized (while possibly drifting) by gradients in the magnetic field. During this process, the particle pitch angle changes gradually, as would be expected due to the quasi-conservation of the adiabatic invariant: Because the field changes along the trajectory (at most by the δB/B rms calculated at the scale of the Larmor radius of the particles) and v ⊥ p ⊥ /B must stay constant, then p ⊥ must change, hence the pitch angle changes (Voelk 1975) by about δB/B rms . This leads to a spread in the PDF of the pitch-angle distribution of the most energetic particles.

DISCUSSION OF THE PHYSICAL PROCESSES AT WORK
Here we discuss our results and their implications for astrophysical systems. It has been found that relativistic particles moving in stationary turbulent electric and magnetic fields, constituted by a snapshot in time of the incompressible MHD simulations, experience a complex energization process owing to the turbulent inductive electric field. The character of the process reveals at least a twofold nature. The secondorder, stochastic acceleration affecting the whole ensemble of particles goes hand in hand with a first-order mechanism, impacting only a few particles that are trapped in large-scale coherent structures of turbulence. These trapped particles experience exponential energy growth until the time when the gyroradius becomes comparable with the transverse size of the structure.
The phase of exponential increase in particle energy is of special interest for its potential implications for astrophysical systems. As discussed in previous sections, this stage is rather complex: It requires that a particle enter a region where two large-scale structures seem to be interacting, thereby leading to gradients in the induced electric field and most likely in the magnetic field as well. The number of particles showing evidence of this exponential energy increase is very small, only a few out of 10 5 (0.001%). All these particles seem to have a very small pitch-angle cosine, namely µ ∼ 0. They also exhibit a similar phenomenology in terms of characteristic times: The trapping time is quite large (∼ 5 − 10l c /v A ), and various oscillations, as described in the previous section, are recovered. The formation of such structure in a generic turbulence box is difficult to characterize in a quantitative way, also due to intermittency: The standard Kolmogorov scalings could be locally violated due to spatial inhomogeneities (Matthaeus et al. 2015). In particular, magnetic fluctuations could become weaker in large-scale coherent structures with respect to the globally averaged values because nonlinearities could be partially suppressed within large-scale structures, and perhaps elevated near flux-tube boundaries. In fact, as discussed earlier in this article, the structure in which the exponential acceleration takes place does not show evidence for special activity: it might have originated from a reconnection event, but it is not a site of reconnection. Indeed, reconnection plays no role in causing the rapid energy increase that we observe in the simulations. Figure 8 illustrates very clearly how complex the region is where the particle energy is seen to increase exponentially. The only properties that we can confidently associate to this region are (1) the presence of a rather organized large-scale flux tube that seems to be interacting with another structure, and (2) a gradient in the local electric field, due to plasma motion. Moreover, it is reasonable to speculate that there may be a gradient in the magnetic field due to the interaction between structures.
Given the complexity of the situation, it may be helpful to use a toy model to illustrate the different effects and check whether the qualitative picture is reproduced. The toy model gives us the opportunity to comment on the different physical processes at work.
We show the geometry of our toy model in Figure 12: The magnetic field B in that region is assumed to be oriented in theẑ direction, but it is assumed to have a gradient in the direction of the center of the tube, within a ring (green re-gion) of size L grad . We also assume that at least on a fraction of the surface of the tube, the plasma velocity has a gradient along thex direction (see zoom-in in the lower part of Figure  12). In the figure, this gradient is shown in the form of a sign reversal of the velocity but this is not required. We will comment below on the implications of a fluid velocity crossing the value u = 0.
Let us start by discussing the role of a gradient in the magnetic field, although the existence of a strong gradient does not emerge in an evident way from the simulations: A magnetic field, with or without gradients, cannot make work on charged particles. Hence, the exponential energy increase is certainly not related to such a gradient. On the other hand, the gradient introduces a drift, whose direction depends on the direction of the gradient. With reference to the situation illustrated in the lower part of Figure 12, this grad-B drift leads the particle to move in the y direction, namely to stay in the green-colored region. The drift velocity is proportional to the gradient: namely, the lower the gradient, the lower the drift velocity, so that the particle can stay longer in the island. Note that Equation (9) includes the grad-B drift, while the curvature drift is subdominant because we showed that the motion occurs at µ ∼ 0, i.e. p perp p. Inspired by the results of our simulations, we first assume that the particle has a very small pitch-angle cosine µ, namely, its motion is constrained to be in the x−y plane (blue curve in Figure 12). A small value of µ means that the particle can travel in theẑ direction with velocity ∼ cµ (all of the particles in the simulation are relativistic) and will eventually leave the island in a time ∼ L isl /cµ. For µ ∼ 1 (red curve in Figure 12), the escape time would be exceedingly fast and no appreciable acceleration can take place (see below).
At the same time that the particle rotates in the x − y plane, it is advected with the local plasma, which is expected to move at speeds close to the local Alfvén speed v A . The region of the gradient in the plasma speed (which corresponds to the gradient in the induced electric field) can then be crossed in a time L grad /v A ∼ ηL isl /v A , where we assumed that the gradient develops on a fraction η of the size of the island, L grad ∼ ηL isl . The time to escape the island alongẑ exceeds the time scale of advection if µη v A c . This apparently simple and rather constraining conclusion on the pitch angle of the particles in the acceleration region is in fact affected by several phenomena that can possibly enhance the trapping: one such phenomenon is diffusion, but it is hard to imagine that it may be effective. In fact, the particles that we are considering have an initial Larmor radius of only one order of magnitude smaller than the coherence scale, which in turn is of the same of order as the size of the island, L isl . Hence, the escape time is a fraction of l 2 c /(cl c /3) ∼ 3l c /c, where we assumed that the diffusion coefficient is not too far from that expected at r L ∼ l c (see Figure 3). The time estimated in this way is very close to the ballistic time scale (within an order of magnitude), too short to be of any relevance for CR trapping. Moreover, the fact that the pitch angle is close to zero makes diffusion ineffective, because the resonance condition would require modes with a very large wavenumber k, which is not accessible in our simulations, and not abundant in the standard turbulence spectrum anyway, especially after accounting for anisotropy in the cascade.
A second phenomenon that is instead expected to be more effective is associated with the existence of perturbations in the magnetic field along theẑ direction: Because on the scale of a few gyrations it is a good approximation to assume that the quantity v ⊥ p ⊥ /B is conserved, if there are fluctuations of order ∆ = δB(1/r L )/δB(1/L isl ) in the local field, the conservation of this quantity also implies that p ⊥ must be changing, namely µ must be changing. Oscillations in B imply oscillations in the pitch angle of magnitude ∝ ∆ 1/2 (Voelk 1975). Unfortunately, for the parameters of our simulation, the quantities ∆ 1/2 and v A /c are too close to identify this effect unambiguously. However, for more realistic values of the ratio v A /c, the effect of oscillations associated with longitudinal gradients in the magnetic field should dominate and force particle trapping provided that µ is smaller than ∆.
We finally come to the effect of the gradient in the plasma velocity. If there were a homogeneous plasma velocity u in thex direction, this would result in an induced electric field E y = u(x) c B, directed in theŷ direction. Notice that as long as the plasma speed is spatially constant, one can always move into a reference frame in which this velocity is absent: Hence, the presence of this induced electric field only introduces a drift velocity in thex direction and no net increase in the energy of the particles can exist. In other words this drift is only the E × B drift and leads to the particle advection with the background plasma in thex direction. The case where the velocity of the plasma is not constant in the x direction is more interesting. Notice that the simulations used here are incompressible, hence ∇ · u = 0. This does not contradict the assumption that a gradient du/dx = 0 exists, for instance, if two structures are moving against each other at u ∼ v A .
Below we show that the existence of this drift is the very source of particle acceleration. We will reach this conclusion in two independent ways.
The evolution in time of the particle distribution function f of the particles, under the effect of the du/dx gradient alone, can be written as: where ξ = ln(p), and we used du/dx ≈ v A /L grad . If the two structures are moving against each other, a more likely estimate for the gradient would be 2v A /L grad , but the difference is only quantitative, not qualitative. Using the method of characteristics, one gets The particle momentum is expected to grow exponentially due to the presence of a difference in velocity felt by particles during gyration around the magnetic field. One can look at this phenomenon in at least two independent ways: One way is to imagine that in each particle orbit the electric field in one half is not exactly compensated by the second half, hence there is a net electric field that can energize the particle. The second way to look at this as a first-order Fermi process, in which the particle bounces (each half-orbit) on a fluid moving with a different speed, qualitatively similar to what happens close to a shock front, where however the gradient is much larger. The acceleration has to come to an end when the gyroradius becomes comparable with the size of the trapping region unless more constraining phenomena occur on shorter time scales. Notice that the factor of 1/3 in Equation (10) is derived under the assumption of particle isotropy, which clearly does not apply here. We expect that the time scale of acceleration of particles with µ ∼ 0 must be somewhat shorter than the 3L grad /v A suggested by Equation (11).
In order to prove that our conclusion is solid, we also derive an analogous result in a more formal way, starting from the equation of motion of a particle moving in a setup as in Figure 12. In this derivation, we ignore the effect of particle drift due to the gradient of the magnetic field, hence where we used the expression for the local induced electric field derived above. Recalling now that p = mγv, where γ is the Lorentz factor of the particle, and multiplying Equation (12) by v x and Equation (13) by v y , we can deduce that Summing all terms and recalling some basic relations of relativistic kinematics, one obtains dp dt = q c Bv y u(x).
At this point we assume that the particle velocity is already relativistic, v c, and that the particle trajectory is weakly modified by the gradient in u(x). This means that the gradient is assumed to be weak on the scale of the Larmor gyration of the particles. In this case, Equation 16 can be averaged over many gyrations and recalling that x(t) ∝ sin(Ωt) and v y (t) ∝ sin(Ωt), the average leads to dp dt = 1 2 ap, where a ≈ v A /L grad is the gradient of u(x). It follows that This result is formally the same as in Equation (11), but it correctly shows that the time for exponential increase of the momentum is slightly shorter than expected for an isotropic particle distribution. In the absence of escape mechanisms, the trapping is expected to cease when the particle gyroradius becomes comparable with the island size, r L ∼ L isl . This provides a characteristic time for the duration of the energization process, where r L,0 = p 0 c/eB is the initial particle Larmor radius and Λ = ln L isl r L,0 . In order to have effective acceleration, the acceleration time must be shorter than all the time scales of escape parallel to B and due to drifts. By imposing that the time of escape alongẑ, τ cross,|| = L isl /cµ, be longer than the acceleration time, we obtain the constraint This result provides us with a simple explanation of the reason why only particles with a small value of µ actually show evidence of exponential momentum increase. The issue of drifts is more subtle: the E × B drift, as discussed above, simply leads to particle advection along the x direction. The region of the gradient is then crossed in a time of order τ esc,x ∼ L grad /v A , comparable to, though shorter than, the acceleration time. However, if the two islands move against each other with roughly the same speed, ∼ v A , the guiding center of the particle may be advected at speeds much smaller than v A , as a result of the fact that u(x) crosses zero.
The drift along theŷ direction is due to the gradient in the magnetic field and in principle can be very fast because the drift velocity is given by Equation (9), so that the time to drift out of the region where the gradient exists is of order where we assumed that r L ∼ 0.1L isl . For the parameters adopted in the simulation, this time is τ esc,y ∼ 0.3L isl /c, much shorter than all other times scales. In the presence of such an effect, no appreciable acceleration should be expected. The evidence of an exponential increase in particle energy suggests that either the gradient in the magnetic field is a small fraction of B/L grad , or that it is present all along the surface of the flux tube (green region in the top part of Figure 12), so that the particle drifts around the tube while retaining its pitch angle, as discussed above. It would be very useful to further investigate the phenomena discussed above, but because in our simulations v A = c/20, such time scales are too close to each other to reach a definitive conclusion on the hierarchy of drift timescales.
We conclude our in-depth discussion of the physical processes at work in the acceleration region by commenting on the role of reconnection. As we discussed in §4.3, our simulations suggest that the region where particle energization is fast is not very active: It may be the result of a reconnection event, but it is not a reconnection region. This is also to be expected: The current sheet is very thin compared with the Larmor radius of the relativistic particles considered here, hence the resistive electric fields cannot be responsible for particle acceleration. In this sense, a model like the one of de Gouveia dal Pino and Lazarian (2005), in which the escape from the region is regulated by the speed of the exhaust of the reconnection phenomenon, appears to be not well justified. On the other hand, Kowal et al. (2011) correctly pointed out that interacting islands away from reconnection events can energize particles, even exponentially. This result was, however, obtained in simulations that were optimized to create reconnection regions. We show that even in a generic 3D MHD simulation box there are structures that lead to the same physical phenomenon, and we provide a physical explanation of the phenomenon and a recipe of the conditions required for the acceleration to take place. It is also worth noticing that a first-order particle acceleration in converging islands in relativistic turbulence was also found by Comisso and Sironi (2018). The case of relativistic particles moving in a relativistic plasma v A ∼ c may be considerably different from the one described above.
The final part of this section is devoted to a discussion of the possible relevance of these phenomena for astrophysical turbulence, for instance, in the Galaxy as a whole.
In principle the exponential increase associated with phenomenon of particle trapping in interacting islands may be of great importance: For instance, if we take the turbulence in the Galaxy as our laboratory, one would expect v A ∼ 10 km s −1 and L isl ∼ few tens of parsec, with a typical magnetic field of 3 µG. If a CR particle enters one such structure and suffers an exponential increase in energy up to the point where the Larmor radius equals L isl , a maximum energy of ∼ 20 PeV would be reached, tantalizingly close to the energy of the knee. This simple numerical estimate stimulates some additional questions: What is the time scale for such acceleration? And what is the probability that a CR particle may encounter such a peculiar structure before escaping the Galaxy?
The time scale for exponential increase is as in Equation (19): One can see that for realistic values of the parameters, the acceleration time is ∼ (1 − 10)l c /v A , where we assumed that L isl ∼ l c . This seems to be in accordance with the numerical results shown in Figure 7 (top panel). The fact that this time is comparable to or exceeds the eddy turnaround time l c /v A is a source of concern because both the simulations and the toy model discussed above assume that the tur-bulence is static (the propagation of test particles was carried out in a snapshot at a given time of the MHD simulation).
For turbulence in the Galaxy, the acceleration time estimated above is of order ∼ 1 million years. For particles at the knee, the escape time from the Galaxy can be deduced from an extrapolation to high energies of the low-energy ( 1 TeV) confinement time as inferred from secondary/primary ratios (Evoli et al. 2019) and from the Be/B ratio, as discussed recently by , using AMS-02 data. A reasonable estimate for such escape time at E ∼ PeV is of ∼ 0.5 Myr, comparable with the acceleration time. The situation can be considered somewhat more promising if one thinks that lower-energy particles (with longer confinement time) are the ones required to be trapped in the interacting islands and eventually getting energized.
Assuming that in a time of the order of the confinement time a CR particle can probe the statistical properties of the turbulence, the question arises of how many particles can potentially interact with a region where trapping occurs and particle energy increases exponentially. In the simulations we ran, only a few particles out of the 10 5 (0.001%) experienced exponential energy increase. The probability of order 10 −5 ÷10 −4 that one of the particles, while carrying out a diffusive motion, may encounter a region of size ∼ l c in which there are the right conditions for trapping to occur is hard to quantify in that it is the convolution of turbulence properties (volume filling factor of islands that interact in the right way, in the presence of intermittency) and properties of the particle trajectory at the time of entering the island: As we discussed above, only particles that approach the region with very small µ can be trapped in the region for long enough to experience the exponential energy increase. Because diffusion isotropizes the particle distribution function (outside the island), the fraction of particles that at any point in the box have a pitch-angle cosine µ ∼ v A /c is ∼ 1 2 v A /c = 0.025, where the numerical value refers to the conditions of the simulation, while in the Galaxy that number would be ∼ 1.7 × 10 −5 . The fact that only a few particles out of 10 5 suffer exponential energy increase implies that, in the simulation, the filling factor of the island that allows such phenomenon is very small, of the order of ∼ 10 −3 .
A dedicated investigation of these issues would be most important as a future development of the concepts discussed in the present article and is crucial to assess the importance of these phenomena for particle acceleration in nature.

CONCLUSIONS
The use of high-resolution magnetohydrodynamic simulation data in conjunction with orbit calculations of a large number of charged test particles enables detailed examination of both spatial transport and energization in a generic 3D incompressible turbulence simulation not specifically de-vised to study reconnection. The limited number of scales that can be simulated in 3D at this time makes it difficult to go much beyond the state of the art in terms of investigating particle transport: Nevertheless we determined the diffusion coefficient of relativistic particles in a range of scales where resonances are numerically accessible in the simulation. This bounds us to about one decade in energy below the energy for which the Larmor radius equals the correlation length l c of the turbulence. At such energies, it is difficult to spot the effect of anisotropic cascade with respect to the local magnetic field, typical of MHD, although such anisotropy is certainly visible when studied in terms of statistical indicators (Matthaeus et al. 2012). In this sense, our results on the diffusion coefficient are compatible with those of Cohet and Marcowith (2016) but do not add to it. They are also compatible with the results previously obtained by Subedi et al. (2017) and Dundovic et al. (2020) using synthetic isotropic turbulence, rather than MHD turbulence. This latter result confirms that the effects of anisotropic cascades are not yet visible on the scales accessible to particles.
Our results on particle energization are much more interesting: We find that the bulk of the test particles simulated here are subject to a secular, second-order acceleration process, due to the interaction of the particles with the random electric fields induced through plasma motion in the simulation. A few out of the 10 5 test particles for which we simulate the trajectories happen to suffer very fast accelerationin fact exponential in time. This process is seen to end when the Larmor radius becomes comparable with the size of the magnetic structure in which the particles reside.
The second-order stochastic acceleration process has been analyzed by employing standard techniques including computation of the running diffusion coefficient in energy space. We find empirical values for energy diffusion D EE 0.01v A E 2 0 /l c , and for the associated characteristic time τ diff,E ∼ 10 2 l c /v A . These quantities are in qualitative agreement with naive expectations: The energy of a particle is expected to change in a random way (increase or decrease) due to the interaction with induced electric fields, so that their energy changes as ξv A c 2 E, with ξ 1, in a time that is approximately ∼ l c /c (see Figure 3). It follows that the diffusion coefficient in energy can be estimated as ξlc . For the v A /c = 1/20 adopted in our calculations D EE ≈ ξ 2 0.05E 2 /l c .
We characterize further the second-order acceleration process by evaluating the time-dependent probability distribution functions, for both particle energy and time increments of particle energy. These additional tests show the secondorder nature of the phenomenon.
As mentioned above, a few of the test particles in the simulation appear to have a quite different behavior, in that in addition to the slow second-order process, they happen to encounter special locations in the box where the energy is seen to increase exponentially, for times exceeding or about 10 Alfvén crossing times of the magnetic coherence length.
Closer examination reveals that these particles maintain a pitch angle not far from 90 • during this period, while their spatial trajectory is highly confined and differs greatly from its more typical random walk nature. In essence, these particles have become temporarily trapped within particular turbulence structures that are characterized by strong gradients. The structures are not directly associated with regions of activity in the plasma, and in fact, if any, they show less than normal activity, although visual inspection of these regions suggests the existence of extended interaction areas with at least another structure.
We built a toy model that includes in a simplified way the main drifts and gradients in the local electric fields. The toy model shows that the main reason for rapid particle energization is a sort of first-order Fermi process in which the energy of the particles grows due to a local gradient in the plasma velocity field. Both the temporal evolution of the energization process and the typical time scales for particle acceleration, as well as the maximum energy, are qualitatively reproduced in a correct way.