Martin Lemoine

Particle acceleration in strong MHD turbulence

Nonthermal acceleration of particles in magnetohydrodynamic (MHD) turbulence plays a central role in a wide variety of astrophysical sites. This physics is addressed here in the context of a strong turbulence, composed of coherent structures rather than waves, beyond the realm of quasilinear theory. The present description tracks the momentum of the particle through a sequence of frames in which the electric field vanishes, in the spirit of the original Fermi scenario. It connects the sources of energy gain (or loss) to the gradients of the velocity of the magnetic field lines, in particular the acceleration and the shear of their velocity flow projected along the field line direction, as well as their compression in the transverse plane. Those velocity gradients are subject to strong intermittency: they are spatially localized and their strengths obey power law distributions, as demonstrated through direct measurements in the incompressible MHD simulation of the Johns Hopkins University database. This intermittency impacts the acceleration process in a significant way, which opens up prospects for a rich phenomenology. In particular, the momentum distribution, which is here captured through an analytical random walk model, displays extended power law tails with soft-to-hard evolution in time, in general agreement with recent kinetic numerical simulations. Extensions to this description and possible avenues of exploration are discussed.

I. INTRODUCTION

The remarkable wealth of physical phenomena that inhabit the cascades of magnetohydrodynamic (MHD) turbulence [START_REF] Beresnyak | Turbulence in Magnetohydrodynamics[END_REF][START_REF] Schekochihin | MHD turbulence: a biased review[END_REF][START_REF] Matthaeus | Turbulence in space plasmas: Who needs it?[END_REF], from the large scales of stirring motions down to the microscopic dissipative layers, the ubiquity of magnetized turbulence in space plasmas, and the decisive roles that it plays in various astrophysical settings [START_REF] Balbus | Instability, turbulence, and enhanced transport in accretion disks[END_REF][START_REF] Elmegreen | Interstellar turbulence I: observations and processes[END_REF][START_REF] Chen | Recent progress in astrophysical plasma turbulence from solar wind observations[END_REF][START_REF] Mohapatra | Turbulence in the intracluster medium: Simulations, observables, and thermodynamics[END_REF], have turned its study into a field of research of its own right, with broad interdisciplinary connections. In highly conducting, collisionless astrophysical plasmas, MHD turbulence indeed governs the transport of particles, including cosmic rays, it regulates the thermodynamical properties of plasmas, and of more direct concern to the present study, it promotes particle acceleration, thereby opening a connection between the fundamentals of plasma turbulence and nonthermal astronomy.

In his seminal works [START_REF] Fermi | On the Origin of the cosmic radiation[END_REF][START_REF] Fermi | Galactic magnetic fields and the origin of cosmic radiation[END_REF], E. Fermi described particle acceleration in highly conducting, turbulent plasmas through the kinematics of discrete scattering events: a bounce of the particle on a region of high magnetic field strength in motion, or a turn-around of the trajectory in a curved magnetic field configuration, which, in a more modern context of guiding center theories of particle motion, can be pictured as "grad-B drift" and "curvature drift" energization.

Most of the literature in this field nevertheless hinges on the notion of wave particle resonance in the context of quasilinear theory [START_REF] Kennel | Velocity space diffusion from weak plasma turbulence in a magnetic field[END_REF][START_REF] Hall | Diffusion, scattering, and acceleration of particles by stochastic electromagnetic fields[END_REF][START_REF] Lerche | Quasilinear theory of resonant diffusion in a magneto-Active, relativistic plasma[END_REF][START_REF] Kulsrud | The relativistic quasilinear theory of particle acceleration by hydromagnetic turbulence[END_REF] and its more recent nonlinear extensions [START_REF] Völk | Nonlinear perturbation theory for cosmic ray propagation in random magnetic fields[END_REF][START_REF] Matthaeus | Nonlinear collisionless perpendicular diffusion of charged particles[END_REF][START_REF] Yan | Cosmic-Ray propagation: non-linear diffusion parallel and perpendicular to mean magnetic field[END_REF][START_REF] Shalchi | Nonlinear Cosmic Ray Diffusion Theories[END_REF]. In this picture, the turbulent plasma is described as a linear superposition of uncorrelated waves, e.g. Alfvén and magnetosonic modes in ideal magnetohydrodynamics (MHD). Energy gain then derives from phase locking of the particle trajectory with those waves. A truly remarkable aspect of quasilinear theory is its ability to provide explicit calculations of the acceleration process, up to the very shape of the accelerated spectrum, from a limited number of parameters, most notably the turbulence strength and the power spectrum of magnetic fluctuations [START_REF] Schlickeiser | An explanation of abrupt cutoffs in the optical-infrared spectra of non-thermal sources -A new pile-up mechanism for relativistic electron spectra[END_REF][START_REF] Achatz | Charged particle transport in a turbulent magnetized plasma -A reassessment and extension of quasilinear theory[END_REF][START_REF] Jaekel | The Fokker-Planck coefficients of cosmic ray transport in random electromagnetic fields[END_REF][START_REF] Schlickeiser | Quasi-linear theory of cosmic-ray transport and acceleration: The role of oblique magnetohydrodynamic waves and transit-time damping[END_REF][START_REF] Schlickeiser | Cosmic Ray Astrophysics[END_REF][START_REF] Stawarz | On the momentum diffusion of radiating ultrarelativistic electrons in a turbulent magnetic field[END_REF][START_REF] Schlickeiser | A new cosmic ray transport theory in partially turbulent space plasmas: Extending the quasilinear approach[END_REF]. This has triggered a wealth of applications in astrophysics, on all scales, from the Sun [START_REF] Miller | Stochastic electron acceleration by cascading fast mode waves in impulsive solar flares[END_REF][START_REF] Petrosian | Stochastic acceleration of electrons and protons. I. Acceleration by parallelpropagating waves[END_REF][START_REF] Selkowitz | Stochastic Fermi acceleration of subrelativistic electrons and its role in impulsive solar flares[END_REF][START_REF] Yan | Particle acceleration by fast modes in solar flares[END_REF][START_REF] Bian | A classifica-tion scheme for turbulent acceleration processes in solar flares[END_REF] to the remote Universe [START_REF] Liu | Electron acceleration around the supermassive black hole at the Galactic center[END_REF][START_REF] Brunetti | Compressible turbulence in galaxy clusters: Physics and stochastic particle reacceleration[END_REF][START_REF] Petrosian | Heating and acceleration of intracluster medium electrons by turbulence[END_REF][START_REF] Brunetti | Acceleration of primary and secondary particles in galaxy clusters by compressible MHD turbulence: From radio haloes to gammarays[END_REF][START_REF] Brunetti | Stochastic reacceleration of relativistic electrons by turbulent reconnection: A mechanism for cluster-scale radio emission?[END_REF][START_REF] Eckert | On the connection between turbulent motions and particle acceleration in galaxy clusters[END_REF], including high-energy sources harboring relativistic plasmas in extreme environments [START_REF] Bykov | Electron acceleration and efficiency in nonthermal gamma-ray sources[END_REF][START_REF] Tramacere | Stochastic acceleration and the evolution of spectral distributions in synchro-self-Compton sources: A self-consistent modeling of blazars' flares[END_REF][START_REF] Lemoine | A corrugated termination shock in pulsar wind nebulae?[END_REF][START_REF] Asano | Ultrahigh-energy cosmic ray production by turbulence in gamma-ray burst jets and cosmogenic neutrinos[END_REF][START_REF] Xu | Adiabatic non-resonant acceleration in magnetic turbulence and hard spectra of gammaray bursts[END_REF][START_REF] Asano | Blazar spectra with hardsphere-like acceleration of electrons[END_REF][START_REF] Xu | On the broadband synchrotron spectra of pulsar wind nebulae[END_REF].

However, the realization [START_REF] Chandran | Scattering of energetic particles by anisotropic magnetohydrodynamic turbulence with a Goldreich-Sridhar power spectrum[END_REF][START_REF] Yan | Scattering of cosmic rays by magnetohydrodynamic interstellar turbulence[END_REF][START_REF] Shalchi | Cosmic ray transport in anisotropic magnetohydrodynamic turbulence. III. Mixed magnetosonic and Alfvénic turbulence[END_REF] that wave particle resonances disappear in the context of modern anisotropic MHD turbulence theories [START_REF] Goldreich | Toward a theory of interstellar turbulence. 2: Strong alfvenic turbulence[END_REF][START_REF] Goldreich | Magnetohydrodynamic turbulence revisited[END_REF][START_REF] Cho | The anisotropy of magnetohydrodynamic Alfvénic turbulence[END_REF][START_REF] Maron | Simulations of incompressible magnetohydrodynamic turbulence[END_REF][START_REF] Cho | Compressible magnetohydrodynamic turbulence: Mode coupling, scaling relations, anisotropy, viscosity-damped regime and astrophysical implications[END_REF][START_REF] Boldyrev | Spectrum of magnetohydrodynamic turbulence[END_REF][START_REF] Beresnyak | Spectral slope and Kolmogorov constant of MHD turbulence[END_REF] has cast shadows on this program, and more generally on the role of waves in turbulent dissipation. Accounting for nonlinear effects, such as resonance broadening related to the finite lifetime of the eigenmodes, can at best lead to scattering timescales of the order of the coherence scale of the turbulence, at least for for those modes that are subject to anisotropy. Consequently, this renders hazardous the extrapolation of predictions based on quasilinear calculations to the regime of strong turbulence1 ; this point will be detailed in the following Sec. II.

Meanwhile, in situ measurements in the solar wind [START_REF] Kiyani | Global scale-invariant dissipation in collisionless plasma turbulence[END_REF][START_REF] Osman | Evidence for Inhomogeneous Heating in the Solar Wind[END_REF][START_REF] Osman | Intermittency and local heating in the solar wind[END_REF][START_REF] Bruno | Intermittency in solar wind turbulence from fluid to kinetic scales[END_REF] as well as numerical simulations [START_REF] Wan | Intermittent dissipation at kinetic scales in collisionless plasma turbulence[END_REF][START_REF] Karimabadi | Coherent structures, intermittent turbulence, and dissipation in high-temperature plasmas[END_REF][START_REF] Osman | Magnetic reconnection and intermittent turbulence in the solar wind[END_REF] have revealed that dissipation is closely associated to coherent structures with sharp gradients, e.g. current sheets and flux tubes. It has become apparent that MHD turbulence must be described as a collection of structures rather than linear uncorrelated waves, or possibly as a combination of the two [START_REF] Grošelj | Kinetic turbulence in astrophysical plasmas: waves and/or structures?[END_REF][START_REF] Oughton | Solar wind turbulence: Connections with energetic particles[END_REF].

Recent numerical simulations suggest that particles can draw energy from their interaction with the coherent structures of a turbulent bath [START_REF] Pecora | Ion diffusion and acceleration in plasma turbulence[END_REF][START_REF] Trotta | Fast acceleration of transrelativistic electrons in astrophysical turbulence[END_REF], and it is known that particles exploring a sheared velocity field with a finite mean free path can gain energy in a nonresonant fashion [START_REF] Brunetti | Stochastic reacceleration of relativistic electrons by turbulent reconnection: A mechanism for cluster-scale radio emission?[END_REF][START_REF] Xu | Adiabatic non-resonant acceleration in magnetic turbulence and hard spectra of gammaray bursts[END_REF][START_REF] Bykov | Cosmic rays acceleration by large scale compressible and incompressible motions of plasma[END_REF][START_REF] Ptuskin | Cosmic-ray acceleration by long-wave turbulence[END_REF][START_REF] Bykov | Theory of charge particle acceleration by a shock wave ensemble in a turbulent medium[END_REF][START_REF] Webb | Diffusive-compression acceleration and turbulent diffusion of cosmic rays in quasi-periodic and turbulent flows[END_REF][START_REF] Chandran | Acceleration of energetic particles by large-scale compressible magnetohydrodynamic turbulence[END_REF][START_REF] Jokipii | Compression acceleration in astrophysical Plasmas and the production of f (v) ∝ v -5 spectra in the heliosphere[END_REF][START_REF] Ohira | Turbulent shear acceleration[END_REF][START_REF] Zank | Particle acceleration via reconnection processes in the supersonic solar wind[END_REF][START_REF] Le Roux | A kinetic transport theory for particle acceleration and transport in regions of multiple contracting and reconnecting inertial-scale flux ropes[END_REF][START_REF] Beresnyak | First-order particle acceleration in magnetically-driven flows[END_REF][START_REF] Lemoine | Generalized Fermi acceleration[END_REF]. However, a definite model able to relate the physics of particle acceleration with the characteristics of the structures of a turbulent bath remains lacking.

The aim of this paper is to make progress along those lines, in particular to study the physics of particle acceleration in strong turbulence, through nonresonant interactions with random velocity structures. To do so, it relies on a formalism that connects the sources of energy gain/loss to the shear of the E × B velocity field. That formalism was introduced in Ref. [START_REF] Lemoine | Generalized Fermi acceleration[END_REF], which considered the case of a particle exploring in an isotropic manner a random flow characterized by a single length scale. In a turbulent context, this can represent a particle whose gyroradius r g is comparable or greater than the coherence scale c of the turbulent bath. By contrast, the present paper focuses on particles whose gyroradius is much smaller than this length scale c , i.e. deep in the inertial range. It thus pays attention to the anisotropy imposed by the magnetic field line, to the approximate gyromotion of the particle around that field line, and it considers the influence of all modes of scales larger than r g .

As will be made explicit, such nonresonant acceleration has important consequences which distinguishes it from linear wave particle interactions. Notably, particle energization can be seen as a form of shear acceleration; hence energization scales with the gradients of the magnetic energy density, rather than the turbulent magnetic energy density. More importantly, the physics of acceleration becomes strongly affected by intermittency because the shear of the velocity field peaks on small scales. Different particles thus experience vastly different histories, which gives rise to power law tails of the momentum distribution function. This effect will be made manifest through an analytical random walk model that captures the effect of an inhomogeneous repartition of structures.

While most of the discussion is analytical in nature, contact will be made with large-scale simulations of turbulence for specific points. The statistics of the random velocity structures, in particular, will be extracted from the forced MHD simulation of the Johns Hopkins Turbulence database2 , a 1024 3 simulation of driven, incompressible MHD simulation, hereafter referred to as the JHU-MHD simulation [START_REF] Li | A public turbulence database cluster and applications to study Lagrangian evolution of velocity increments in turbulence[END_REF][START_REF] Eyink | Fluxfreezing breakdown in high-conductivity magnetohydrodynamic turbulence[END_REF].

The discussion is ordered as follows. Section II offers a more in-depth perspective on the difference between resonant wave particle interactions and nonresonant acceleration in a random velocity flow. Section III sets up and describes the present model of nonresonant acceleration in a bath of strong turbulence. Section IV provides a quantitative assessment of particle energization in the present framework; in particular, it evaluates the rates of mean and diffusive energy gains, and it calculates the spectral shape in the frame of a microscopic picture that considers the influence of turbulence intermittency. A summary and conclusions are drawn in Sec. V.

II. MOTIVATIONS: WAVE PARTICLE INTERACTIONS IN STRONG TURBULENCE

Whether collisionless MHD turbulence should be described as a bath of waves, or structures, is a longstanding debate, see [START_REF] Grošelj | Kinetic turbulence in astrophysical plasmas: waves and/or structures?[END_REF][START_REF] Oughton | Solar wind turbulence: Connections with energetic particles[END_REF] for recent discussions and comparison to simulation and observational data. The general point of view of the present study is that, even if wave packets are present on some scales l in a bath of strong turbulence, they more likely appear as random velocity structures of ill-defined polarization to particles of gyroradius r g l, than as linear wave packets. This view relies on the following arguments.

Quasilinear theory (QLT) has been devised for weak turbulence, where the fluctuations are assumed to be much smaller in magnitude than the mean quantities. It also rests on the premise that the turbulent bath can be described as the sum of noninteracting linear waves, and that the particle collects their influence as it propagates unperturbed along the mean field [START_REF] Kennel | Velocity space diffusion from weak plasma turbulence in a magnetic field[END_REF][START_REF] Hall | Diffusion, scattering, and acceleration of particles by stochastic electromagnetic fields[END_REF][START_REF] Lerche | Quasilinear theory of resonant diffusion in a magneto-Active, relativistic plasma[END_REF]. It has since been enlarged and generalized to account for interactions with wave packets of finite temporal coherence, and for deviations of the particle trajectory from its unperturbed orbit, in the frame of so-called extended quasilinear calculations or nonlinear guiding center theories [START_REF] Völk | Nonlinear perturbation theory for cosmic ray propagation in random magnetic fields[END_REF][START_REF] Matthaeus | Nonlinear collisionless perpendicular diffusion of charged particles[END_REF][START_REF] Yan | Cosmic-Ray propagation: non-linear diffusion parallel and perpendicular to mean magnetic field[END_REF][START_REF] Shalchi | Nonlinear Cosmic Ray Diffusion Theories[END_REF].

The extrapolation of quasilinear theory to the realm of large-amplitude turbulence -δB ∼ B -faces essential difficulties, however. A first one is that one can no longer regard the waves as uncorrelated and non-interacting, which invalidates the random phase approximation central to QLT calculations. This is notably illustrated in Ref. [START_REF] Maron | Simulations of incompressible magnetohydrodynamic turbulence[END_REF] and for what regards acceleration, in Ref. [START_REF] Arzner | The effect of coherent structures on stochastic acceleration in MHD turbulence[END_REF]. A second one is related to critical balance [START_REF] Chandran | Scattering of energetic particles by anisotropic magnetohydrodynamic turbulence with a Goldreich-Sridhar power spectrum[END_REF][START_REF] Yan | Scattering of cosmic rays by magnetohydrodynamic interstellar turbulence[END_REF]: eddy anisotropy implies that the interaction between waves and particles becomes a nonresonant phenomenon 3 , to such a degree that all scales r g contribute equally to the diffusion of particles (see below). However, the mean 3 Fast magnetosonic modes seem to preserve an isotropic cascade and thus their efficiency in scattering particles through gyroresonant and transit-time damping resonances [START_REF] Yan | Scattering of cosmic rays by magnetohydrodynamic interstellar turbulence[END_REF][START_REF] Cho | Compressible magnetohydrodynamic turbulence: Mode coupling, scaling relations, anisotropy, viscosity-damped regime and astrophysical implications[END_REF][START_REF] Cho | Particle acceleration by magnetohydrodynamic turbulence[END_REF]. If these modes survive in large-amplitude turbulence, they may play an important role in particle scattering and acceleration. Note that transit-time damping is in itself a form of nonresonant acceleration [START_REF] Teraki | Particle energy diffusion in linear magnetohydrodynamic waves[END_REF], and as such, it will be captured in the nonresonant formalism described in this paper.

field direction on a given scale does not coincide with that seen from another scale, because in strong turbulence, this mean field is defined through spatial coarsegraining, scale by scale [START_REF] Cho | The anisotropy of magnetohydrodynamic Alfvénic turbulence[END_REF]. Consequently, the polarization of a wave on some scale l 1 becomes lost if viewed from another smaller scale l 2 l 1 , which thus invalidates the predictions of QLT as extrapolated to the regime of strong turbulence. This point is made more explicit in the paragraphs that follow.

To understand why all scales contribute to particle energization in the quasilinear view of anisotropic turbulence, consider the QLT predictions for the momentum diffusion coefficient,

D pp ∼ β 2 A e 2 dk ⊥ dk k ⊥ R k J 1 (z ⊥ ) 2 z 2 ⊥ S k , (1)
where only the dominant term has been retained, and numerical prefactors of the order of unity have been ignored; β A ≡ v A /c, v A the Alfvén velocity. A detailed expression for the above can be found in, e.g. [START_REF] Demidem | Particle acceleration in relativistic turbulence: A theoretical appraisal[END_REF]. The wave numbers k ⊥ and k are understood to be oriented perpendicular and parallel to the (scaledependent) mean field, while the function S k represents the power spectrum. For anisotropic turbulence [START_REF] Goldreich | Toward a theory of interstellar turbulence. 2: Strong alfvenic turbulence[END_REF][START_REF] Goldreich | Magnetohydrodynamic turbulence revisited[END_REF][START_REF] Cho | The anisotropy of magnetohydrodynamic Alfvénic turbulence[END_REF][START_REF] Maron | Simulations of incompressible magnetohydrodynamic turbulence[END_REF][START_REF] Cho | Compressible magnetohydrodynamic turbulence: Mode coupling, scaling relations, anisotropy, viscosity-damped regime and astrophysical implications[END_REF][START_REF] Boldyrev | Spectrum of magnetohydrodynamic turbulence[END_REF][START_REF] Beresnyak | Spectral slope and Kolmogorov constant of MHD turbulence[END_REF],

S k ∝ δB 2 k -10/3 ⊥ k 1/3 min g k / k 2/3 ⊥ k 1/3 min
, where g(x) is a function peaked around 0 and spread around [-1, +1], of integral unity, which characterizes the anisotropy. J 1 (z ⊥ ) is a Bessel function, with z ⊥ = k ⊥ r g⊥ , r g⊥ the particle Larmor radius. The quantity of interest, here, is the resonance function R k , which characterizes the strength/duration of the interaction between waves and particles. Accounting for wave damping at a rate | ω| ∼ k v A broadens the resonance and confers to R k a Lorentzian or Gaussian shape. Performing the integration over k then brings the integral under the form

D pp ∝ p 2 k min c +∞ ρ dz ⊥ z -3 ⊥ J 1 (z ⊥ ) 2 , (2)
and ρ ≡ k min r g (k min the minimum wave number of the cascade) is assumed to be much smaller than unity, for particles interacting with modes in the inertial range of the turbulent spectrum. Most of the contribution comes from the interval ρ ≤ z ⊥ ≤ 1, corresponding to spatial scales larger than r g , yet smaller than c ≡ 2π/k min , the integral scale of the turbulent spectrum. The small argument limit of the Bessel function implies that the integrand scales as z -1 ⊥ . This confirms the above claim that, in this QLT description of anisotropic turbulence, all modes of scale larger than r g contribute equally to the scattering of the particle. Now, the anisotropy of the cascade constrains the angle θ k between the direction of the longitudinal axis of an eddy and that of the mean field to θ k k /k.

On small scales, meaning large k, k ∝ k

2/3 ⊥ implies θ k ∼ (k ⊥ /k min ) -1/3 . However, on scale l ∼ 1/k ∼ 1/k ⊥ , the field direction itself undergoes angular excursions of magnitude δθ ∼ δB l /B ∼ η 1/2 (k ⊥ /k min ) -1/3
, where η = δB 2 /B 2 denotes the ratio of the turbulent magnetic energy density to the total magnetic energy density and δB

2 l ∼ k δB 2 k | k∼1/l ∼ k 2 ⊥ k S k | k∼1/l .
Consequently, in large-amplitude turbulence, η ∼ 1, those excursions are, on all scales, comparable to θ k and the wave polarization is ill-defined. As anticipated above, this invalidates the QLT calculation, which rests on Eq. [START_REF] Beresnyak | Turbulence in Magnetohydrodynamics[END_REF].

The above compares the directions of the magnetic field after coarse-graining on two different scales, but at the same physical point, whereas the particle travels along the line as it cumulates the effect of waves. Given the intrinsic stochasticity of magnetic field lines in strong turbulence [START_REF] Eyink | Fluxfreezing breakdown in high-conductivity magnetohydrodynamic turbulence[END_REF], the correlation between the direction of the magnetic field, scale to scale, is bound to be lost even more rapidly. Small fluctuations in the magnetic field line direction can be incorporated within the frame of extended quasilinear theories as an additional resonance broadening attached to the randomization of the particle pitch angle, see e.g. [START_REF] Demidem | Particle acceleration in relativistic turbulence: A theoretical appraisal[END_REF][START_REF] Voelk | Cosmic ray propagation in interplanetary space[END_REF][START_REF] Xu | Resonance-broadened transit time damping of particles in MHD turbulence[END_REF], but this does not change the above conclusions.

Yet another implication of strong turbulence is that a wave packet on a given scale l 1 is nonlinearly distorted on a timescale comparable to its wave period [START_REF] Maron | Simulations of incompressible magnetohydrodynamic turbulence[END_REF], hence during the interaction itself it loses its polarization and characteristics of linear waves.

In summary, therefore, from the point of view of a particle propagating in a strongly turbulent environment, large-scale wave packets -assuming that they represent a fair representation of the turbulence -are more likely seen as imprinting velocity and magnetic fluctuations, rather than as waves contributing to wave particle interactions that are anyhow nonresonant given their inherent anisotropy.

Numerical simulations provide essential tools to test such models, but their interpretation deserves special care. In Ref. [START_REF] Lynn | Acceleration of relativistic electrons by magnetohydrodynamic turbulence: Implications for nonthermal emission from black hole accretion disks[END_REF], for instance, it was shown that the above predictions of (extended) quasilinear theory, viz. a momentum diffusion coefficient D pp ∼ v 2 A p 2 /( c c), could account for the energization of test-particles in MHD simulations of mildly sub-Alfvénic incompressible turbulence. However, that agreement cannot be interpreted as a confirmation of quasilinear predictions because the scaling of that diffusion coefficient is, as a matter of fact, generic to any nonresonant acceleration mechanism. The diffusion coefficient indeed scales, for dimensional reasons, as D pp ∝ p 2 /t int , where t int defines some interaction timescale. In the absence of resonant interactions, the integral length scale c sets this interaction time, viz. t int ∼ c /v (with v the particle velocity), while actual resonances would introduce another scale in the problem (r g ) and lead to the estimate t int ∝ r 2-q g , with q the index of the power spectrum of magnetic fluctuations [START_REF] Schlickeiser | Cosmic Ray Astrophysics[END_REF]. The nonresonant acceleration mechanisms discussed thereafter, which do not involve waves but velocity structures, will also yield estimates that match those numerical simulations. Finally, recent kinetic numerical simulations have contradicted, in a rather direct way, the generic predictions of a quasilinear picture of a homogeneous bath of waves: Refs. [START_REF] Trotta | Fast acceleration of transrelativistic electrons in astrophysical turbulence[END_REF][START_REF] Nättilä | Radiative turbulent flares in magnetically-dominated plasmas[END_REF] have observed that particles can be accelerated in localized regions, while Refs. [START_REF] Zhdankin | Kinetic turbulence in relativistic plasma: From thermal bath to nonthermal continuum[END_REF][START_REF] Zhdankin | System-size convergence of nonthermal particle acceleration in relativistic plasma turbulence[END_REF][START_REF] Zhdankin | Numerical investigation of kinetic turbulence in relativistic pair plasmas -I. Turbulence statistics[END_REF][START_REF] Comisso | Particle acceleration in relativistic plasma turbulence[END_REF][START_REF] Wong | First-principles demonstration of diffusive-advective particle acceleration in kinetic simulations of relativistic plasma turbulence[END_REF][START_REF] Comisso | The interplay of magnetically dominated turbulence and magnetic reconnection in producing nonthermal particles[END_REF][START_REF] Zhdankin | Electron and ion energization in relativistic plasma turbulence[END_REF][START_REF] Nättilä | Runko: Modern multi-physics toolbox for simulating plasma[END_REF] have demonstrated that the particle distribution function exhibit non-thermal power law tails, in stark contrast with the standard lognormal spectrum predicted by QLT. Such power law spectra appear to be related to the intermittent nature of the acceleration process [START_REF] Lemoine | Power-law spectra from stochastic acceleration[END_REF]. The prevalent role of intermittency and the development of power law spectra will emerge as natural consequences of the nonresonant acceleration scenario that is discussed next.

III. PARTICLE ACCELERATION IN RANDOM VELOCITY STRUCTURES

A. General framework

Within the frame of ideal MHD, a particle can gain energy and diffuse in momentum space as it explores the shear of the motional electric field in the plasma. It is the shear, rather than the magnitude of the electric field that matters, because in any region devoid of shear, one can boost to a frame in which this electric field vanishes, and there the particle undergoes helicoidal motion around the magnetic field, at constant energy.

To describe such random interactions, it proves convenient to define the frame R / E in which this motional electric field vanishes. With respect to the lab-frame, R / E moves at velocity (three-velocity, measured in units of c) β E ≡ E × B/B 2 (Gaussian units are adopted all throughout). More specifically, following Ref. [START_REF] Lemoine | Generalized Fermi acceleration[END_REF], the idea is to track the motion of the particle in physical space in the lab-frame, and the history of its momentum in this R / E frame. This provides a generalized picture of the classical Fermi mechanism, which is usually described as elastic interactions in the frame of pointlike, discrete magnetized scattering centers. The present picture is well suited to describe particle transport in a continuous random flow of velocity structures.

An important difference with respect to quasilinear theory is that the calculation for energy gain is performed in this frame R / E , not in the lab frame. For subrelativistic turbulence, in which the rms velocity fluctuation δu 2 p 1/2 c, this amounts to a simple velocity shift. Nevertheless, it will be shown in the following that the statistics of the field line velocity field u E differ from those of the plasma velocity.

The formalism that follows is fully relativistic because relativity offers the requisite tools to describe the history of particles through this continuous sequence of (noninertial) reference frames. This relativistic description is furthermore mandatory to make contact with kinetic numerical simulations of relativistic turbulence, which have revealed broad momentum spectra of accelerated particles [START_REF] Zhdankin | Kinetic turbulence in relativistic plasma: From thermal bath to nonthermal continuum[END_REF][START_REF] Zhdankin | System-size convergence of nonthermal particle acceleration in relativistic plasma turbulence[END_REF][START_REF] Zhdankin | Numerical investigation of kinetic turbulence in relativistic pair plasmas -I. Turbulence statistics[END_REF][START_REF] Comisso | Particle acceleration in relativistic plasma turbulence[END_REF][START_REF] Wong | First-principles demonstration of diffusive-advective particle acceleration in kinetic simulations of relativistic plasma turbulence[END_REF][START_REF] Comisso | The interplay of magnetically dominated turbulence and magnetic reconnection in producing nonthermal particles[END_REF][START_REF] Zhdankin | Electron and ion energization in relativistic plasma turbulence[END_REF][START_REF] Nättilä | Runko: Modern multi-physics toolbox for simulating plasma[END_REF], and which thus offer ideal experiments for the present model. Subrelativistic approximations will nevertheless be provided where opportune.

As β E is a random quantity, the R / E frame is not globally inertial, but the electric field is zero everywhere along the particle trajectory. Consequently, energy gains and losses are entirely expressed as inertial forces felt by the particle as the R / E frame accelerates or decelerates. One clear advantage of this point of view is to relate in a direct way the statistics of energy gain to those of the velocity field.

A key ingredient is how the particle sees, or travels through, the velocity structures. While Ref. [START_REF] Lemoine | Generalized Fermi acceleration[END_REF] assumed that the particle explored in an isotropic manner a random velocity field whose variations were characterized by a single length scale, the present work considers the more realistic case of a particle subject to a crossed electricmagnetic field configuration with perturbations spanning a broad range of scales. The former study thus applies to particles whose gyroradius is comparable or larger than the coherence scale c , while the present generalizes it to particles in the inertial range, that is r g c . Those particles are sensitive to fluctuations on all scales larger than their gyroradius as they gyrate around the local mean field. On the other hand, fluctuations on scales smaller than r g can be neglected here, as commonly done, because they carry a small fraction of the power if r g c and because their influence averages out over a gyroperiod.

The acceleration of particles is thus described here as a continuous sequence of interactions with moving magnetized structures on scales larger than the gyroradius. This implies that all quantities -velocity flow, magnetic field in particular -must be understood as coarse-grained over the scale r g . Formally, this means that we consider a velocity field

β E l = E l × B l /B 2 l
, where E l and B l are as before defined through coarse-graining on scale l, which differs from the "true" β E by O δB 2 l /B 2 . Note that the field entering the "true" contribution is obtained without coarse-graining here, although this limit is somewhat ambiguous in the frame of an MHD description, which breaks down on kinetic scales (see [START_REF] Eyink | Cascades and dissipative anomalies in nearly collisionless plasma turbulence[END_REF], in particular). Let us stress at this point that ideal MHD is not a requisite of the present approach, which only needs to assume that parallel electric fields play a negligible role in particle acceleration. This means, in particular, that the electric field characterizing the velocity field β E can incorporate other terms of Ohm's law beyond the ideal one, such as, for instance, the Hall MHD term.

For the sake of clarity and simplicity, all indices l referring to coarse-graining will be suppressed when not necessary. All quantities expressed in the frame R / E will be primed. Furthermore, all four-velocities u in R / E are from now on expressed in units of c, so that u = p /(mc) in terms of the momentum p , and similarly for the fourvelocity of the magnetic field lines,

u E = β E / √ 1 -u E 2 (unless otherwise noted).
Given the above scheme of approximation, the R / E frame is mainly characterized by two four-vectors: u E α , the (timelike) four-velocity of R / E , and b α , the unit spacelike four-vector indicating the direction of the magnetic field line. This construction is made more explicit in Appendix A, which also provides the details of the calculations.

In R / E , the particle Lorentz factor γ and four-velocity component parallel to the mean field, u , are then found to evolve as follows:

1 c dγ dτ = -u b u c e β b e γ c ∂ ∂x γ u E β , 1 c du dτ = -u b u c e β b e γ c ∂ ∂x γ b β . (3)
In this equation, u b , with b ∈ {0, . . . , 3}, represents a component of the particle four-velocity and e β b is a component of the vierbein that connects lab-frame quantities to quantities evaluated in R / E , i.e. which provides the instantaneous Lorentz transform between those two frames. Note, in particular, that the spatiotemporal derivative ∂ γ ≡ ∂/∂x γ in the above equation is expressed in labframe quantities. On the other hand, the derivative on the lhs is taken with respect to proper time τ which parametrizes the particle trajectory. The connection between the statistics of the velocity flow of the magnetic field lines and particle energization is manifest in Eq. ( 3).

On scales larger than r g , we can assume that, at each point along the trajectory, the particle undergoes a helicoidal orbit around the (curved) magnetic field line. In this case, the above equation can be further simplified, see Appendix A,

1 c dγ dτ = -γ u a E • b -u 2 Θ - 1 2 u ⊥ 2 Θ ⊥ , 1 c du dτ = -γ 2 a E • b -γ u Θ - 1 2 u ⊥ 2 b • ∇ ln B , (4)
introducing the components of the shear of the field line four-velocity, parallel (Θ ) and perpendicular (Θ ⊥ ) to the magnetic field line, as well as its acceleration a E ,

Θ ≡ b α b β ∂ α u E β , Θ ⊥ ≡ η αβ -b α b β ∂ α u E β , a E α ≡ u E β ∂ β u E α . (5)
The first line of Eq. ( 4) represents the main equation governing the acceleration of particles in random velocity structures in our current scheme of approximation. The force contributions that enter this equation find a natural interpretation. The first term ∝ a E • b describes the effective gravity that the particle feels in the R / E frame in the direction of the magnetic field line as this line accelerates or decelerates. The second term ∝ Θ captures energy gain or loss through shear acceleration along the direction of motion of the particle (along the field line). Similarly, the particle can gain (or loose) energy if the field lines, to which it is attached, are compressed (or dilated) in their transverse plane; this gives rise to the third contribution.

In the sub-relativistic regime, β E 1, u E β E , and the above expressions can be approximated as

Θ b • (b • ∇) u E Θ ⊥ ∇ • u E -Θ a E (1/c)∂ t u E . ( 6
)
The acceleration is thus one order higher in β E than the shear components because of the time derivative; consequently, its contribution can be neglected when β E 1.

In the present scheme, no particular difference arises between particles of different energies. Although the rhs of Eq. ( 4) scales as the square of the particle four-velocity in R / E , one power is absorbed by the fact that the derivative is taken with respect to proper time τ on the lhs, because dτ = dt/γ in the lab frame (here, γ denotes the particle Lorentz in the lab frame). The remaining power implies that energy gain/loss scales in proportion to the particle energy, as in the original Fermi model.

Figure 1 provides a contour map of the quantity Θ 2 ⊥ in the JHU-MHD simulation of incompressible turbulence, superimposed on a density plot of the magnetic field energy density. The distribution of Θ 2 is not shown here, but comparable to that of Θ 2 ⊥ . The contribution of the gravitational term ∝ a E • b is smaller than that of the other two; for the JHU-MHD turbulence simulation, u E 2 1/2 0.41 (units of c). Note that, for linear Alfvén waves, Θ ⊥ = Θ = 0 exactly, while for linear magnetosonic modes, Θ = 0 but Θ ⊥ ∼ β MS k δB/B = 0, where β MS denotes the phase velocity of magnetosonic modes (units of c), k the wave number and δB/B the relative mode amplitude.

In the context of quasilinear theory, particle acceleration scales with the power spectrum of magnetic fluctuations, hence with the turbulent magnetic energy density, while in the present nonresonant description, acceleration scales with the shear of the field line velocity flow. ⊥ tends to take large values in regions of strong gradients of the magnetic energy density, with some preference for regions with lower than average magnetic energy density. Such regions are subject to the magnetic pressure from the neighbouring regions of strong B; at the same time, small scale modes exert a larger (relative) influence on the magnetic field direction in zones of magnetic depressions, which can also induce strong gradients.

This correlation with the gradients of the fields is of particular importance in the context of modern turbulence theories, as it directly connects particle acceleration with turbulence intermittency, and with the coherent structures that characterize this intermittency.

B. Connection to plasma heating/cooling

The two terms Θ ⊥ and Θ appear in different combinations in Eq. ( 4), because of the nature of particle gyromotion, but they generalize the idea that a particle trapped in a flow undergoing compression or dilation gains or looses energy. Assume indeed, for the time being, that the particle four-velocity is isotropic at all times: the acceleration contribution then vanishes because u = 0, while both longitudinal and transverse shear contributions add up into the isotropic combination

1 c dγ dτ = - 1 3 u 2 ∂ α u E α , (7)
as in [START_REF] Lemoine | Generalized Fermi acceleration[END_REF] for isotropic scattering. In the subrelativistic limit,

∂ α u E α ∇ • β E .
In a fluid description, adiabatic heating/cooling through compression/dilation is expressed through a similar operator, although it is expressed in terms of the fluid velocity u p , not the velocity of field lines u E . The fluid approximation tacitly assumes that particles are isotropized in the fluid rest frame through some effective collision, while in the collisionless case, what governs the energization is the velocity of the field lines to which the particle is attached.

The quantities ∇ • u E and ∇ • u p differ from one another. To express this difference and relate the gradients of u E to the sources of heating and cooling, we assume here u E 1 to make contact with subrelativistic MHD; similar developments can be found in [START_REF] Le Roux | A kinetic transport theory for particle acceleration and transport in regions of multiple contracting and reconnecting inertial-scale flux ropes[END_REF][START_REF] Li | The roles of fluid compression and shear in electron energization during magnetic reconnection[END_REF]. Note first that

u E = u p -b (u p • b) , (8)
so that

∇•u E = ∇ ⊥ •u p + u p • B B 4 B •∇B 2 - 1 B 2 u p •(B • ∇) B , (9)
which indicates that the difference between the compression of u E and that of u p mainly relates to field line curvature.

From Poynting's theorem to first order in the electric field contribution, 1 2c

∂ ∂t B 2 = -E • (∇ × B) -∇ • (E × B) , (10)
we obtain

∇ • u E = - 1 2B 2 c d dt B 2 - 1 B 2 u E • (B • ∇) B , (11)
with d/dt = ∂ t + cu E • ∇. This equation relates the compression of the field lines to the change of magnetic energy density along the field line trajectory and to the field line tension. The last term on the rhs of the above equation, in particular, can be recast as Θ , given that

u E • B = 0. Therefore, Θ ⊥ = ∇ ⊥ • u E = - 1 2B 2 c d dt B 2 . ( 12
)
This relation can also be obtained directly from the advection equation. In the fully relativistic regime, its analog is provided in Appendix A,

1 2B 2 1 c ∂ ∂t B 2 = -Θ ⊥ , (13)
providing a direct connection between the tranverse compression of field lines and the temporal evolution of the magnetic energy density in the R / E frame. The terms that contribute to energization, up to the inertial one, are thus connected to the sinks to magnetic energy. In this regard, the acceleration process which is depicted here does not differ much from the physics of heating in a plasma of infinite conductivity and, to a first approximation, those particles that get accelerated to high energies can be regarded as those that have had the chance to encounter more zones of heating than the rest of the population. The shape of the spectrum at high-energy is then determined by sequences of outlier scattering events. Those will be captured in Sec. IV B through the use of large deviation theory.

C. Connection to the guiding center description Equation ( 4) has been obtained under the assumption that the particle experiences inertial forces in the R / E frame, on spatiotemporal scales larger than its gyration around the magnetic field line, with proper averaging over this gyromotion and in the frame of ideal MHD. For this reason, it captures to this order the contribution of guiding center theories (at their first order).

The detailed relationships connecting the terms of Eq. ( 4) to those of guiding center theories are provided in Appendix A. It is shown, in particular, that the term involving Θ represents the curvature drift, while the term Θ ⊥ captures the grad-B drift, i.e the mirror force responsible for transit-time damping acceleration, and the gravitational term is connected to the inertial drift. Equation (13) above also expresses the betatron nature of the grad-B drift.

These two schemes of approximation remain of a different nature, however, and differences should be seen at higher orders. Theories of guiding center motion are perturbative schemes based on expansions in r g /L, where L ∼ c represents the typical length scale on which the electromagnetic fields evolve. In the framework of Hamiltonian guiding center models, one defines at each order an adiabatic invariant conjugate to a gyrofrequency, its invariance expressing the fact that the short timescales of gyromotion have been averaged out of the dynamics. Such theories are difficult to expand, however, given their inherent sophistication, e.g. [START_REF] Boghosian | Covariant Lagrangian methods of relativistic plasma theory[END_REF][START_REF] Beklemishev | Covariant descriptions of the relativistic guiding-center dynamics[END_REF][START_REF] Cary | Hamiltonian theory of guiding-center motion[END_REF], and the equations of motions, written for quantities that are redefined order by order, rapidly take a complicated form.

By contrast, Eq. ( 3) is exact, as it simply expresses the equation of motion in a frame in which the electric field vanishes. The equation that we use in this paper, Eq. ( 4), is itself an approximation of that equation, to lowest order in r g / c : namely, the derivatives ∂ γ u E β and ∂ γ b β have been calculated at the position of the guiding center and the perpendicular velocity components u 2 and u 3 have been averaged out over this gyromotion. To go to higher orders, one would nevertheless start from Eq. (3) again and include different approximations for those terms.

In this sense, the present formalism may allow for a more phenomenological approach to the problem than guiding center theories, which are based on a specific hierarchy in r g / c . For instance, Sec. IV will account for particle diffusion along the field line, so that the particle position in Eq. ( 4) will be treated as a random variable, in the spirit of extensions to quasilinear theories [START_REF] Völk | Nonlinear perturbation theory for cosmic ray propagation in random magnetic fields[END_REF][START_REF] Matthaeus | Nonlinear collisionless perpendicular diffusion of charged particles[END_REF][START_REF] Yan | Cosmic-Ray propagation: non-linear diffusion parallel and perpendicular to mean magnetic field[END_REF][START_REF] Shalchi | Nonlinear Cosmic Ray Diffusion Theories[END_REF]. A similar improvement over the present description would be to account for perpendicular transport in space, through field-line wandering combined with particle jump from field line to field line.

One should also consider the influence of drifts beyond the classical E × B drift: such terms, spurred by large-scale modes, would modify the lowest order form for the particle four-velocity u b and therefore generate additional force terms, of higher order in the power of the gradients of the velocity field. As mentioned earlier, one may also consider introducing terms beyond ideal Ohm's law in the definition of the electric field, to capture physics at smaller length scales. Regarding smallscale modes, one extension that deserves investigation, is their contribution to the velocity of the field line, as alluded to earlier. This effect, of order O δB 2 l /B 2 , would also introduce higher powers of the velocity fluctuations.

D. Violation of the adiabatic invariant and pitch-angle scattering

The adiabatic invariant M ≡ u ⊥ 2 /(2B ) is conserved to the present order of approximation, as can be explicitly verified using Eq. ( 4). Such conservation represents an obvious barrier to stochastic acceleration to large momenta, because B can be regarded as a bounded quantity, just as the pitch angle cosine in R / E , µ . In this regard, higher order effects, such as those discussed above, are a necessity; see also Refs. [START_REF] Pecora | Ion diffusion and acceleration in plasma turbulence[END_REF][START_REF] Cerri | On stochastic heating and its phase-space signatures in lowβ kinetic turbulence[END_REF] for a related discussion. An important issue, therefore, is the hierarchy between the rate at which M is violated and that at which particles can gain momentum according to the model described by Eq. ( 4). We argue here that those two rates are likely of the same order of magnitude.

Note first that the inertial forces that control energization also contribute to scattering, as in any Fermi process. In particular, Eq. ( 4) implies

1 c dµ dτ = - γ 2 u 1 -µ 2 a E • b -γ µ 1 -µ 2 Θ - 1 2 u 1 -µ 2 b • ∇ ln B + 1 2 γ µ 1 -µ 2 Θ ⊥ . (14)
In a turbulence of small r.m.s. velocity fluctuation, the third term, which corresponds to the mirror force, dominates, because all others scale with u E . This suggests that particle scattering would be dominated in that regime by magnetic mirroring effects. In strong turbulence, such mirrors may differ from their linear magnetosonic mode counterparts, as discussed in Sec. II.

However, what matters to particle acceleration, here, is chaoticity, and in that respect, one must distinguish the violation of M , which at fixed u and B corresponds to a stochastic excursion of µ , from the regular, nonchaotic evolution of the pitch-angle cosine in a M -conserving environment.

Interactions that violate the conservation of M likely arise from small-scale perturbations and/or from regions of magnetic depressions where a magnetic fluctuation δB can strongly distort the mean field direction. At one extreme, large-scale modes (l r g ) renormalize the adiabatic invariant by contributions of order δM/M ∼ O [(r g /l) (δB /B )] [START_REF] Littlejohn | Hamiltonian formulation of guiding center motion[END_REF]. At the other extreme, small scale perturbations (l r g ) impart a deflection δu ⊥ /u ⊥ ∼ O [(l/r g ) (δB /B )] over their length scale l. This brings a change in M to the order δM/M = 2δu ⊥ /u ⊥ -δB /B .

At a fixed value of δB /B , the comparison of the above large-scale and small-scale influences suggest that the dominant contribution comes from scales comparable to the gyroradius. However, intermittency implies that large excursions in δB /B become more frequent at smaller length scales, which means that, overall, order unity violations of M are more likely to occur from small scale modes. Shock waves or reconnection layers provide examples of microscopic phenomena that can lead to substantial variation in δB .

Magnetic perturbations δB /B , when coarse-grained on spatial scales r g , also set the magnitude of the (M -conserving) forces that control pitch-angle scattering, see Eq. ( 14), hence acceleration, given that ∆u /u ∼ µ ∆µ /(1-µ 2 ) at fixed M and B . Consequently, one can reasonably expect that the rate of M -violating events is comparable -or even larger, if one takes into account the contributions of modes l < r g that can violate M but do not contribute to energization in the present modelto that of acceleration. Finally, acceleration itself takes place on timescales larger than the gyrotime t g because it results from the influence of modes on scales l > r g ; consequently, approximate gyromotion remains a valid approximation at all times. This validates the assumptions that underlie the main Eq. ( 4).

Another way to picture the energization process is thus the following. At fixed M , the phase space variables u and µ evolve in a banana-shaped island defined by u 2 (1 -µ 2 ) = 2M B , whose width is controlled by the excursions of B around its mean. Islands corresponding to two neighboring values of M are therefore not disconnected, and the junction of all islands for all values of M span the whole phase space. Hence, provided that the rate of evolution of M is comparable or larger to that of u , acceleration takes place at the rate fixed by Eq. (4).

IV. PARTICLE ENERGIZATION

The objective of this section is to quantify particle energization in the frame of the present model. It is split into three main parts. The first Sec. IV A calculates the one-and two-point statistics of energy gain, that is the mean energy gain ∆p and the diffusive term ∆p 2 .

Emphasis is placed here on the ambiguity that arises in the calculation of the latter quantity, given that the force terms that control energization display highly non-Gaussian features characteristic of intermittency. Section IV B then studies the acceleration process in a microscopic picture, which describes the momentum trajectory as a discrete random walk in a set of intermittent structures. This will clarify the role of intermittency and explain how it shapes the momentum distribution of accelerated particles. Finally, Sec. IV C provides some closing remarks.

A. Drift vs diffusion

Mean drift

In a purely kinematic description of Fermi acceleration, particles gain energy through head-on collisions with moving magnetized scattering centers, lose energy through tail-on collisions. On average, particles suffer more head-on than tail-on collisions, hence Fermi acceleration comes with a net (log-)momentum gain ∆ ln p /∆t (in the laboratory frame); the latter is comparable, in order of magnitude, to the stochastic diffusive term ∆ ln 2 p /∆t.

The present model departs from the Fermi scenario in one fundamental way: as the particle interacts with a region in which Θ = 0 or Θ ⊥ = 0, it does not gain or loose energy according to whether the interaction is head-on or tail-on, but according to whether Θ (or Θ ⊥ ) is negative or positive; see Eq. ( 4) 4 . Inside any such region, the acceleration process is regular, akin to a Fermi type-1 process, and it is rather through its encounter with many different sites of alternating polarity that the particle undergoes a random walk in momentum space. Consequently, the particle population gains (or loose) energy at a net rate depending on the sign of Θ and Θ ⊥ , where the average is computed over the spatial volume (assuming a homogeneous density distribution of particles).

As discussed earlier, those mean values relate to turbulent plasma heating (in the MHD approximation) through large scale motions. In the JHU-MHD turbulence simulation, Θ +0.03v A / c , while Θ ⊥ -0.1v A / c , and a E • b 0. Assuming a homogeneous and isotropic distribution of particles, this implies net heating at a rate

d ln p γ dτ = 1 3 Θ + Θ ⊥ 0.02 v A / c . (15)
Those values do not depend on the coarse-graining spatial scale ∆x that is adopted in calculating the spatial average, provided ∆x < c . On scales ∆x > c , all averages rapidly tend to zero, which is expected insofar as energy is continuously injected on scales c in that simulation. In a steady state, energy is thus spread over the cascade, down to microscopic scales where it is dissipated through nonideal processes. In the frame of the present model, the energy contained in sheared motions on scales l contributes to net heating for particles with gyroradius r g l.

The above heating rate indicates a rather small drift of the particle distribution on an eddy turnover timescale c /v A . In the context of particle acceleration, what truly matters, however, is the acceleration timescale, which will be found to be of the order of c c /v 2 A for that same simulation. Extrapolating this to other regimes suggest that, at small values of v A , heating may dominate over particle acceleration, while at larger values of v A (relativistic range), acceleration may dominate over the drift. A word of caution is needed to specify what is understood by v A here: in the JHU-MHD simulation, it corresponds to u 2 E 1/2 , meaning that it represents a rms velocity fluctuation, not a standard Alfvén velocity defined with respect to a coherent magnetic field. As Eq. ( 4) makes clear, indeed, what controls the heating (and the acceleration) is the magnitude of velocity fluctuations.

Those remarks illustrate another difference of the present scheme with respect to quasilinear theory, where the velocity that enters the diffusion coefficient is the phase velocity of the linear waves, defined relatively to the mean field.

Second-order moment

The diffusion coefficient can be derived as usual (in R / E ) from

D p p = lim ∆t → +∞ 1 2∆t dτ 1 dτ 2 dp dτ 1 dp dτ 2 . ( 16
)
with p the particle momentum in R / E . For reference, the equation of evolution of the momentum reads

1 c dp dt = -γ mc µ a E •b-p µ 2 Θ - 1 2 p 1 -µ 2 Θ ⊥ . (17)
The average in Eq. ( 16) involves the correlation function at two different space-time points of the three force terms in Eq. ( 4), including their cross-correlation functions. It also involves a limit ∆t → +∞, which, less formally, means that the time interval should be much larger than the coherence time of the random force that is exerted on the particle. This point will be discussed in greater detail further below.

The main contribution to the integral is rewritten using

dt 1 dt 2 Q 1 Q 2 = 2∆t dk R k S Q k , (18)
with the following notations: Q 1 (respectively Q 2 ) denotes one of the three random forces, that is the inertial, the parallel shear or the perpendicular compressive term, calculated at the particle position at time t 1 (resp. t 2 ). On the rhs, S Q k = Q k 2 represents the (one-dimensional in k-space) power spectrum of Q. We have neglected here the cross-correlations between those terms, for the sake of clarity, and because those crosscorrelations are smaller in magnitude (as checked in the JHU-MHD simulation). We have also assumed, as is customary, a random-phase approximation, which allows us to write the rhs of Eq. ( 18) as a one-dimensional integral over wave numbers, thus separating the contributions mode by mode. Our goal here is not to provide a detailed calculation of the diffusion coefficient, but to illustrate the ambiguity that arises in this calculation when non-Gaussian features arise, and to motivate the model developed in the next sub-section. Those assumptions are thus not restrictive. The resonance function R k captures as in Sec. II the duration of the interaction with the structure. We assume here ballistic trajectories, in which case R k ∼ 1/(k v). Finally, we will also assume for definiteness a Goldreich-Sridhar type spectrum, in which case k → k ⊥ in Eq. ( 18) and k → k

2/3 ⊥ k 1/3 min , with k min ∼ 1/ c .
Because Q involves a gradient, one generally expects

S Q k ∝ k 2 S u E k , with S u E k
∝ k -qu the (one-dimensional in k-space) power spectrum of fluctuations of u E . However, Θ ⊥ represents a perpendicular gradient, while Θ is a gradient measured parallel to the mean field, hence

S Θ ⊥ ∼ k ⊥ 2 S u E k and S Θ ∼ k 2 S u E k .
The above calculation of the diffusion coefficient would then result in the following scalings (for q u = 5/3):

D Θ p p ∼ p 2 δu E 2 k min c ln k k min , D Θ ⊥ p p ∼ p 2 δu E 2 k min c k k min 2/3 , ( 19
)
where δu E 2 provides the overall amplitude (squared) of u E velocity fluctuations. The momentum k represents the maximum wave number up to which the integral can be performed. This wave number is set by the coarsegraining procedure, which considers only scales above r g , i.e. k ∼ r -1 g , or by the maximum wave number of the cascade k max , whichever is smaller. While the parallel shear contribution involves all scales above r g , the perpendicular compressive one peaks on small scales.

The above calculation recovers the quasilinear prediction for transit-time damping in fast mode waves (described by the Θ ⊥ contribution in its linear limit) if one The ambiguity that arises in the above calculation is related to the scaling of S Q k with k, once Q departs from a Gaussian random field. This problem is illustrated by the statistics of the parallel shear and perpendicular compressive term extracted from the JHU-MHD turbulence simulation. Those are presented in Fig. 2, expressed as Q 2 1/2 ∆x , meaning the r.m.s. coarse-grained over scale ∆x. The correspondence with the spectra in k-space is obtained as usual, through Θ 2

⊥ ∆x ∼ kS Θ ⊥ k k∼∆x -1 at k ∼ k ⊥ . We thus expect Θ 2
⊥ ∆x ∝ ∆x -3+qu ∝ ∆x -4/3 for q u = 5/3, but Θ 2 ∆x ∝ ∆x -7/3+qu ∝ ∆x -2/3 for q u = 5/3.

The scaling Q 2 1/2 ∆x ∝ ∆x -2/3 (expected for Θ ⊥ and a E •b) is overplotted as a dotted line in Fig. 2. It appears to represent fairly well the statistics of the inertial term, but not those of Θ ⊥ , which flatten on small scales to an apparent Q 2 1/2 ∆x ∝ ∆x 0 . The parallel contribution Θ , as well as the mirror term (b • ∇ ln B ) 2 1/2 ∆x , also flattens out, instead of the expected Q 2 1/2 ∆x ∝ ∆x -1/3 . Such flat scalings would correspond to S Q k ∝ k -1 , hence to diffusion coefficients strongly dominated by the large scales,

D Θ p p ∼ p 2 δu E 2 k min c , D Θ ⊥ p p ∼ p 2 δu E 2 k min c . (20)
Out of the two results for the diffusion coefficients, Eqs. ( 19) and [START_REF] Jaekel | The Fokker-Planck coefficients of cosmic ray transport in random electromagnetic fields[END_REF], the latter appears more reliable, because based on explicit measurements in the JHU-MHD simulation. Neither is however able to capture the physics of particle acceleration, for about the same reason that a second-order moment cannot suffice to characterize a non-Gaussian distribution. This is nicely illustrated in Fig. 3, which plots the probability distribution function of those force terms for various smoothing scales, as extracted from the JHU-MHD turbulence simulation.

This figure reveals that the mean of the inertial contribution shifts to larger values as the coarse-graining scale decreases, in general agreement with the scaling observed in Fig. 2, and with the expectations for a Gaussian random field. By contrast, the distribution functions dN/d ln Q of the three other quantities, namely Θ , Θ ⊥ and b•∇ ln B , develop power law tails up to large values, with slopes close to (yet slightly larger than) -1. The mean of those three quantities hardly shifts as the coarsegraining scale decreases, but the slope of the power law becomes harder 5 . This explains, at least qualitatively, the tendency of Q 2 1/2 ∆x to fall below the predicted scaling in Fig. 2.

In this context, note that the variances shown in Fig. 2 have been evaluated as the median of the squared fluctuations, instead of their mean, in order to avoid an excessive source of noise associated to the finite statistics given the frequent, large excursions of the quantities. Those peculiar statistics are directly related to the scarcity of the structures, quite apparent in Fig. 1: the force terms can take large values in isolated regions of small spatial extent. This, in turn, directly impacts the acceleration process; in particular, it shapes the power law of accelerated particles, as discussed next.

B. Random walk in intermittent structures

Overview

The classical description of turbulence intermittency comes through the fractal β-model, which ascribes to daughter eddies a probability of being turned on as the parent cascades down [START_REF] Frisch | A simple dynamical model of intermittent fully developed turbulence[END_REF]. Said otherwise, active daughter eddies do not span the full volume of their parent eddy, but only a filling fraction, which is characterized by a fractal Hausdorff dimension D < 3. It is defined by D = ln N/ ln 2, if each parent eddy creates N daughter eddies out of 2 3 possible ones, when cascading from scale l to scale l/2. Correspondingly, at scale l, the filling fraction of active eddies is f l ∼ f c (l/ c ) (3-D) , with f c the filling fraction of active regions on the integral scale c . Note that Fig. 1 already suggests that f c takes values below unity. The value of D can take different values depending on the shape of the dissipative structures, for instance D = 1 for filaments, D = 2 for sheets. Its value, hence the filling fraction of active eddies, may depend on the general characteristics of the turbulence [START_REF] Padoan | Structure function scaling in compressible super-Alfvénic MHD turbulence[END_REF].

The mean free path to interaction with random, uncorrelated structures of extent l and filling fraction f l in three-dimensional space is l int ∼ f -1 l l or, for f l as calculated above,

l int ∼ f -1 c c (l/ c ) D-2
, which can take values much larger than l if D takes values below 3, and even larger than c if D 2. The probability of interaction per distance traveled in a fractal structure has been calculated and studied in detail in Ref. [START_REF] Isliker | Random walk through fractal environments[END_REF]. It leads to a slightly more complex form for the mean free path to interaction. Interestingly, it demonstrates that the probability of "no interaction" remains nonzero up to infinite distances if D < 2. In such a case, the random walk can be modeled as a continuous time random walk with a Lévy type distribution in waiting time. This situation has been considered in Ref. [START_REF] Lemoine | Power-law spectra from stochastic acceleration[END_REF], where it is shown that the momentum distribution then takes on a (near) power law shape at all times. We will ignore this possibility in the following.

The intermittent nature of acceleration sites introduces a new timescale t int = l int /v, whose magnitude (relative to other timescales) controls the phenomenology, if large enough. Here, t int sets the coherence time of the random force that acts on particles. Consequently, a diffusion coefficient becomes meaningful only on time intervals ∆t t int . On timescales ∆t < t int , different particles undergo vastly different histories, some interacting many times, others none. As discussed further below, this gives rise to the emergence of power law spectra. Such spectra have been observed recently in simulations of test-particle transport in synthetic fractal turbulent environments [START_REF] Sioulas | Stochastic turbulent acceleration in a fractal environment[END_REF].

On large timescales, ∆t t int , particles have had time to interact with many acceleration sites, hence one expects, in the spirit of the central limit theorem, to recover the predictions of a standard random walk in momentum space. However, this convergence will be found to be asymptotic, in the sense that the spectrum will preserve a power law tail at large momenta at all times. This will be made more explicit further below.

Random walk in intermittent turbulence

We determine the time-dependent distribution function in a microscopic picture, which describes a particle undergoing a random walk in a bath of isolated structures, of extent l and filling fraction f = f + + f - 1. Here, f + (respectively f -) represents the filling fraction of those structures giving rise to energy gain ∆p /p = +g (respectively energy loss ∆p /p = -g). In the frame of the present model, those structures have a negative (respectively positive) sign of Θ or Θ ⊥ . We drop the primes in this section, for clarity.

At each step of length l, the particle thus has a probability f -of loosing a fraction g of its momentum, a probability f + of gaining g, and a probability f ∅ = 1 -f --f + of crossing a void region, which corresponds to an inactive eddy in the framework of the β-model. The time-dependent solution to this random walk is captured thanks to the methods of large deviation theory [START_REF] Touchette | The large deviation approach to statistical mechanics[END_REF], which is ideally suited to model outlier statistics before central limit convergence.

with

S p ≡ f 2 ∅ ln 2 p + 4f -f + ln 2 pmax -ln 2 p . ( 22
)
This makes use of the following short-hand notations: p = p/p inj , where p inj denotes the momentum at which particles are injected in the turbulent bath; t = t/τ , where τ denotes the time between two steps, so that, correspondingly, t indicates the number of steps taken; ln pmax = g t, which represents the maximum logmomentum that the particle can achieve in t. We can assume τ l/v (v the particle velocity), which corresponds to ballistic trajectories.

Two combinations of the parameters play a particular role:

t int = τ / f + + f --(f + -f -) 2 ∼ l/(f + + f -)v,
the interaction time with active regions, and f + -f -, which controls the amount of net heating. The momentum distribution function displays the following properties:

1. The mean log-momentum ln p scales in proportion to f + -f -; see Fig. 4 for an illustration. Defining this mean log-momentum as the peak of p 3 f (p, t) ∼ dN/d ln p, we find ln p = g t (f + -f -). This result agrees with the discussion of Sec. IV A, where it is explained that the drift scales with the spatial average of positive minus negative contributions, that is f + -f -.

2. The distribution possesses a log-normal core with an extended power law tail. For the log-normal core, the variance can be expressed as

∆ ln 2 p = ln 2 p -ln p 2 = g 2 t f + + f --(f + -f -) 2 , (23)
Consequently, at any time, the particle distribution appears to follow a diffusive process characterized by the diffusion coefficient, introducing the effective interaction timescale

D pp = g 2 p 2 f + + f --(f + -f -) 2 τ = g 2 p 2 t int , (24)
t int = τ / f + + f --(f + -f -) 2 . Note that t int
approximately matches the time between two interactions with active eddies, l/(f + + f -)c, at small values of f + and f -. For t int ∝ p 0 and g ∝ p 0 , as discussed below, we recover the scaling D pp ∝ p 2 observed in numerical simulations.

3. This diffusion coefficient only represents an average value, measured over the particle population. It characterizes the general broadening ∆ ln 2 p of the core of the spectrum, but not its shape at large momenta, which takes on a power law shape. This general spectral shape -thermal log-normal core FIG. 5. Spectral index of the power law tail, as derived from Eq. ( 25) and plotted as a function of time in units of tint/g 2 = l/ g 2 (f+ + f-) v , for g spanning the range from 0.3 to 3 in steps of 0.3 (from bottom to top, or dark to light). This figure assumes f+ = f-. This spectral index is measured here at a momentum p = 10pinj.

with power law extension -matches that observed in recent kinetic numerical simulations [START_REF] Zhdankin | Kinetic turbulence in relativistic plasma: From thermal bath to nonthermal continuum[END_REF][START_REF] Zhdankin | System-size convergence of nonthermal particle acceleration in relativistic plasma turbulence[END_REF][START_REF] Zhdankin | Numerical investigation of kinetic turbulence in relativistic pair plasmas -I. Turbulence statistics[END_REF][START_REF] Comisso | Particle acceleration in relativistic plasma turbulence[END_REF][START_REF] Wong | First-principles demonstration of diffusive-advective particle acceleration in kinetic simulations of relativistic plasma turbulence[END_REF][START_REF] Comisso | The interplay of magnetically dominated turbulence and magnetic reconnection in producing nonthermal particles[END_REF][START_REF] Zhdankin | Electron and ion energization in relativistic plasma turbulence[END_REF].

4. The spectral index s, defined through p 2 f (p, t) ∝ dN/dp ∝ p -s , can be derived as

s = 1 - d d ln p ln dN d ln p = 1 - 1 g ln f ∅ ln p + S p 2f + (ln pmax -ln p) . ( 25
)
The power law is thus approximate, as the index runs logarithmically with ln p. This noteworthy result reveals that the extrapolation of the spectral index measured in numerical simulations -which are restricted in dynamic range -to a turbulent cascade extending over many orders of magnitude, is not straightforward. The power law tail extends here from the r.m.s. ∆ ln p 2 1/2 ∼ g (t/t int ) 1/2 to ln pmax = gt/τ ; its extent in log-space scales as (t/t int )

1/2 /f . 5. The power law shape hardens with time, just as it hardens with increasing energy gain per interaction. Those features are illustrated in Fig. 4. Those properties also agree with the measurements of power law spectra in kinetic simulations of relativistic turbulence [START_REF] Zhdankin | Kinetic turbulence in relativistic plasma: From thermal bath to nonthermal continuum[END_REF][START_REF] Zhdankin | System-size convergence of nonthermal particle acceleration in relativistic plasma turbulence[END_REF][START_REF] Zhdankin | Numerical investigation of kinetic turbulence in relativistic pair plasmas -I. Turbulence statistics[END_REF][START_REF] Comisso | Particle acceleration in relativistic plasma turbulence[END_REF][START_REF] Wong | First-principles demonstration of diffusive-advective particle acceleration in kinetic simulations of relativistic plasma turbulence[END_REF][START_REF] Comisso | The interplay of magnetically dominated turbulence and magnetic reconnection in producing nonthermal particles[END_REF][START_REF] Zhdankin | Electron and ion energization in relativistic plasma turbulence[END_REF], which show that the spectrum hardens with time, with increasing magnetization and with increasing fluctuation level δB/B. Note that the energy gain per scattering event is related to these latter two quantities, magnetization and turbulence amplitude.

6. The spectral index mainly depends on g and t int , for f + ∼ f -. In the power law range, indeed, S p f ∅ ln p and 2f + (ln pmax -ln p) ∼ 2f + ln pmax ∼ gt/t int . More specifically, at time t ∼ t int /g 2 , s takes values close to 3 for g ∼ 1, see Fig. 5. The quantity t int /g 2 represents the average acceleration timescale, hence this spectral index matches nicely the values s ∼ 3 observed in kinetic simulations of relativistic turbulence if t int ∼ c /c.

Finally, the spectrum converges asymptotically to the log-normal shape, as announced earlier. This is manifest from Eqs. (B2) and (B6) in the limit N → +∞ (with N = t the number of steps taken). This property of large deviation functions derives from the central limit theorem. On asymptotic timescales, diffusion then takes place as in homogeneous turbulence because the particles can probe all structures many times over, although the effective interaction timescale is now set by t int , and the acceleration timescale by t int /g 2 . To be more precise, the distribution function preserves at all times a nonthermal power law tail, which is pushed to larger and larger values of the momentum, because the Gaussian core extends farther in ln p, given that ∆ ln p 2 1/2 ∼ g (t/t int )

1/2 . The convergence to the standard log-normal shape is thus asymptotic.

In kinetic numerical simulations, the power law slows down its evolution in time once the maximum energy has reached the confinement limit, viz. r g (p max ) ∼ c , which takes place on a timescale t ∼ t int /g 2 ln [ c /r g (p inj )] [START_REF] Zhdankin | System-size convergence of nonthermal particle acceleration in relativistic plasma turbulence[END_REF]. Acceleration slows down at larger momenta, because the particle decouples from the turbulence (its scattering timescale increases fast with energy) [START_REF] Lemoine | Generalized Fermi acceleration[END_REF]. This bottleneck results in a piling-up of the momenta close to this limit. It is not captured by the present random walk, which considers a momentumindependent energy gain per scattering event. However, the above timescale is significantly longer than the acceleration timescale t int /g 2 , all the more so if the dynamic range of the cascade c /r g (p inj ) is large.

Estimates for g and tint

The above analytical model assumes a given value for g, t int and the scale l. Those can be estimated from the values that maximize the acceleration rate. Here, we extract them from the JHU-MHD simulation, using two different methods. We first consider a given coarsegraining scale l r g , then compare the contributions of various scales l.

In the presence of a broad distribution of values of Q (with Q symbolizing the force contributions as before), the maximum contribution to the diffusion rate at a given l maximizes Q 2 dN/d ln |Q|, see Eq. ( 24) above, because g ∝ Q. In this first exercise, we thus extract this optimal value of |Q| for each of the three force terms. We then derive the optimal relative energy gain as g = |Q| l, see Eq. ( 17), and extract the corresponding value of t int as l/(f l v), where f l (> |Q|) represents the filling fraction of FIG. 6. Estimates of the effective acceleration timescale tint/g 2 (upper curves, in units of c c/v 2 A , assuming v c) and of the average energy gain per interaction g (lower curves, in units of vA/c) in the JHU-MHD simulation, for different coarse-graining scales. The contributions of the three force terms are colored as follows: red for Θ , blue for Θ ⊥ and magenta for a E • b. Solid lines are estimated using the first method discussed in the text, dashed lines using the second method. The gray areas delimit the dominant contributions (meaning larger g, smaller tint/g 2 ), which are associated to the Θ term.

regions giving a contribution larger than |Q| (thus, of either sign),

f l (> |Q|) = 1 N +∞ ln |Q| d ln |Q | dN d ln |Q | l . (26)
The results are shown in Fig. [START_REF] Chen | Recent progress in astrophysical plasma turbulence from solar wind observations[END_REF], which suggests that g ∼ 0.4v A /c and t int /g 2 ∼ O 20 c 2 c /v 2 A v are reasonable estimates for those parameters. The measured values do not depend strongly on l as long as l c , suggesting that all scales contribute about equally to the acceleration process. Finally, this exercise also suggests that, at least for the JHU-MHD simulation, the parallel shear contribution Θ seemingly provides the dominant contribution to energy gain.

Our second method for determining those parameters takes into account the possibility of particle trapping inside structures, which happens at values of Q sufficiently large to impart pitch-angle deflections by order of unity. Once inside a structure, the particle gains (or loses) energy in a regular way, as explained in Sec. IV A. Consequently, trapping may give rise to momentum gains in excess of Ql for some particles. Such effects deserve a careful investigation, which lies beyond the scope of the present paper. Here, we seek a first-order estimate of t int from the structures that are strong enough to give rise to particle trapping, and for such structures, we derive a first-order estimate of g as the minimum (most conservative) value, set by g = Ql. The actual gain g is determined by that acquired over the lifetime of the structure, or by the time the particle takes to escape from the structure, whichever happens first. As for any such Fermi process, g is likely distributed as a power law whose index is controlled by the competition between the rate of energy gain Q vs the rate of escape ∼ 1/l. At values Ql below unity, this slope is large enough that the average value of g is set by its minimum value, g ∼ Ql.

Physically, this describes a situation in which the particle suffers a pitch-angle deflection of the order of unity through its crossing of the structure, which enhances its energy gain by a factor of unity over the ballistic estimate, but generally not more.

To determine t int , we use the distribution of the mirror force term Q µ ≡ b • ∇ ln B that dominates pitch-angle scattering and search for the filling fraction f l (> |Q µ |) such that |Q µ | = 1/l (guaranteeing pitch-angle deflection by order unity). We then set t int = l/f l (> |Q µ |) as before. Finally, we determine, for each of the three force terms Q contributing to particle energization, the energy gain g as Ql where Q is the value associated to the above value of the filling fraction. The results are overplotted in Fig. 6. They happen to be rather consistent with those found through the previous method, even though the two differ in their spirit. It so happens, however, that the value of Q µ , as determined through the previous method, is already close to unity, meaning that the structures that maximize the product Q 2 dN/d ln |Q| also tend to trap the particles.

More in-depth studies appear needed to establish a rigorous connection between the properties of the intermittency and the statistics that govern the acceleration rate as a function of scale. Nevertheless, the consistent estimates that emerge from the above exercises suggest that t int c /v, g v A /c provide reasonable estimates of the parameters governing acceleration in the JHU-MHD simulation. In turn, this suggests an average diffusion coefficient D pp ∝ p 2 (v A /c) 2 v/ c , as announced earlier, just as it predicts very steep spectra for particle acceleration in this JHU-MHD simulation configuration at a time of several t int c 2 /v 2 A , see Fig. 5.

C. Further remarks and consequences

The model developed in the previous Sec. IV B, which describes the momentum gain through a random walk in intermittent structures shares some interesting features with the observations made in dedicated numerical simulations: it agrees on the general scaling of the diffusion coefficient; it predicts a spectral shape similar to that observed, that is a log-normal core slowly drifting in time to larger momenta, with a power law tail; notably so, the spectral index agrees with that measured in those simulations. It would certainly benefit from further detailed analysis aimed at establishing a closer link between the spectral shape and the properties of the turbulence intermittency, as one may notably expect those properties to depend on various external parameters, such as the physics of stirring, the overall magnetization, the amplitude of the turbulent component, the plasma beta pa-rameter etc, e.g. [START_REF] Makwana | Properties of magnetohydrodynamic modes in compressively driven plasma turbulence[END_REF]. Such studies can involve analytical models, or numerical simulations of particle transport in synthetic intermittent settings, such as undertaken in Refs. [START_REF] Sioulas | Stochastic turbulent acceleration in a fractal environment[END_REF][START_REF] Vlahos | Particle acceleration in an evolving network of unstable current sheets[END_REF][START_REF] Turkmani | Particle acceleration in stochastic current sheets in stressed coronal active regions[END_REF][START_REF] Pucci | Energetic particle transport in the presence of magnetic turbulence: Influence of spectral extension and intermittency[END_REF][START_REF] Isliker | Fractional transport in strongly turbulent plasmas[END_REF][START_REF] Sioulas | Superdiffusive stochastic Fermi acceleration in space and energy[END_REF]. A more detailed analysis of the statistics of interactions between particles and structures in full-blown kinetic (or MHD) simulations of magnetized turbulence would also bring direct constraints and input for the present picture.

We also recall an important prediction of the above model: the (logarithmic) running of the spectral index of the power law tail with momentum. Such a running is not apparent in kinetic numerical simulations, and it should not be apparent given the restricted dynamic range that even the most massive simulations can probe. However, this running clearly deserves further investigation, as it may have a direct impact on phenomenology, for turbulent cascades extending over many orders of magnitude.

To conclude, it proves instructive to recast the above model in the frame of usual quasilinear or Fermi-like treatments of stochastic acceleration and mention a few noteworthy consequences of the impact of intermittency. In contrast to the present microphysical description of particle transport, most treatments and phenomenological applications of stochastic acceleration make use of the Fokker Planck equation to determine the particle distribution function f (p, t). This equation,

∂f ∂t = 1 p 2 ∂ ∂p p 2 D pp ∂f ∂p , (27)
here written without a drift (other than the noise-induced one) is entirely characterized by the momentum diffusion coefficient D pp ≡ ∆p 2 /(2∆t). This equation derives from the original Fermi model provided the interaction time between two scattering events does not depend on momentum [START_REF] Duffy | Cosmic ray transport and acceleration[END_REF], corresponding to the scaling D pp ∝ p 2 observed in numerical simulations [START_REF] Lynn | Acceleration of relativistic electrons by magnetohydrodynamic turbulence: Implications for nonthermal emission from black hole accretion disks[END_REF][START_REF] Zhdankin | Kinetic turbulence in relativistic plasma: From thermal bath to nonthermal continuum[END_REF][START_REF] Zhdankin | System-size convergence of nonthermal particle acceleration in relativistic plasma turbulence[END_REF][START_REF] Zhdankin | Numerical investigation of kinetic turbulence in relativistic pair plasmas -I. Turbulence statistics[END_REF][START_REF] Comisso | Particle acceleration in relativistic plasma turbulence[END_REF][START_REF] Wong | First-principles demonstration of diffusive-advective particle acceleration in kinetic simulations of relativistic plasma turbulence[END_REF][START_REF] Comisso | The interplay of magnetically dominated turbulence and magnetic reconnection in producing nonthermal particles[END_REF][START_REF] Zhdankin | Electron and ion energization in relativistic plasma turbulence[END_REF].

The same Fokker-Planck equation also derives from the quasilinear transport equation, when rewritten in the diffusion approximation, after averaging over spatial variables [START_REF] Skilling | Cosmic ray streaming -I. Effect of Alfvén waves on particles[END_REF][START_REF] Schlickeiser | Cosmic-Ray transport and acceleration. I. Derivation of the kinetic equation and application to cosmic rays in static cold media[END_REF]. For D pp ∝ p 2 , the solution for monoenergetic injection at momentum p inj takes a Gaussian shape in ln p,

p 3 f = (4πν diff t) -1/2 exp - [ln (p/p inj ) -3ν diff t] 2 4ν diff t (28
) with ν diff = D pp /p 2 . The mean of this distribution increases as fast as the variance, indicating substantial drift in momentum, in contrast with the present description; see Sec. IV A.

Formally, the Fokker-Planck equation represents a truncation to the second order of the Kramers-Moyal expansion of the more fundamental equation that governs the evolution of f . This truncation is exact if the corresponding microscopic dynamics can be described by a Langevin equation with Gaussian noise of infinitesimally small coherence time [START_REF] Risken | The Fokker-Planck Equation. Methods of Solution and Applications[END_REF]. It is thus not surprising that this solution fails to reproduce a power law extension in the present context because the intermittency effectively enlarges the coherence time of the random force to the macroscopic timescale t int . The above Gaussian shape is rather recovered as the asymptotic limit of the random walk in intermittent turbulence at t/t int → +∞, as discussed earlier. Several models have been proposed to generalize the Fokker-Planck equation to incorporate higherorder terms, e.g. effective resummations [START_REF] Bykov | Electron acceleration and efficiency in nonthermal gamma-ray sources[END_REF][START_REF] Bykov | Theory of charge particle acceleration by a shock wave ensemble in a turbulent medium[END_REF], fractional transport equations [START_REF] Isliker | Fractional transport in strongly turbulent plasmas[END_REF], or more pragmatically, a Fokker-Planck equation with a momentum-dependent advection coefficient, tuned in such a way as to reproduce the results of kinetic simulations [START_REF] Wong | First-principles demonstration of diffusive-advective particle acceleration in kinetic simulations of relativistic plasma turbulence[END_REF].

The role played by intermittency may have other, direct implications for phenomenology. For one, strong intermittency is associated to inhomogeneity, which directly implies nonuniformity in space and time of the radiation emitted by the accelerated particles [START_REF] Nättilä | Radiative turbulent flares in magnetically-dominated plasmas[END_REF][START_REF] Zhdankin | Kinetic turbulence in shining pair plasma: Intermittent beaming and thermalization by radiative cooling[END_REF][START_REF] Zhdankin | Production and persistence of extreme two-temperature plasmas in radiative relativistic turbulence[END_REF]. It can also impact the spatial transport of particles, just as it impacts their transport in momentum space. In particular, if the interaction length between two scattering events is ∼ c , as suggested by the analysis of the previous section, then the fraction of unscattered particles must remain of the order of unity on length scales c . This implies that significant anisotropies can be maintained on long length scales, with obvious consequences for cosmic-ray physics.

V. SUMMARY AND CONCLUSIONS

This paper has discussed the physics of particle acceleration in strong MHD turbulence characterized by an ensemble of structures, rather than by the linear superposition of random waves of quasilinear theory. To this effect, it has made use of a novel formalism, applicable from the subrelativistic to relativistic limits, that describes the history of the particle in a mixed phase space, with spatial variables expressed in the lab-frame, and momenta in the R / E frame that moves at velocity β E = E ×B/B 2 (in units of c; its four-velocity is written u E ), where the motional electric field vanishes. All energy gains or losses can then be expressed as the influence of inertial forces that derive from the noninertial nature of that frame R / E . These energy gains/losses are thus directly related to the gradient of the velocity flow of the magnetic field lines, in particular their acceleration, their compression and their shear. This allows us to connect in a direct manner the physics of particle acceleration to the description of the turbulent bath in terms of velocity structures. To quantify the terms of energy gain/loss, we have assumed that the particle undergoes local gyromotion around the perturbed field lines and retained the lowest order contribution in r g / c (ratio of particle gyroradius to coherence length scale). The main results can be summarized as follows:

1. Three terms contribute to particle acceleration: the effective gravity projected along the field line (a E • b), associated to the acceleration field a E of u E , the shear of u E along the field line (Θ ), and its compression in the transverse plane (Θ ⊥ ); the last two dominate over the first in subrelativistic turbulence. In ideal MHD, and to the present order of approximation, those forces capture the classical drifts of guiding center models. The present formalism is amenable to extensions that can capture terms of higher order in r g / c or nonideal contributions to Ohm's law.

2. The present model is closer in spirit to the original Fermi scenario of stochastic particle acceleration than to quasilinear theory (QLT). It differs from this original Fermi scenario in one crucial aspect: in each region of nonzero gradient Θ or Θ ⊥ , the particle gains (or loose) energy in a regular (nonstochastic) manner that depends on the sign of the gradient, not whether the collision with that region is head-on or tail-on. The acceleration process becomes diffusive once the particle has encountered several uncorrelated energizing sites.

3. On average, the particles suffer net energy gain or loss, whose amount is determined by the volume average of those gradients (weighted by the particle distribution). This mean energy gain can be directly related to the net heating of the plasma through large-scale motions. This also differs from the usual quasilinear treatment, which implicitly assumes an infinite reservoir of energy for particle energization, and which predicts similar amounts of net and stochastic energy gains.

4. As demonstrated by measurements carried out in the JHU-MHD 1024 3 numerical simulation of incompressible MHD, the acceleration sites are localized in isolated regions that occupy a small filling fraction of the overall turbulent volume, unlike QLT, which assumes a homogeneous space-filling background of waves. Furthermore, those acceleration sites are preferentially found in regions of strong gradients of the magnetic energy density, whereas QLT predicts an acceleration rate that scales with the turbulent magnetic energy density.

5. The circumscription of accelerating sites to sparse regions is a hallmark of intermittency. Correspondingly, the statistics of the velocity gradients exhibit highly non-Gaussian features characterized by power law tails leading to large excursions, all the more so as the coarse-graining scale decreases.

6. Intermittency plays a key role in shaping the distribution function of the accelerated particles. The spectrum has been derived here in the frame of a microscopic model describing the transport of particles in a bath of intermittent acceleration sites.

This spectrum exhibits a Gaussian core that broadens in time, as in a diffusive process, together with a (near-)power law at large momenta, which hardens with time. Overall, the spectral shape will thus differ from that evaluated assuming a homogeneous turbulence, which opens interesting prospects for a rich phenomenology.

7. In such a situation, the average diffusion coefficient characterizes the broadening of the Gaussian core, but not the power law tail. The reason why the standard Fokker-Planck treatment fails here is because the intermittency introduces a macroscopic timescale t int , defined as the typical interaction time with an accelerating site. In the JHU-MHD simulation, t int ∼ c /v (v particle velocity). The corresponding diffusion coefficient D pp ∝ p 2 δu 2 E v/ c agrees qualitatively well with that measured in test particle simulations of subrelativistic, weakly sub-Alfvénic MHD turbulence [START_REF] Lynn | Acceleration of relativistic electrons by magnetohydrodynamic turbulence: Implications for nonthermal emission from black hole accretion disks[END_REF], as well as that measured in fully kinetic simulations of relativistic turbulence [START_REF] Zhdankin | Kinetic turbulence in relativistic plasma: From thermal bath to nonthermal continuum[END_REF][START_REF] Zhdankin | System-size convergence of nonthermal particle acceleration in relativistic plasma turbulence[END_REF][START_REF] Zhdankin | Numerical investigation of kinetic turbulence in relativistic pair plasmas -I. Turbulence statistics[END_REF][START_REF] Comisso | Particle acceleration in relativistic plasma turbulence[END_REF][START_REF] Wong | First-principles demonstration of diffusive-advective particle acceleration in kinetic simulations of relativistic plasma turbulence[END_REF][START_REF] Comisso | The interplay of magnetically dominated turbulence and magnetic reconnection in producing nonthermal particles[END_REF][START_REF] Zhdankin | Electron and ion energization in relativistic plasma turbulence[END_REF].

8. The particles that populate the nonthermal power law are those that have been able to interact with sufficiently many sites in a given time. At time t ∼ t int /g 2 , with g the characteristic relative energy gain per interaction with an accelerating site, the spectral index s (dN/dp ∝ p -s ) is predicted to be of the order of 3 for g ∼ O(1). This agrees rather well with the aforementioned kinetic numerical simulations of relativistic turbulence. Smaller values of g (as expected in subrelativistic turbulence) lead to softer spectra.

As discussed in Sec. III, the energization in regions of large shear of the velocity field of magnetic field lines is closely related to the sources and sinks of magnetic energy in the MHD description, and more generally with intermittent structures characterized by large gradients. This should be put in perspective with numerous observations and simulations that have revealed that plasma heating in turbulent environments is indeed associated with intermittent structures [START_REF] Kiyani | Global scale-invariant dissipation in collisionless plasma turbulence[END_REF][START_REF] Osman | Evidence for Inhomogeneous Heating in the Solar Wind[END_REF][START_REF] Osman | Intermittency and local heating in the solar wind[END_REF][START_REF] Bruno | Intermittency in solar wind turbulence from fluid to kinetic scales[END_REF][START_REF] Wan | Intermittent dissipation at kinetic scales in collisionless plasma turbulence[END_REF][START_REF] Karimabadi | Coherent structures, intermittent turbulence, and dissipation in high-temperature plasmas[END_REF][START_REF] Osman | Magnetic reconnection and intermittent turbulence in the solar wind[END_REF], albeit generally interpreted as the consequence of reconnection.

Those two pictures enjoy intriguing connections. For one, reconnecting currents sheets are associated with violent bulk plasma motions that can be described in an MHD approximation over most of the volume, outside the diffusion regions. In such environments, particle energization is dominated by Fermi-type interactions, at least in the absence of a strong guide field, and all the more so at high energies, e.g. [START_REF] Dahlin | Parallel electric fields are inefficient drivers of energetic electrons in magnetic reconnection[END_REF][START_REF] Guo | Recent progress on particle acceleration and reconnection physics during magnetic reconnection in the magnetically-dominated relativistic regime[END_REF]. More precisely, the perpendicular electric field that governs those Fermi-type processes can be related to the shearing and compressive motions of the flow, as in the present formalism [START_REF] Li | The roles of fluid compression and shear in electron energization during magnetic reconnection[END_REF]. Zooming out to larger scales, particles could also be accelerated through repeated encounters with the exhausts of many, randomly oriented reconnection regions [START_REF] Hoshino | Stochastic particle acceleration in multiple magnetic islands during reconnection[END_REF],

or more generally with the nonuniform stirring motions around reconnection regions, if reconnection takes place in a turbulent fashion on all scales [START_REF] Lazarian | Reconnection in a weakly stochastic field[END_REF]. Therefore, what truly governs plasma heating in the vicinity of intermittent structures deserves closer scrutiny, from observational, theoretical and numerical viewpoints.

One merit of the present model is that it can be directly tested against observations or numerical simulations by comparing the observed amount of particle energization with the expected contributions for the inertial, the parallel shear and the transverse compression terms. To calculate those expected contributions, it suffices to measure the statistics of the gradients of u E and use Eq. ( 4). In that regard, Fig. 1 provides a direct tomography of the acceleration regions in the JHU-MHD simulation.

This contrasts somewhat with the quasilinear view of energization through wave particle resonant interactions, which are difficult to observe; see however Ref. [START_REF] Klein | Diagnosing collisionless energy transfer using field-particle correlations: Gyrokinetic turbulence[END_REF] for recent progress in that direction. For that reason, one commonly tests models of particle acceleration (or transport) through the comparison of the predicted diffusion coefficient with that measured in simulations. However, this does not suffice, because various theoretical models can lead to the same overall scaling with particle momentum, and furthermore, the diffusion coefficient itself does not fully characterize the acceleration process, as discussed above.

The role of turbulence intermittency regarding particle transport, a novel aspect of the present model, deserves further investigation. One of its most obvious signatures to look for, in observations or numerical simulations, is inhomogeneity in the rates of particle diffusion; such inhomogeneity has actually been observed in Refs. [START_REF] Trotta | Fast acceleration of transrelativistic electrons in astrophysical turbulence[END_REF][START_REF] Nättilä | Radiative turbulent flares in magnetically-dominated plasmas[END_REF]. Intermittency is likely to entail important consequences for phenomenological applications in space plasma physics. For instance, it may lead to significant anisotropies in phase space on scales comparable to the coherence length, for the same reason that it affects the momentum distribution function. One should thus aim at establishing a closer relationship between the properties of the turbulence (nature of the forcing, magnetization, amplitude etc.) and the phenomenology of particle transport and acceleration. This may involve test-particle propagation in synthetic turbulence as well as full numerical simulations of turbulence with advanced diagnostics. This choice of tetrad tacitly implies a rotation of the three spatial axes with respect to the lab frame. In particular, the magnetic four-vector in R / E is written B a = e a µ B µ = {0, B , 0, 0}. In the R / E frame, the four-velocity of the particle is written u a = γ , u , u 2 , u 3 . Noting that the Lorentz force in Eq. (A2) contributes neither to the evolution of γ nor to that of u , we obtain

1 c dγ dτ = -γ u b β u E γ ∂ γ u E β -u 2 b β b γ ∂ γ u E β -u 2 2
h E αβ = η αβ + u E α u E β , σ E αβ = 1 2 h E αµ h E βν ∂ ν u E µ + ∂ µ u E ν - 1 3 Θ E h E αβ , ω E αβ = 1 2 h E αµ h E βν ∂ ν u E µ -∂ µ u E ν , Θ E = ∂ α u E α , a E αβ = -u E β u E µ ∂ µ u E α . ( A8
)
The acceleration four-vector is a E α = u E β a E αβ = u E µ ∂ µ u E α . Given that the magnetic field vector defines a preferred direction, we split Θ E = Θ + Θ ⊥ with

Θ = b α b β ∂ α u E β , Θ ⊥ = η αβ -b α b β ∂ α u E β , (A9)
which characterize the shear along and the compression perpendicular to the magnetic field line.

This eventually leads to the main equation,

1 c dγ dτ = -γ u a E • b -u 2 Θ - 1 2 u ⊥ 2 Θ ⊥ , 1 c du dτ = -γ 2 a E • b -γ u Θ - 1 2 u ⊥ 2 b • ∇ ln B .
(A10)

The first two terms on the rhs of each equation derive directly from the decomposition of ∂ γ u E β and ∂ γ b β . The third are obtained as averages over a gyroperiod, assuming that on such short timescales, the gradients do not change significantly. Then, u 2 = u 3 = u 2 u 3 = 0 but u 2

h ⊥ βγ ∂ γ b β = 1 B h ⊥ βγ ∂ γ B β
The equations of motion can be derived using the Poisson brackets, Combining both equations of Eq. (A18) with the fourvelocity normalization u a u a + 2M B = -1 leads to dM/dτ = 0, expressing conservation of the adiabatic moment.

x α , u β = 1 m -u E α u E β + b α b β , (
We can establish a connection between the two formalisms and rewrite Eq. (A18) in the form of Eq. (A10), as follows. First note that the inertial contribution, viz. the first term on the rhs of Eq. (A18) breaks down into the first two terms on the rhs of Eq. (A10), noting that, in the guiding center approximation, 1 c

d dτ ≡ γ 1 c ∂ ∂t + u b • ∇ + γ β E • ∇ , ≡ γ u E α ∂ α + u b α ∂ α , (A19)
because

u E α ∂ α = γ E c -1 ∂ t + u E • ∇ = c -1 ∂ t , b α ∂ α = b • ∇ = b • ∇ . (A20)
Those two terms on the rhs of Eq. (A19) lead to the the effective gravity force projected along the field line direction (a E α b α ) and the shear acceleration or compression along the field line direction (Θ ). The latter characterizes the curvature drift part of the energy change. This is because the curvature drift velocity scales as

v d, curv ∝ B × (b • ∇b), so that E • v d, curv ∝ (E × B) • (∇b) ∝ b • (b • ∇u E ), given that b • u E = 0.
The third terms on the rhs of Eqs. (A10) and (A18) match one another in the expression for u ; they describe the grad-B force. In the equation for γ , this grad-B term can also be rewritten in terms of derivatives of the velocity field through the projection of the advection equation on b β ,

b β ∂ α u E α B β -u E β B α = 0 ⇒ u E α ∂ α B = -B η αβ -b α b β ∂ α u E β = -Θ ⊥ B . (A21)
Hence, the grad-B term appears here under a variety of forms: it can be depicted as a betatron process, as in Eq. (A22), as a mirror force, or as the result of compression in the plane transverse to the field lines. Finally, one can check that these equations match the guiding center equations of motion of Vandervoort and Northrop [START_REF] Vandervoort | The relativistic motion of a charged particle in an inhomogeneous electromagnetic field[END_REF][START_REF] Northrop | The guiding center approximation to charged particle motion[END_REF],

d γ dτ = q mc 2 E • v d + M ∂ ∂t B , (A22)
where γ = γ E γ corresponds to the Lorentz force of the particle in the lab-frame, averaged over a gyroperiod, with drift velocity,

v d = mc 2 q u B 2 B × db dτ + γ B 2 B × dβ E dτ + M c B 2 B × ∇B + M B 2 E ∂B ∂t .
(A23)

Rewriting Eq. (A22) above for γ instead of γ , it can be shown that the inertial drift contribution reproduces the first term of Eq. (A18), while the polarization drift vanishes in the proper frame, and the grad-B drift terms can be reorganized as the second term of Eq. (A18). Similar operations can be conducted to recover the evolution of u .

This gives

I(ŷ) = ŷ ln   f ∅ ŷ + f 2 ∅ ŷ2 + 4f -f + (1 -ŷ2 ) 2f + (1 -ŷ)   -ln   f ∅ + f 2 ∅ ŷ2 + 4f -f + (1 -ŷ2 ) (1 -ŷ2 )   . (B5)
The mean value for ŷ is ŷm ≡ ŷ = f + -f -, as is obvious from Eq. (B1). At ŷm , we have

I(ŷ m ) = 0 , I (ŷ m ) = 0 , I (ŷ m ) = 1 f -+ f + -(f + -f -) 2 .
(B6)

Consequently, for ŷ ∼ ŷm , corresponding to ln(p/p 0 ) ∼ N g(f + -f -), the distribution is nearly Gaussian,

φ(ŷ) ∼ N →+∞ exp - N 2 I (ŷ m ) (ŷ -ŷm ) 2 , (B7)
but, on intermediate timescales N < +∞, it exhibits a non-Gaussian tail at large values of ŷ.

FIG. 1 .

 1 FIG. 1. Density plot of the magnetic energy density (black and white, in units of the average energy density) in a twodimensional slice of the JHU-MHD simulation cube, with overlayed contours of Θ 2 ⊥ at (from thin to thick, in blue, solid lines): 0.5, 2., 8., 32. in units of (vA/ c) 2 , and similarly for negative values in yellow. Top panel: 1024 2 resolution, extending over about 10 c × 10 c; bottom panel: zoom over a region of 1.3 c × 1.3 c.

  Figure 1 thus provides a vivid illustration of the difference between these two pictures. It shows in particular that Θ 2

FIG. 2 .

 2 FIG. 2. Scaling of the r.m.s. of the terms contributing to particle energization and scattering, as indicated: inertial a E • b, parallel shear Θ and perpendicular compression Θ ⊥ (all three in units of vA/ c), and for scattering, the mirror term b • ∇ ln B in units of 1/ c, as extracted from the JHU-MHD turbulence simulation and plotted as a function of smoothing scale ∆x (units of c). The dotted line shows the scaling ∝ ∆x -2/3 ; see the text for details.

FIG. 3 .

 3 FIG. 3. Normalized distribution of the absolute values of the main contributions to the energization and scattering, for different smoothing scales: ∆x/ c = (0.01, 0.09, 0.18, 0.26, 0.43, 1.1), ordered from blue to red, or equivalently right to left. From top to bottom: |Θ |; |Θ ⊥ |; |a E • b|; |b • ∇ ln B |. The first three are expressed in units vA/ c, while the last is expressed in units 1/ c. Each subpanel plots dN/d ln Q, where Q represents one of the above quantities. The departure from Gaussian statistics, in particular the existence of power law tails is a clear signature of intermittency. The dotted line represents a power law of index -1.

FIG. 4 .

 4 FIG. 4. Spectra dN/d ln p from the stochastic acceleration of particles in intermittent turbulence, as obtained from Eq. (21), see also the discussion in Appendix B. Top panel: spectra plotted at different (normalized) times t = 12., 25., 50., 100. (ordered from red to blue or left to right), for g = 1 and f+ = f-= 0.01. Note that t = t/τ characterizes the number of steps taken, hence the average number of structures encountered is (f+ + f-) t, which takes values between 0.25 and 2 here. Middle panel: spectra plotted at time t = 50., for different values of g = 0.3, 0.6, 1.2, 2.4 (ordered from red to blue or left to right), and f+ = f-= 0.01. Bottom panel: spectra plotted at time t = 50., for different values of f+ = 0.02, 0.03, 0.045, 0.68 (ordered from red to blue or left to right), g = 1 and f-= 0.01.

B

  µ = {0, B } in the lab-frame. The dual field strength tensor is writtenF µν = 1 2 µναβ F αβ , with 0123 = -1.The last two spacelike vectors e 2 µ and e 3 µ span the plane orthogonal to u E µ and b µ . Their exact definition is not important in what follows. They can be defined through linear combinations of the eigenvectors of F a b with eigenvalues +iB and -iB .

e 2 β e 2 γ

 2 ∂ γ u E β -u 3 2 e 3 β e 3 γ ∂ γ u E β , β u E γ ∂ γ b β + γ u u E β b γ ∂ γ b β + u 2 2 e 2 β e 2 γ ∂ γ b β + u 3 2 e 3 β e 3 γ ∂ γ b β . (A6)Crossed terms containing only one power of u 2 or u 3 have been neglected, for reasons discussed further below. The force terms in Eq. (A6) find a clear meaning after being rewritten in terms of derivatives of the velocity field. To do this, decompose ∂ β u E α as tensor, Θ E the expansion scalar, h E α β the orthogonal projector and a E α β the acceleration tensor, defined through

2 = u 3 2 = u ⊥ 2 / 2 .E β +b α b β +e 2 α e 2 β +e 3 α e 3

 22223 All terms containing only one power of u 2 or u 3 have already been excluded from Eq. (A6). The decomposition of the identity,η αβ = -u E α u β , combined with u E α ∂ β u E α =0 then leads to the last term in the equation for γ . Defining h ⊥ αβ ≡ e 2 α e 2 β + e 3 α e 3 β , the corresponding term ∝ b α ∂ α B in the equation for u is obtained through

  u E β u E γ + b β b γ ∂ γ B . (A17)The latter equation can be rewritten for γ and u ,

  The details are provided in Appendix B. At time t, this distribution can be approx-

	imated as
	dN d ln p	∼	1 g t exp -	ln g	p	ln	f ∅ ln p + S p 2f + (ln pmax -ln p)
			+	ln pmax g	ln	f ∅ ln pmax + S p ln 2 pmax -ln 2 p	,

By "strong turbulence", it is meant here the large-amplitude regime in which the magnetic fluctuations are comparable to the mean field, δB B 0 . It also implies that, all along the cascade, the turbulence will be strong in the sense defined by Goldreich and Sridhar, i.e. governed by critical balance[START_REF] Goldreich | Toward a theory of interstellar turbulence. 2: Strong alfvenic turbulence[END_REF].

available from: http://turbulence.pha.jhu.edu/Forced MHD turbulence.aspx.

The inertial term a E • b also depends on the sign of the particle momentum; we neglect this (generally subdominant) term here, for the sake of the discussion.

Similar power laws have been recently reported for the curvature of magnetic field lines in Refs.[START_REF] Yang | Role of magnetic field curvature in magnetohydrodynamic turbulence[END_REF][START_REF] Yuen | Curvature of magnetic field lines in compressible magnetized turbulence: Statistics, magnetization predictions, gradient curvature, modes, and self-gravitating media[END_REF].

ACKNOWLEDGMENTS

The possibility to use the resources of the JH Turbulence Database (JHTDB), which is developed as an open resource by the Johns Hopkins University, under the sponsorship of the National Science Foundation, as well as the tools made available for public use, is gratefully acknowledged. This work has been supported in part by the Sorbonne Université DIWINE Emergence-2019 program. It is a pleasure to thank V. Bresci, A. Bykov, L. Comisso, C. Demidem, L. Gremillet, M. Malkov, L. Sironi, and L. Vlahos for insightful discussions.

Appendix A: PARTICLE ACCELERATION IN R / E

The following conventions are used. Greek indices refer to spacetime indices in the laboratory frame, while latin indices (a, b, . . .) are spacetime indices in the frame R / E in which the comoving electric field vanishes. Correspondingly, tensors indexed with greek indices are understood to be evaluated in the lab frame, while those indexed with latin indices are expressed in the R / E frame. Where confusion can arise, a prime is used to label comoving frame quantities. Hence B denotes the magnetic field in the comoving frame, while B (as the electric field E) that in the lab frame. These two are related through the usual (instantaneous) Lorentz transforms,

The quantities γ E = B/B and β E = E × B/B 2 represent the Lorentz factor and threevelocity of R / E in units of c; its four-velocity is written

Following [START_REF] Lemoine | Generalized Fermi acceleration[END_REF], the four-velocity of the particle in the locally inertial R / E frame, u a , evolves according to

with F ab the Maxwell field strength tensor in R / E (thus, purely magnetic), and Γ a b c the connection that captures the inertial forces,

here expressed as a function of the tetrad e β b (and its inverse) that connects the lab to the R / E frame, see e.g. [START_REF] Webb | Relativistic transport theory for cosmic rays[END_REF][START_REF] Webb | The diffusion approximation and transport theory for cosmic rays in relativistic flows[END_REF]. Note that the derivative is taken with respect to proper time τ which parametrizes the particle worldline.

The choice of this tetrad is not unique and, in the present case, it proves useful to define it as follows [START_REF] Boghosian | Covariant Lagrangian methods of relativistic plasma theory[END_REF][START_REF] Fradkin | Covariant electromagnetic projection operators and a covariant description of charged particle guiding centre motion[END_REF]

with u E µ the (timelike) four-velocity of R / E in the lab frame. The spacelike four-vector b µ is the unit vector defined by b µ = B µ / √ B α B α , with

the magnetic four-vector, which represents the magnetic field measured by an observer moving at velocity u E ν :

The differential operator D α denotes the orthogonally projected derivative relatively to the timelike direction set by u E α , see [START_REF] Tsagas | Electromagnetic fields in curved spacetimes[END_REF][START_REF] Barrow | Cosmology with inhomogeneous magnetic fields[END_REF],

and, in the MHD approximation, D α B α = 0 identically; that latter equation generalizes ∇ • B = 0 to accelerated frames. So far, everything has been written in terms of partial, noncovariant derivatives, but general covariantization is straightforward.

Connection to the guiding center approximation

The guiding center approximation relies on the conservation of the magnetic moment,

and it is best described (to lowest order) by the covariant relativistic Hamiltonian [START_REF] Littlejohn | Hamiltonian formulation of guiding center motion[END_REF],

in terms of u µ the four-velocity of the guiding center, defined as the particle four-velocity deprived of gyromotion,

The random walk described in Sec. IV B can be modeled as the adimensioned process

with the short-hand notations: y = ln(p/p inj )/ (N g), g the mean energy gain, N = t/τ the number of jumps, each taking time τ ; xi a random variable that takes value +1 with (energy gain) probability f + , -1 with (energy loss) probability f -and 0 with probability

The probability distribution function φ(ŷ) can be approximated using the techniques of large deviation theory [START_REF] Touchette | The large deviation approach to statistical mechanics[END_REF], as φ(ŷ) ∼