Particle acceleration in relativistic turbulence: A theoretical appraisal
Camilia Demidem Martin Lemoine Fabien Casse







 demidem@apc.in2p3.fr
Université  de Paris, CNRS/IN2P3, Astroparticule et Cosmologie, F-75006, Paris, France
Institut d'Astrophysique de Paris, CNRS – Sorbonne Université, 

98 bis boulevard Arago, F-75014 Paris, France
Nordita, KTH Royal Institute of Technology and Stockholm University,
 Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden
 
Institut d'Astrophysique de Paris, CNRS – Sorbonne Université, 

98 bis boulevard Arago, F-75014 Paris, France

Université  de Paris, CNRS/IN2P3, Astroparticule et Cosmologie, F-75006, Paris, France




  


We discuss the physics of stochastic particle acceleration in relativistic magnetohydrodynamic (MHD) turbulence, combining numerical simulations of test-particle acceleration in synthetic wave turbulence spectra with detailed analytical estimates. In particular, we study particle acceleration in wavelike isotropic fast mode turbulence, in Alfvén and slow Goldreich-Sridhar type wave turbulence (properly accounting for 
anisotropy effects), including resonance broadening due to wave decay and pitch-angle randomization. At high particle rigidities, the contributions of those three modes to acceleration are comparable to within an order of magnitude, as a combination of several effects (partial disappearance of transit-time damping for fast modes, increased scattering rate for Alfvén and slow modes due to resonance broadening). Additionally, we provide analytical arguments regarding acceleration beyond the regime of MHD wave turbulence, addressing the issue of nonresonant acceleration in a turbulence comprised of structures rather than waves, as well as the issue of acceleration in small-scale parallel electric fields. Finally, we compare our results to the existing literature and provide ready-to-use formulas for applications to high-energy astrophysical phenomenology.

DOI: https://doi.org/10.1103/PhysRevD.102.02300310.1103/PhysRevD.102.023003





Particle acceleration in relativistic turbulence: A theoretical appraisal
    Fabien Casse
    Received 27 September 2019; accepted 25 May 2020; published 1 July 2020
===========================================================================






§ INTRODUCTION



The physical mechanisms that govern energy dissipation and conversion in astrophysical sources and the acceleration of particles, from subrelativistic momenta to the most extreme observed cosmic-ray energies, are long-standing problems in space plasma physics and high-energy astrophysics. Magnetized turbulence plays a central role in this field, whether directly or indirectly. At the least, it provides the essential scattering agent at the core of diffusive shock acceleration <cit.>, and reconnection itself is intimately linked to turbulence, either because turbulence generates reconnection, e.g., Ref. <cit.>, or the converse <cit.>. Electromagnetic turbulence indeed provides an efficient source of particle acceleration, following Fermi's original idea <cit.> that a particle interacting with randomly moving magnetic mirrors of typical velocity dispersion β_m c  gains energy in a stochastic manner at a rate ∝β_m^2c^2. Such stochastic acceleration has been invoked to explain high-energy emission in a variety of astrophysical environments, from solar system plasmas to the remote high-energy Universe, e.g., impulsive solar flares <cit.>, the Galactic center <cit.>, accretion disks <cit.>, pulsar wind nebulae <cit.>, galaxy clusters <cit.>, active galactic nuclei <cit.>, gamma-ray bursts <cit.>, etc.

Astrophysical collisionless turbulence is commonly described as an energy cascade spanning orders of magnitude in length scales, from a large stirring scale down to the small wavelengths of dissipative physics, with, in general, most of the fluctuation power in velocity and electromagnetic fields being carried by the larger scales. The usual order-of-magnitude timescale for particle acceleration reads t_acc∼ t_scatt/β_m^2, where t_scatt represents the scattering timescale, i.e., the time needed for the particle to start diffusing in the turbulence. Hence, for the purpose of accelerating particles to high energies in cosmic plasmas, which is the main theme of this paper, one is interested in particles interacting with large-scale modes of a fast-moving turbulence spectrum: large scale, because the scattering timescale is a growing function of energy, and fast moving because of the scaling of the acceleration timescale. Here, we will thus be interested in the physics of particle interactions with relativistic ideal magnetohydrodynamic (MHD) turbulence, the ideal Ohm's law being applicable on the scales of interest. 

How particles gain energy in a turbulent setting can be addressed using a variety of theoretical or numerical tools. Quasilinear theory (QLT), e.g., Ref. <cit.>, provides an analytical estimate of the various diffusion tensor components to first order in the spectrum of electromagnetic fluctuations, which are commonly described as a sum of linear eigenmodes of the plasma. In the ideal MHD approximation, those are the incompressible Alfvén modes (thereafter indexed with A) as well as the fast and slow magnetosonic modes (indexed by F and S), while entropy modes, pure density perturbations advected with the medium, do not play any role. Recent nonlinear extensions allow to include part of the perturbative expansion in turbulent fluctuations by introducing corrections to the particle trajectory, e.g., Ref. <cit.> and references therein. An important outcome of quasilinear  calculations is the existence of resonant particle-wave interactions which provide possibly fast scattering rates, hence short acceleration timescales. However, whether collisionless turbulence can be realistically described as a sum of waves is a long-standing debate; see, e.g., Ref. <cit.> for a recent appraisal. Furthermore, the intrinsic anisotropic nature of modern MHD turbulence theories <cit.> prohibits  particle-wave resonances <cit.>. Hence, nonresonant phenomena, corresponding to the interaction of particles with nontrivial velocity structures, are likely to play a role, and our investigation will confirm this point of view.

On the numerical front, the physics of transport and acceleration can be probed i) by following test particles in a synthetic turbulence generated from a sum of plane waves, e.g., Refs. <cit.>; ii) by following test particles in full three-dimensional (3D) MHD simulations <cit.>; or iii), more recently, from 3D particle-in-cell (PIC) simulations <cit.>. While the latter method offers a fully kinetic picture of the collisionless turbulence, from plasma length scales upward, thus allowing, in particular, a self-consistent treatment of the early injection and acceleration stages, MHD simulations provide a useful representation of the largest length scales with a potentially large dynamic range. Although the first method i) remains subject to the criticism of describing the turbulence as a superposition of linear waves, it allows to relate in a direct way the theoretical predictions for acceleration to the assumptions of the model and to probe effects beyond quasilinear theory. It therefore represents an interesting tool to interpret the results of more evolved simulations, which remain expensive if there are to cover a reasonable dynamic range, and are, besides, possibly sensitive to some degree, to how the turbulence is initialized. 

Hence, the above analytical and numerical methods nicely complement each other and we adopt this stance in the present paper. More specifically, we combine theoretical arguments borrowed from quasilinear and extended quasilinear theories as applied to modern turbulence theories to derive predictions for acceleration in relativistic MHD wave turbulence, which we compare to test-particle simulations. With respect to previous work on this topic, our work improves on Ref. <cit.>, which considered an isotropic bath of Alfvén waves and neglected resonance broadening effects; it improves on Ref. <cit.>, which considered isotropic turbulence of pure fast, slow and Alfvén waves, by taking into account the effect of anisotropy and mode decay in the turbulence spectrum; it also improves on Ref. <cit.> by paying attention to the alignment of small-scale modes to the local large-scale magnetic field, which has dramatic consequences on the acceleration physics, as we show; we also improve on the previous semianalytical studies of Refs.<cit.> by considering the relativistic regime of wave phase velocities. Finally, we also provide original theoretical estimates for turbulent acceleration in the absence of resonant wave-particle interactions.

Our paper is laid out as follows. Section <ref> addresses the case of a pure fast mode turbulence, Sec. <ref> that of a pure slow mode turbulence, and Sec. <ref> that of pure Alfvén turbulence. For each, we compare analytical predictions to numerical simulations. In Sec. <ref>, we discuss some aspects of stochastic acceleration outside the realm of ideal MHD wave turbulence, i.e., we provide estimates for nonresonant acceleration and characterize the influence of small-scale violations of Ohm's law. Our numerical and analytical results are brought together in Sec. <ref> and compared to recent ab initio simulations of turbulent acceleration. Finally, we provide a summary in Sec. <ref>. Unless otherwise noted, we assume ultrarelativistic particles and use units such that c = 1.



§ FAST MODE TURBULENCE

Simulations of subrelativistic compressible MHD turbulence seem to indicate that Alfvén and slow magnetosonic modes are the main contributors to the kinetic and/or magnetic energy spectra, with little contribution from fast magnetosonic modes <cit.>, although recent work <cit.> indicates that this very result depends on how turbulence is driven at the outer scale. In any case, this does not exclude that fast modes can influence particle acceleration, depending on their scattering efficiency. Besides, simulations by <cit.> suggest that in the relativistic regime, fast modes are strongly coupled to Alfvénic modes and maintain a substantial fraction of energy, which furthermore tends to increase with magnetization.   

Given the variety of numerical setups used to simulate MHD turbulence, the size of the  parameter space to be probed (e.g., the β-parameter of the plasma, the magnetization, the strength of the turbulence, etc.) and the large uncertainty on the scaling laws obtained from simulations or observations, it is still not clear which of the available phenomenology, e.g., Refs. <cit.>, appears best suited to describe the properties of the cascade for these different modes. Here, we rely on the general result that emerges from theoretical considerations and simulations in the subrelativistic limit, which defines the cascade of fast modes as isotropic, with a power-spectrum index around 1.5–1.7, and an anisotropic cascade of Alfvén modes and (passive) slow modes, à la Goldreich-Sridhar.



 §.§ Theoretical predictions

For fast modes, we introduce the power spectrum

    𝒮_k^F = η/1-ηB_0^2 | 1-q_F|/8 k_min^3 (k/k_min)^-q_F-2,

properly normalized over the range of wave vectors k ∈ [k_min, k_max] with k_max ≫ k_min, to the amplitude of the magnetic perturbations contributed by fast modes, 2∫ d^3k 𝒮_k^F = η_F B_0^2/(1-η) = η B_0^2/(1-η), where the last equality holds for pure fast mode turbulence and the factor 2 accounts for the summation over positive and negative frequencies. Here, η=⟨ B^2⟩/(B_0^2 + ⟨ B^2⟩) characterizes the relative magnitude of the turbulent magnetic energy density.

The polarization of fast modes in special relativistic MHD is given by, e.g., Ref. <cit.>,

    ω_F   =±k/√(2){β_F⊥^2 +β_A^2 β_s^2μ_k^2   
         + [ ( β_F⊥^2 +β_A^2 β_s^2μ_k^2 )^2 -4 β_A^2 β_s^2μ_k^2   ]^1/2}^1/2,
    
    B_k^F   = b_k (-k_∥/k_⊥k_⊥/k +k_⊥/k_∥k_∥/k), 
    u_k^F   =ω_F/k  b_k/B_0 ( k_⊥/k_⊥ + β_s^2 k_⊥/ω_F^2-β_s^2k_∥^2k_∥), 
    E_k^F   =ω_F/k  b_k B_0×k_⊥/B_0 k_⊥,

where β_A, β_s and β_F⊥=(β_A^2+β_s^2-β_A^2 β_s^2)^1/2 denote the Alfvén speed, the sound speed and the fast speed, respectively. The wave vector k, with modulus k, is decomposed over a parallel and a perpendicular component with respect to B_0: k_∥=k_∥B_0/B_0, k_⊥=k-k_∥, while μ_k = k_∥/k. The (complex) perturbation amplitude b_k is related to the power spectrum through 𝒮_k^F = ⟨| b_k|^2⟩.

In our analytical predictions, we neglect the dependence on μ_k of the phase speed, β_F; in particular for β_Aβ_s≪β_F⊥, we use the simplified dispersion relation ω_F ≃ ±β_F⊥ k, which provides a good approximation over the range of μ_k ∈ [-1,1], while for β_A≈β_s, the value of β_F is set to the average over μ_k of the phase speed. These analytical predictions are obtained through quasilinear theory, which evaluates the second-order moments
⟨Δμ^2⟩ and ⟨Δ p^2⟩ by direct integration over the electromagnetic fields that the particle experiences over a trajectory which is described by the zeroth-order gyration around the background magnetic field; see Ref. <cit.> for details. Everywhere in this paper, μ represents the pitch-angle cosine of the particle with respect to the background magnetic field, and p represents its momentum; hence μ = p_∥/p.

A detailed application of quasilinear theory then provides the following form for the diffusion coefficients:

    D̃_μμ^F   =Ω^2 ^2(1-μ^2)/2B_0^2∑_±∑_n=-∞^+∞∫_k_min^k_maxdk
          ×∫_-1^1dμ_k ℛ_k[J_n+1(z_⊥)-J_n-1(z_⊥)]^2
          ×[(kμ_k - μω_F)^2+γ_ d^2μ^2ω_F^2] 𝒮_k^F,
    
    D̃_pp^F   =Ω^2^2p^2(1-μ^2)/2B_0^2∑_±∑_n=-∞^+∞∫_k_min^k_maxdk  
          ×∫_-1^1dμ_k ℛ_k[J_n+1(z_⊥)-J_n-1(z_⊥)]^2
          ×ω_F^2(1+γ_ d^2) 𝒮_k^F,

with the following notations: the sum over ± sums over positive and negative real frequencies ω_F, while that over n runs over the harmonics of the gyrofrequency; γ_d>0 accounts for the possible finite lifetime of modes, i.e., the full mode pulsation ω = ω_F - iγ_d|ω_F|; Ω = c/r_g = e B_0/p and  z_⊥ ≡ k_⊥Ω^-1√(1-μ^2). 

Finally, Eqs. (<ref>) and (<ref>) can be averaged over μ to obtain the quasilinear predictions for the diffusion coefficients for an isotropic population of particles, D_μμ and D_pp, linked by the relation

    D_pp∼ p^2 β_F^2 D_μμ.

Below, we provide the theoretical scalings for D_μμ; those for D_pp can be directly obtained through the above relation.

The resonance function ℛ_k characterizes the interaction between the particle and a given mode. Various forms for this function in standard and extended quasilinear theories, which account for resonance broadening through wave decay and partial randomization of the particle pitch-angle cosine, are introduced and discussed in Appendix <ref>. By wave decay, it is meant that modes are assigned a finite lifetime, which effectively implies a temporal correlation of finite extent. Here, we do not include this effect for fast modes, but we will comment on its possible influence. By contrast, the partial randomization of the pitch-angle cosine of the particle in a turbulent bath is guaranteed, because the direction of the total magnetic field does not coincide exactly with that of B_0; hence μ, which is defined relative to the latter, is effectively a random quantity (see Appendix <ref> for details).

Neglecting this effect for the time being, and considering furthermore γ_d→0, we have ℛ_k=δ(kμ_kμ - ω_F + nΩ),
with n∈ℤ,  bringing out the infinite discrete set of resonances of standard quasilinear theory. The n=0 resonance describes the Landau, also called transit-time damping (TTD) resonance, which results from the interaction of the particle with the compressive modes moving along the background magnetic field <cit.>. The n≠0 resonances are described as Landau-synchrotron, or gyroresonances, between the particle motion around the background field and the motion of the mode along the field. 

For the above Dirac resonance function, we derive the following quasilinear scalings in the limit r_gk_min≪1 (r_g=c/Ω the particle gyroradius),

    D_μμ^F-TTD   ∼ (1-β_F)^α η/1-η (r_gk_min)^q_F-2 k_min,
    
        D_μμ^F-Gyr   ∼0.1 η/1-η (r_gk_min)^q_F-2 k_min,

for the quasilinear TTD and gyroresonant contributions of the pitch-angle diffusion coefficient, where α≈ 4 is a (β_A-dependent) effective exponent. This scaling reveals that the transit-time damping contribution is strongly suppressed in the relativistic limit β_F→1. This is easily understood from the resonance condition,  kμ_kμ ± β_Fk = 0, which requires |μ|≥β_F, hence a superluminal parallel phase velocity when β_F≃1, thus effectively shutting down the resonance. 

Accounting for resonance broadening restores part of the transit-time damping contribution in the relativistic limit. To see this, consider the resonance broadening that results from pitch-angle randomization, with ⟨μ⟩∼⟨Δμ^2⟩^1/2∼η^1/4, as discussed in Appendix <ref>, for strongly magnetized particles, meaning r_gk_min≪1. Then ℛ_k∝exp[-(μ_kη^1/4 - β_F)^2/2η^1/2]. At large values of β_F, more specifically β_F> η^1/4 where the term in the exponential cannot vanish, the resonance function can be approximated by its value at μ_k=1 and a width β_F/η^1/4 in μ_k. One then obtains an estimate for  D_μμ^F-TTD that is comparable to, albeit sligthly smaller than, D_μμ^F-Gyr at β_F∼1. 

For β_F<η^1/4, the resonance kμ_k⟨μ⟩ - β_Fk = 0 can be met at values |μ_k|<1; hence, one recovers the standard quasilinear result, in which the transit-time damping contribution exceeds the gyroresonant ones by about 1.5 orders of magnitude.

For particles of large rigidity, meaning r_gk_min≫1, the scattering is dominated by gyroresonances because D_μμ^F-TTD∼η(r_gk_min)^-3k_min while 
D_μμ^F-Gyr∼η(r_gk_min)^-2k_min. The latter scaling is typical of a particle of gyroradius r_g interacting with a turbulence on coherence scales k_min^-1≪ r_g.



 §.§ Numerical simulations and discussion

Appendix <ref> presents in detail our numerical procedure for the Monte Carlo simulations of test-particle transport and acceleration in a given synthetic turbulence, described as a sum of linear waves. It appears important to recall here that the electric field is calculated using the ideal Ohm's law, see Eq. (<ref>), rather than as the sum of the E_k^F modes detailed above, in order to avoid the emergence of non-MHD effects <cit.>. 

Besides those parameters defining the turbulence, q_F, β_A, β_s, and η, the main (dimensionless) parameter that pilots the interaction of particles with waves is r_gk_min, with r_g larger than the smallest wavelength of the turbulence (r_g>k_max^-1), unless specified otherwise. 
We stress that these Monte Carlo simulations are, by default, carried out in the rest frame of the unperturbed background plasma, which does not generally coincide with that of the plasma including perturbations. This means, in particular, that the turbulence generally carries a nonzero bulk velocity, hence a net bulk electric field in the simulation frame. By contrast, quasilinear predictions are (implicitly) carried out in a frame in which there is no bulk electric field. In order to gauge the influence of this difference of frames, we have carried out additional simulations for which the initial data of our Monte Carlo simulations are instead specified in the local fluid rest frame, in which the motional electric field vanishes by construction (see Appendix <ref>). 
We will come back to this subtlety later.

To limit the size of parameter space, we restrict ourselves to the fiducial values: q_F=5/3, β_s≈10^-2β_A, η=0.3. The value of β_s exerts no influence on our results as long as β_s≪β_A, since β_F≃β_A in this limit. The quasilinear predictions for the diffusion coefficients directly scale with η/(1-η), and we have checked that this scaling holds in our numerical simulations.


We first show in Fig. <ref> the time evolution of various quantities for  β_A = 0.7 and r_gk_min≈ 10^-3.  This figure illustrates several noteworthy points about the acceleration of particles in a relativistic turbulence. The variance of the momentum distribution function σ_p^2 can be described, at early times, by a coherent oscillating pattern of energy gains and losses as particles gyrate around the background magnetic field and collect the influence of the perturbations. This modulation remains coherent over the coherence time of the random force that the particles suffer. Eventually, the stochastic buildup of the perturbations leads to the decorrelation of the pitch angle, after a time of the order of t_scatt, or around 0.3 k_min^-1 in Fig. <ref>. Particle acceleration occurs soon after, since the expected acceleration timescale t_acc∼ t_scatt/β_A^2. Once acceleration takes place, the behavior of σ_p^2/2t becomes strongly superdiffusive, because the diffusion coefficient D_pp is a growing function of energy. In detail,  quasilinear theory predicts D_pp^F∝ p^q_F, and the solution of the corresponding Fokker-Planck equation can be shown to exhibit in this case a momentum evolutionary law ⟨ p^2⟩^1/2∼⟨ p⟩∝Δ t^1/(2-q_F)∝Δ t^3. 

Hence, once acceleration starts to operate, the particle momentum increases at a fast rate until it reaches the point where r_gk_min≃1, at t∼ 100 k_min^-1 in Fig. <ref>. At larger values of the rigidity, quasilinear theory now predicts D_pp^F∝ p^0, which directly stems from the lack of resonances for  particles with gyroradius outside the inertial range of turbulence. Accordingly, at r_gk_min ≫ 1, the particle experiences a turbulence that it sees as small scale on its gyroradius scale, hence t_scatt∝ p^2, guaranteeing that t_acc∝ p^2, in agreement with D_pp^F=p^2/(2t_acc)∝ p^0.
The constancy of D_pp^F then indicates that σ_p^2/(2t) remains fixed at this plateau value, and one recovers normal diffusion of the momentum. 

Interestingly, as the momentum evolves fast in the acceleration process, the spatial diffusion coefficient, in particular the parallel one, also becomes strongly superdiffusive. We note here that quasilinear theory cannot describe this stage, because one of its intrinsic assumption is that particles move on unperturbed orbits at constant energy, which is clearly not the case here. This super-diffusive regime has important consequences for the maximum energy of accelerated particles, because it allows them to escape faster than on a naive diffusive timescale; this is discussed in more detail in Sec. <ref>.

We now compare the quasilinear predictions for the pitch-angle averaged diffusion coefficients D_μμ^F and D_pp^F over a range of rigidities and mode velocities, in Fig. <ref> and <ref>. In our simulations, particles are injected with random initial directions, hence random initial pitch angles. These diffusion coefficients are thus understood as averaged over μ. The details of the procedure used to derive them are given in Appendix <ref>; in short, D_μμ is estimated as the slope of the linear fit of ⟨[Δμ(t)]^2 ⟩/2 while D_pp is determined by a linear fit of σ_p^2(t)/2, as indicated by the dashed horizontal lines in the illustrative Fig. <ref>. Estimations of the acceleration and scattering timescales are also indicated as vertical lines, with t_acc taken as the time at which the average ⟨ p⟩ over the population of Monte Carlo particles is twice the initial value, and t_scatt is defined as the decorrelation time of the pitch-angle cosine, inferred through an exponential fit of ⟨μ(Δ t)μ(0)⟩. This procedure necessarily introduces a degree of arbitrariness, however, the scalings of the diffusion coefficients with respect to rigidity r_gk_min and mode velocity β_A (and with turbulence strength parameter η) are respected.







As discussed above, the transit-time damping contribution largely dominates the gyroresonant contributions at low β_A, but disappears, at least partially, in the relativistic regime β_A→ 1. This feature is illustrated in Fig. <ref>, which shows the (theoretical) relative  contributions of TTD and gyroresonances at two extreme values of the mode velocity, β_A=0.01 and β_A=0.9. In standard quasilinear theory, the TTD contribution has almost completely disappeared at β_A=0.9. 


Here, we present the theoretical predictions for quasilinear theory including only the resonance broadening associated to the uncertainty in pitch angle because our simulations indeed consider fast modes of an infinite lifetime. Our theoretical calculations reveal that the inclusion of resonance broadening does not affect much the diffusion coefficients in fast mode turbulence, although it will play an important role in the case of slow mode anisotropic turbulence. More precisely, its effect is negligible at values β_A≲ 0.1, because the TTD resonance is then strongly dominant; at values β_A≳ 0.7, it increases the standard QLT predictions by a factor of approximately 3–4, notably due to a partial restoration of the otherwise vanishing TTD contribution (Fig. <ref>). Nevertheless, we wish to point out that because these predictions with resonance broadening implicitly rely on the assumption that r_g k_min≪ 1, they tend to overestimate the diffusion coefficients for r_g k_min≳ 1. For these high rigidities, standard QLT is therefore better suited.

Overall, our numerical results for D_μμ, and for D_pp, appear in good agreement with the theoretical expectations of extended quasilinear theory. We, however, need to come back to some corrections introduced for simulations at β_A=0.7 and even more at β_A=0.9 because of the presence of a net bulk electric field which makes the comparison with our QLT calculations less direct for the following reasons.

This electric field is of typical strength ⟨δ E^2⟩^1/2≃𝒪(η^1/2β_A) and it is directly associated with a net bulk motion of the turbulence in the simulation frame, with Lorentz factor γ_u ≃(1+η u_A^2)^1/2. When particles are initialized with a random pitch angle and a momentum p_0 in the simulation frame, their energy changes by Δ p/p_0≃γ_u^2-1 once they interact with the turbulence, as a result of a first-order Fermi process. Consequently, their rigidity in the simulation frame is not the one that was initialized, r_g,0k_min, but rather r_gk_min=[1+Δ p/p_0] r_g,0k_min. In the relativistic regime, this effect becomes significant for D_pp, whose scaling with r_g is strong. More specifically, the values of the initial energy jump Δ p/p_0 measured in the simulations are approximately 0.6 for β_A=0.9 and 0.4 for β_A=0.7. We used these values to report the results in Fig. <ref> as a function of the effective initial r_gk_min. These values are represented by the lighter symbols in Fig. <ref> and they lie above the theoretical predictions.

Given that there is no bulk electric field involved in our quasilinear computations, one can expect the local rest frame of the turbulence to be better suited to set up the initial conditions of our Monte Carlo simulations. We therefore also conducted simulations, which initialize the energy distribution as a Dirac delta function in that frame, and their results are represented by opaque symbols in Figs. <ref> and <ref>. The corresponding values of D_pp appear to be in slightly better agreement with the quasilinear predictions, which, as mentioned earlier, also account for resonance broadening (although at high rigidities, standard QLT is better justified). 

Finally, it is noteworthy that in the relativistic limit, the D_pp values that we measure here appear to scale as a function of the 3-velocity (rather than the 4-velocity), D_pp∼β_F^2.




§ SLOW MODE TURBULENCE




 §.§ Theoretical predictions


From a formal standpoint, the analysis for slow modes is the same as for fast modes. The expressions for the perturbation polarizations, Eqs. (<ref>)–(<ref>) remain valid, and one has simply to replace ω_F with 

    ω_S   =±k/√(2){β_F⊥^2 +β_A^2 β_s^2μ_k^2  
         - . [ ( β_F⊥^2 +β_A^2 β_s^2μ_k^2 )^2 -4 β_A^2 β_s^2μ_k^2   ]^1/2}^1/2.

Likewise, the diffusion coefficients adopt the same form as before, Eqs. (<ref>) and (<ref>), except that we now adopt a Goldreich-Sridhar power spectrum,


    𝒮_k^S =
     η/1-ηB_0^2  C_g (k_⊥/k_min)^-q_S-1g[k_∥/k_⊥^2/3k_min^1/3],

with g(x) a function peaking at -1≲ x ≲ 1, imposing the anisotropic scaling between the parallel and perpendicular components of modes and with C_g, a normalization constant ensuring that ∑_±∫ d^3k 𝒮_k^S= η_S B_0^2/(1-η)=η B_0^2/(1-η). We have considered two types of anisotropy functions, g(x)=Θ(1-| x|), for which  C_g=(q_S-1)(3q_S-5)/(16 k_min^3), and g(x)=δ(√(2)x-1) for which C_g=√(2)(3q_S-5)/(24 k_min^3) and found that both functions lead to similar results. To simplify our analytical calculations, an approximate dispersion relation is used again, ω_S≈±β_S k_∥.


We now derive the predicted analytical scalings, starting with the case of standard quasilinear theory, without resonance broadening. We first recall that even though slow modes carry compressive perturbations, the TTD contribution is null in standard quasilinear theory as the linear resonance condition reads ± k_∥β_S-k_∥μ=0 and is virtually never met. Using Eqs. (<ref>) and (<ref>) and inserting the above power spectrum (with q_S=7/3), one derives the following scalings for gyroresonant interactions at small rigidities, i.e., r_gk_min≪1,

    D_μμ^S-Gyr∼ 0.1η/1-η(r_gk_min)^3/2 k_min,

and D_μμ^S-Gyr∼ 0.1η(r_gk_min)^-2 k_min as usual, in the high-energy limit r_gk_min≫1. The strong scaling D_μμ^S-Gyr∝(r_gk_min)^3/2 stems from the anisotropy of the Goldreich-Sridhar spectrum <cit.>, since k_⊥≫ k_∥ implies that resonant particles with r_gk_∥∼ 1 will cross many uncorrelated fluctuations in the perpendicular direction. 

We now turn to the more physically motivated case of damped modes. The condition of critical balance underlying Goldreich-Sridhar model implies that the timescale of nonlinear interactions coincides with the linear propagation timescale; accordingly, we set γ_d=1 in numerical applications. In the presence of time decorrelation, the resonance function takes on a Breit-Wigner form; see Eq. (<ref>) with ω=ω_S. A direct calculation of the integrals then provides the following scalings in the low-rigidity limit r_gk_min≪1:

    D_μμ^S-TTD   ∼0.1η/1-ηγ_dβ_S[1 - 3ln(r_gk_min)^] k_min,
     
    D_μμ^S-Gyr   ∼0.1η/1-ηγ_dβ_S k_min.


In agreement with previous studies carried out in the nonrelativistic limit <cit.>, we find that the resonance broadening allows TTD interactions, which now dominate the transport, and provide for q_S=7/3, an interaction time with a weak logarithmic dependence on the rigidity[Ref. <cit.>, which uses an exponential decay model for the turbulence decorrelation similar to ours, obtains the same logarithmic dependence for the diffusion coefficients; our numerical values are, however, larger by about an order of magnitude.] for r_gk_min≪1.  We also note that D_μμ^S scales with β_Sγ_d which characterizes the scaling of the decay coefficient.

In the high-rigidity limit, meaning r_gk_min≫1, D_μμ^S-TTD∼ηγ_dβ_S(r_gk_min)^-3 k_min while D_μμ^S-Gyr∼ 0.1 η(r_gk_min)^-2 k_min.  

Consider now the effect of the partial randomization of the pitch angle of the particle with respect to the direction of the magnetic field, as discussed in Appendix <ref>. To distinguish this effect from that associated to the finite lifetime of the modes, we assume here γ_d=0. Using Eq. (<ref>),  we thus derive, for r_gk_min≪1,

    D_μμ^S-TTD∼ 0.1κ(η,β_S) η k_min [1 - 1.8ln(r_gk_min)],

with κ(η,β_S)∼exp[-(1-β_S/η^1/4)^2] a function of η and β_S that stems from resonance broadening, and which is such that κ∼𝒪(1) for β_S≲η^1/4. In the opposite limit, κ drops exponentially fast toward zero as η^1/4/β_S→ 0, in agreement with the vanishing TTD contribution for slow modes in the absence of resonance broadening. 

The gyroresonant contribution is here strongly suppressed relative to its TTD counterpart, and can be neglected safely. This difference, with respect to resonance broadening by wave decay, Eqs. (<ref>) and (<ref>), results from the shape of the resonance function, which is exponential in the present case, vs Breit-Wigner in the previous case. One can indeed show that the contribution to D_μμ^S-Gyr in Eq. (<ref>) results from modes with 
k_⊥≲ k_min(r_gk_min)^-3/2, i.e., k_∥ r_g≲ 1. However, imposing this constraint in the resonance function given by Eq. (<ref>) results in an exponential suppression of the form exp{-[n/(k_∥ r_gη^1/4)-…]^2}, where … represent terms of the order of unity, and n≥ 1 for gyroresonant interactions.

At rigidities below the inertial range, i.e., such that r_g k_max<1, the [1-1.8ln(r_gk_min)] prefactor is replaced by ln(k_max/k_min), and the D_μμ tends to a constant, as observed in Ref. <cit.>. 
In the high-rigidity limit, r_gk_min≫ 1, D_μμ^S-TTD ∼ κ(η,β_S) η(r_gk_min)^-3 k_min.



 §.§ Numerical simulations and discussion

Again, we relegate the details of the numerical setup to Appendix <ref>. However, we emphasize here that setting up the Goldreich-Sridhar phenomenology is not a trivial task as one must pay attention to the notion of a local mean field  <cit.>; the direction of anisotropy, which establishes the hierarchy between k_∥ and k_⊥, must be defined relative to the local mean field, which is the total field at the given point averaged over scales larger than the particle gyroradius. In order to account for this effect, we artificially reduce the amplitude of the turbulence down to low values, η∼0.01, which guarantees that everywhere, the field line direction does not depart from that of B_0  by an angle larger than B/B_0∼η^1/2. Then, provided that the eddy anisotropy verifies k_∥/k≳ B/B_0, this eddy can be considered as aligned with respect to the local mean field B, even if in practice, it is aligned along B_0. The interest of setting the parallel direction along B_0 is, of course,  to preserve a standard Fourier decomposition throughout space for our numerical simulations. Since particle-wave resonances require k_∥ r_g∼1, with k_∥/k∼ k_∥^-1/2k_min^1/2 (critical balance), the above constraints now imply r_gk_min≳η.

The above setup with η=0.01 thus allows us to probe the physics of particle acceleration in a realistic anisotropic Goldreich-Sridhar-like configuration for particles with r_gk_min≳ 0.01. Since the analytical scalings predict D_μμ∝η/(1-η) and D_pp∝η/(1-η), we may, in turn, extrapolate these scalings to the regime of larger η provided they match the numerical simulations. To ensure that β_S remains close to unity, we consider a hot plasma with β_s=1/√(3) and β_A=0.5. For β_s≪β_A or β_A≪β_s, we have indeed β_S≃min(β_s, β_A). 





Figures <ref> and <ref> compare the generalized quasilinear predictions for the pitch-angle averaged diffusion coefficients described in Sec. <ref> to the values derived from the Monte Carlo simulations. For dynamical turbulence (γ_d=1), we find good agreement between the simulations (red symbols in the figures) and the semianalytical predictions, with D_μμ∼ cst up to a logarithmic correction ln(r_gk_min).

For undamped modes (γ_d=0), however, the scattering time in Monte Carlo simulations (orange symbols) remains comparable to that for damped modes and several orders of magnitude larger than the standard quasilinear estimates (blue curves), which furthermore predict D_μμ∝ (r_gk_min)^3/2. The momentum diffusion coefficient D_pp reveals a similar discrepancy. These results are, however, in satisfactory agreement with our theoretical predictions given in Eq. (<ref>), which include resonance broadening by pitch-angle randomization. Although this effect lies beyond quasilinear theory, which considers unperturbed orbits, such partial randomization is inherent to our numerical Monte Carlo simulations, because they follow the trajectories of the particles in the exact electromagnetic field configuration. Our test-particle simulations also indicate that the diffusion coefficient depends on the pitch-angle cosine, being more pronounced at values μ≃±0.5 than at μ≃±1, as expected from the broadening of the TTD resonance described above. This suggests that our theoretical scaling given in Eq. (<ref>) captures the physics that the test-particle simulations reproduce. 

Note the scale in the diffusion coefficients, in particular that of D_μμ, which suggests t_scatt∼10^2k_min^-1 at r_gk_min≪1. This low value of D_μμ (equivalently, this large value of t_scatt) is a direct consequence of our choice η=0.01, since D_μμ∝η/(1-η) (hence t_scatt∝η^-1 for η≪ 1); see Eqs. (<ref>), (<ref>) and (<ref>). Hence, extrapolating to η∼1 predicts a scattering timescale of the order of k_min^-1. When adapting the values of D_pp, one must also keep in mind that those shown in Fig. <ref> assume η=0.01 and that D_pp∝η/(1-η).



§ ALFVÉN MODE TURBULENCE




 §.§ Theoretical predictions

We  use the same Goldreich-Sridhar power spectrum as for slow modes, Eq. (<ref>), and repeat the generalized quasilinear calculations of Sec. <ref> with dispersion and polarization relations appropriate for Alfvén modes,

    ω_A   =±β_A k_∥, 
    B_k^A   =  b_k B_0 ×k_⊥/B_0 k_⊥, 
    u_k^A   =∓β_AB_k^A/B_0,
    E_k^A   =∓β_A b_kk_⊥/k_⊥,

to obtain

    D̃_μμ^A   =Ω^2^2(1-μ^2)/2B_0^2∑_±∑_n=-∞^+∞∬d k_⊥d k_∥  k_⊥
          ×ℛ_k[J_n+1(z_⊥)+J_n-1(z_⊥)]^2 k_∥^-2
          ×[(k_∥ - μω_A)^2+μ^2γ_d^2ω_A^2] 𝒮_k^A,
           
    D̃_pp^A   =Ω^2^2p^2(1-μ^2)/2B_0^2∑_±∑_n=-∞^+∞∬d k_⊥d k_∥  k_⊥
          × ℛ_k[J_n+1(z_⊥)+J_n-1(z_⊥)]^2k_∥^-2
          ×ω_A^2(1+γ_d^2) 𝒮_k^A.

The expressions differ from those derived for magnetosonic modes through the relative sign of the Bessel functions, which incidentally ensures that the n=0 term is null (J_-1=-J_1), as well as through the dispersion relation and power spectrum of the modes (which differs from that of fast modes), of course. But besides the fact that there are no TTD interactions, as a result of  the absence of compressible magnetic perturbations, the general remarks made for slow modes in Sec. <ref> remain true.

For the idealized case of undamped waves, γ_d=0, we derive the following theoretical scaling for r_gk_min≪1,

    D_μμ^A-Gyr∼ 0.1η/1-η(r_gk_min)^3/2 k_min,

similar to that obtained for slow mode waves in the same case. At high rigidities, r_gk_min≫1, D_μμ^Gyr∼ 0.1η(r_gk_min)^-2 k_min, as always.

If one now accounts for a finite lifetime of the waves, we obtain, for r_gk_min≪1,

    D_μμ^A-Gyr∼ 0.1 η/1-ηγ_dβ_A[1 - 2.7ln(r_gk_min)^] k_min,

and that for r_gk_min≫1 remains unchanged.

By contrast with slow modes, the resonance broadening that results from the partial randomization of the pitch angle does not play any significant role here, for two essential reasons. First, Alfvén waves do not possess a magnetic perturbation parallel to the mean magnetic field, so that the broadening is of relative order η^1/2, instead of η^1/4 for magnetosonic modes. Consequently, the resonance remains narrow if η≪1. Second, as discussed in the case of slow modes, this source of broadening does not affect as strongly gyroresonant interactions as it does enhance TTD contributions. For this reason, we do not anticipate any effect from the partial randomization of the pitch angle, and our numerical simulations will confirm this result.



 §.§ Numerical simulations and discussion







In our Monte Carlo simulations of Alfvén turbulence, the Goldreich-Sridhar local anisotropy is enforced using the same technique as for slow modes: we artificially reduce the amplitude of the turbulence down to a low level, η=0.01, which allows us to probe gyroresonant interactions down to r_gk_min≃0.01.

Figures <ref> and <ref> compare the values of the diffusion coefficients obtained from numerical estimations of Eqs. (<ref>) and (<ref>) to those derived from our Monte Carlo simulations for the fiducial parameters η=10^-2, β_A=0.9. 


In these figures, the blue symbols (respectively, the blue dashed lines) correspond to the Monte Carlo simulation results (respectively, the quasilinear predictions) for undamped waves. They are found to agree fairly well one with the other, thus confirming the theoretical scalings of this (idealized) Goldreich-Sridhar phenomenology. We note, in particular, the difference with respect to slow mode waves, for which a non-negligible scattering frequency was observed in our test-particle Monte Carlo simulations for this case, because of the broadening of the TTD contribution by pitch-angle randomization.

For damped modes, we find again reasonable agreement for the r_gk_min-scaling, while Figs. <ref> and <ref> support the ∝β_A and ∝β_A^3 scaling of D_μμ and D_pp, respectively.

Finally, we note here as well that the values of D_μμ and D_pp take comparatively low values because of our choice η=0.01. To extrapolate these values to larger values of η, one must keep in mind that both scale linearly with η/(1-η); see Eqs. (<ref>) and (<ref>). In particular, for η∼1, one expects a scattering timescale of the order of k_min^-1 at r_gk_min≪1.



§ BEYOND MHD WAVE TURBULENCE




 §.§ Nonresonant turbulent acceleration


So far, we have discussed particle acceleration through resonant (or quasiresonant) interactions with MHD waves. Yet, it is not obvious that such linear eigenmodes of the plasma provide a faithful description of actual modes in strongly interacting turbulence. In this section, we address some salient features of the nonresonant acceleration of a particle in a (relativistic) turbulent bath.

On a general level, this physics can be described as the interaction of a particle with a fluctuating electric field that is directly associated through ideal Ohm's law to the fluctuating velocity field of the plasma. Particles then gain energy as they experience the compressive, shearing, accelerating and vortical motions of the medium, e.g. <cit.>. 

To characterize the acceleration rate, one needs to specify the statistics of the velocity field and how the particle experiences the resulting electric fields, i.e., how it is transported across the cells of coherence of the turbulence. At large rigidities, meaning r_gk_min≫ 1, particles cross a coherence cell in a near-ballistic manner; hence, the scattering timescale can be directly expressed as t_scatt≃η^-1 k_min^-1(r_gk_min)^2. For magnetized particles, i.e., r_gk_min≪ 1, the problem is, however, much more complex and lies beyond the scope of the present study. Here, we rather aim to discuss how a given scattering timescale t_scatt impacts the acceleration process. This scattering timescale t_scatt may be smaller or larger than k_min^-1, with important consequences for the physics of acceleration.

A key point, indeed, is that the particle experiences in a different way the modes on large scales k^-1, i.e., such that t_scatt≲ k^-1, and the small-scale modes for which k^-1≲ t_scatt <cit.>. Specifically, for a turbulent spectrum strongly peaked on a single scale K^-1, the diffusion coefficient obeys

    D_pp   ≃ p^2⟨ (∂ u)^2⟩_K t_scatt   (t_scattK≪1) 
       ≃ p^2⟨ (∂ u)^2⟩_K/K^2 t_scatt    (t_scattK≫1).

In the above formula, the quantity ⟨ (∂ u)^2⟩_K represents, in a symbolic way, the contribution of the compressive, shearing, accelerating and vortical motions on scale K^-1, to the acceleration process. In detail, the decomposition of the random fluid 4-velocity in terms of acceleration (characterized by a 3-vector a), compressive (characterized by a compression scalar θ), shearing (3-tensor σ^ij), and vortical (3-tensor ω^ij) motions leads to

    ∂_α u_i = δ^0_α a_i+δ^j_α( 1/3 θδ_ij+σ_ji+ω_ji),

with: θ = ∇·u, σ_ij = (∂_i u_j + ∂_j u_i)/2-1/3 δ_ij θ, and ω_ij = (∂_i u_j - ∂_j u_i)/2, i,j∈{1,2,3} and α∈{0,1,2,3}. In the subrelativistic limit ⟨ u^2⟩≪ 1 (u is here understood as a 4-velocity), only the compressive and shearing terms contribute to leading order. In the highly relativistic limit, however, all terms contribute roughly equally; see Ref. <cit.> for details. 

We now generalize these results, in particular Eq. (<ref>), to the case of a broadband spectrum of 4-velocity fluctuations extending over the range [k_min, k_max], characterized by a one-dimensional (1D) spectral index 1<q_u<2. We note that Ref. <cit.> has studied a similar problem, using a different formalism for turbulent shear acceleration in the subrelativistic limit. 

We therefore write

    ⟨ u^2⟩ = ∫_ln k_mindln k ⟨| u_k|^2⟩_k,

with ⟨| u_k|^2⟩_k ∝ k^1-q_u the typical 4-velocity perturbation on scale k^-1. We also approximate: 

    ⟨(∂ u)^2⟩_k ≡α_u k^2 ⟨| u_k|^2⟩_k,

with α_u a factor of the order of unity that depends on the properties of the turbulence. For instance, Alfvén waves do not contain a compressive term, while magnetosonic modes do.

If k_max^-1≪ t_scatt≪ k_min^-1, one can split the turbulence in a large-scale and a small-scale cascade around t_scatt, so that the diffusion coefficient receives two contributions (of comparable magnitude):

    . D_pp|_k^-1≳ t_scatt   ≃     p^2 t_scatt⟨ u^2⟩ k_min^2 
          ×∫^ln t_scatt^-1 dln k (k/k_min)^3-q_u
       ≃    p^2 (t_scattk_min)^q_u-2⟨ u^2⟩ k_min,
          
    . D_pp|_k^-1≲ t_scatt   ≃     p^2 t_scatt^-1⟨ u^2⟩
          ×∫_ln t_scatt^-1 dln k (k/k_min)^1-q_u
       ≃   p^2(t_scattk_min)^q_u-2⟨ u^2⟩ k_min.

Interestingly, the total contribution scales as

    D_pp∝ p^2-ϵ⟨δ u^2⟩,
 
with ϵ positive but small compared to unity, because t_scatt usually has a mild dependence on p, while q_u-2 is negative and small in magnitude. Using for instance the quasilinear scaling for isotropic turbulence, t_scatt∝(r_gk_min)^2-q, we find D_pp∝ p^2(r_gk_min)^-(2-q)(2-q_u).
One may furthermore note that the above result, Eq. (<ref>), departs from the naive scaling t_acc=p^2/D_pp≃ t_scatt/⟨ u^2⟩, because the particle is now sensitive to the detailed structure of the turbulent power spectrum. Indeed, one now obtains

    t_acc≈k_min^-1(t_scattk_min)^2-q_u/⟨ u^2⟩.


If t_scatt≳ k_min^-1, however, the above integration procedure now gives 

    .D_pp|_k_min^-1≲ t_scatt≃ p^2 ⟨ u^2⟩/t_scatt.

In particular, for r_gk_min≫1, we have t_scatt≫ k_min^-1 and t_scatt∝ r_g^2; hence, we recover D_pp∝ p^0.




 §.§ Non-MHD electric fields

Ideal Ohm's law E =-β×B prevents the existence of electric fields aligned with the local magnetic field and guarantees that there exists at every point a reference frame in which the electric field can be screened out (the local plasma rest frame). 

In the presence of a parallel electric field, acceleration could proceed at a much faster rate, unhindered by the magnetic field. On the largest physical scales, however, it is believed that the ideal Ohm's law provides a satisfactory approximation and our simulations have implemented this constraint. Ideal MHD rather breaks down on small length scales, generating parallel electric fields in reconnecting current sheets <cit.>, or, more generally, because of inertial and kinetic effects, e.g., Ref. <cit.>. Of course, if the typical current sheet width is l_⊥, then particles with gyroradius r_g≫ l_⊥ stream through the sheet with small deflection and energy gain; hence, such reconnection processes govern the small-scale physics of dissipation, heating, and injection, but not acceleration to high energies.

Our aim, here, is to quantify the influence of these small-scale violations of Ohm's law on the acceleration process. We first characterize the statistics of these electric field fluctuations through a power spectrum, assuming that they extend on spatial scales greater than k_max^-1 but peak on k_max^-1 with, correspondingly, a 1D index q_E<1 (q_E can possibly take negative values). In Appendix <ref>, we calculate the relevant index from the corrections to Ohm's law for Alfvén, fast and slow wave turbulence respectively, for pair and electron-proton plasmas. 

The general problem can be brought back to the one-dimensional analog that describes the evolution of the particle momentum through

    dp/dt = q μ E_∥,

where E_∥ is understood here as a random field and μ represents the effective (random) velocity of the particle, characterized by a time correlation function of step t_scatt. Note that this scattering timescale may differ from that introduced in the previous sections, as it may be affected by the small-scale electric field itself. 

Equation (<ref>) describes a diffusion process, because over a time interval Δ t, ⟨Δ p⟩=0 but ⟨Δ p^2⟩≠0. 

For explicit calculations, we adopt for the velocity correlation function

    ⟨μ(t)μ(0)⟩=μ(0)^2exp[-/4t^2/t_scatt^2].

As it should, the scattering timescale verifies t_scatt=∫_0^+∞dt ⟨μ(t)μ(0)⟩/μ(0)^2.

In the spirit of the previous sections, we decompose the electric field fluctuations over scales k^-1  as a power-law spectrum and we further assume that at each such scale, the coherence length of the corresponding mode is also k^-1. This generalized decomposition in wave packets mimics the decomposition of the damped modes of the Alfvén and slow modes.  Hence, we use

    ⟨ E_∥(t_1) E_∥(t_2)⟩=α_E ⟨ E_∥^2⟩∫dk k^-q_Eexp[-k | t_1-t_2|].

The exponential term involves a correlation time k^-1, which can be seen as the turnover timescale of the parallel electric field structure in the turbulence.  Note that, using linear or squared expressions in the exponentials, either for Eq. (<ref>) or for Eq. (<ref>), would modify our final results by a factor of the order of unity only. The prefactor α_E=|1-q_E| k_max^q_E-1 provides the correct normalization to the equal-time (equal-position) amplitude ⟨ E_∥^2⟩.

We thus derive the momentum diffusion coefficient as

    ⟨Δ p^2⟩    =
        2Δ t q^2⟨ E_∥^2⟩α_Eμ(0)^2 t_scatt
         ×∫dk k^1-q_E ^k^2 t_scatt^2 / Erfc[k t_scatt/√()].

The complementary error function is here defined as Erfc(x)=(2/√())∫_x^+∞dt ^-t^2. Integrating this contribution over k then gives the approximate diffusion coefficient D_pp^ E_∥=⟨Δ p^2⟩/2Δ t:

    D_pp^ E_∥   ≈ q^2⟨ E_∥^2⟩ t_scatt   (t_scatt≪ k_max^-1) 
       
       ≈ q^2⟨ E_∥^2⟩( t_scattk_max)^q_E/k_max   (t_scatt≫ k_max^-1,  0<q_E) 
       
       ≈ q^2⟨ E_∥^2⟩ k_max^-1   (t_scatt≫ k_max^-1,  q_E<0).

This diffusion coefficient thus increases up as the scattering timescale until the particle starts to see these fluctuations as small scales, meaning t_scatt≳ k_max^-1, at which point it may either level off, if q_E<0 (corresponding to a spectrum that is sharply peaked at k_max), or increase as t_scatt^q_E, if the long-wavelength fluctuations at k^-1>t_scatt retain enough power, meaning 0<q_E<1. We note that, in Eq. (<ref>), the first approximation indeed describes the diffusion coefficient of a particle changing direction every t_scatt while traveling in a roughly coherent electric field, while the third approximation describes how a particle gains energy by traveling ballistically over many incoherent patches of parallel electric field of typical scale k_max^-1. In the intermediate limit, the particle feels the structure of the power spectrum of electric field fluctuations. For reference, we derive in Appendix <ref>, q_E≃ -7/3, -5/3, -1 respectively, for isotropic fast mode, Goldreich-Sridhar slow mode, and Goldreich-Sridhar Alfvén mode turbulences, from the deviations to the ideal Ohm's law in a pair plasma. 

Provided t_scatt^q_E roughly scales less fast with energy than a square law, acceleration by the large-scale MHD turbulence, as described in the previous sections, eventually takes over at some energy which may be easily calculated from the above expressions, once t_scatt(p), q_E, and ⟨ E_∥^2⟩ are specified. To ease the comparison, we note that the units in which we express D_pp in Figs. <ref>, <ref>, and <ref> are q^2B_0^2/(k_minc), while those for D_pp^ E_∥ are  q^2⟨ E_∥^2⟩/(k_maxc).



§ DISCUSSION



The physics of particle acceleration in a turbulent setting is governed by a variety of effects, depending on whether the cascade can be described as isotropic or not, whether it can be approximated as linear waves or not, whether one assumes ideal MHD to hold or not. In this section, we recap and bring together the results obtained in the previous sections and compare them to the results of recent first-principles extensive numerical simulations of turbulent acceleration.



 §.§ Comparison between modes and general results

In an actual turbulent setting, one expects a mixed contribution from various modes. As mentioned earlier, numerical MHD simulations of subrelativistic turbulence generally point to the dominance of Alfvén and slow modes, with a minor contribution from fast magnetosonic modes <cit.>. This result, however, appears to depend on how the turbulence is driven at the outer scale <cit.>, and the phenomenology may differ in the  relativistic regime <cit.>. We can nevertheless combine our various results, adopting the notations η_F, η_A, and η_S for the respective contributions of the various modes to the total magnetic energy density.

Our numerical simulations for acceleration in fast, Alfvén and slow mode turbulences provide the following scalings in the inertial range r_gk_min<1,

    D^F_pp≃ 0.5η_F/1-η(r_gk_min)^q_F(β_A/0.9)^2, 
       D^A_pp≃ 0.07 η_A/1-η(r_gk_min)^2[1-2.7ln(r_gk_min)](β_A/0.9)^3,
       D^S_pp≃ 0.01 η_S/1-η(r_gk_min)^2 [1-1.8ln(r_gk_min)](β_S/0.5)^2,

where all diffusion coefficients are here written in units of m^2 Ω_0^2/k_min, with Ω_0=eB_0/m the cyclotron frequency. Note that m^2 Ω_0^2/k_min=p_conf^2 k_min, with p_conf the confinement momentum such that r_g(p_conf)k_min=1.
We rely here on the (more realistic) picture of damped modes for Alfvén and slow modes, and we have used q_F=5/3, q_S=q_A=7/3 in our simulations. The wave number k_min is related to the outer scale of the turbulence through k_min=2π/L_max. Note that the effective slow mode phase speed β_S cannot be realistically larger than about 0.5 because it is bounded by the sound speed, which is itself bounded by 1/√(3)≃0.58 in a relativistically hot plasma. Note also the scaling D_pp∝β_A^3 for Alfvén modes, to be contrasted with the naive square law. Here, one extra power comes from ω∝ω, which controls the scattering through gyroresonance broadening. We have omitted this contribution for slow modes, since our analysis has revealed that pitch-angle randomization provides an equally strong source of scattering. Finally, in the region r_gk_min≫1, these scalings are each continued into a constant value D_pp∝η m^2k_min.

The acceleration timescale is conveniently written as t_acc=(r_gk_min)^2 (m^2Ω_0^2k_min^-1/D_pp)k_min^-1. We thus derive the following scalings:

    t^F_acc≃ 2.11-η/η_F (r_gk_min)^2-q_F(β_A/0.9)^-2k_min^-1,  
       t^A_acc≃ 14 1-η/η_A[1-2.7ln(r_gk_min)]^-1(β_A/0.9)^-3 k_min^-1, 
       t^S_acc≃ 100 1-η/η_S[1-1.8ln(r_gk_min)]^-1(β_S/0.5)^-2 k_min^-1.

If one rather defines the acceleration timescale as the timescale for the mean energy to double in the Monte Carlo simulations, the prefactors in Eq. (<ref>) are found to be instead approximately equal to 2, 4, and 30 for fast modes, damped Alfvén modes, and damped slow modes, respectively.

In Sec. <ref>, we have also discussed the physics of acceleration of particles in a more generic turbulent setting, which is not described as a sum of linear waves but whose velocity field is decomposed into a sum of compressive, shearing, vortical, and accelerating motions. Although the scattering timescale cannot be predicted on general grounds in such a situation, we have found that, for t_scatt ≲ k_min^-1, D_pp∼ p^2 ⟨δ u^2⟩ k_min (t_scattk_min)^2-q_u in terms of the power spectrum index of the 4-velocity fluctuations, q_u. We have argued there that, since t_scatt is generically a mild (increasing) function of p, and since 2-q_u is small compared to unity, D_pp can be written in the general form 

    D_pp∼ p^2 ⟨δ u^2⟩ k_min,

i.e., D_pp∼⟨δ u^2⟩(r_gk_min)^2 in units of  m^2Ω_0^2/k_min. The above estimate is valid up to a correction of the order (r_gk_min)^-ϵ, with ϵ significantly smaller than unity. Note that the amplitude of the turbulence, η, is included here in the 4-velocity fluctuation amplitude ⟨δ u^2⟩.





These different contributions are brought together in Figs. <ref> and <ref> for the case of a relativistically hot plasma with magnetization of order unity for η_F,S,A=η=0.3. 
Fig. <ref> reveals that, for β_A≃ 1, the contributions of fast and Alfvén  modes are roughly comparable at large rigidities, r_gk_min≳ 0.1,[Especially given that the predictions for fast modes displayed here include resonance broadening due to the pitch-angle uncertainty which tend to slightly overestimate the diffusion coefficients at large rigidities.] but that the fast mode contribution dominates at lower rigidities due to its slightly softer dependence on r_g. Recall however that this dependence scales directly with q_F, and if q_F∼2, the fast mode would no longer dominate over Alfvén modes.

In a relativistically hot, magnetized plasma, the phase speed of the slow modes can reach mildly relativistic values; hence, slow mode acceleration becomes truly efficient as well. In particular, we find that, at equal (relativistic) phase velocities, the diffusion coefficient in slow modes is only a factor of a few below that of Alfvén modes (Fig. <ref>). This hierarchy differs strongly from what is observed in the subrelativistic regime, namely a strong dominance of fast modes over Alfvén modes, with a negligible contribution from slow modes. This arises as a combination of several effects, notably the partial disappearance of transit-time damping for fast modes, the partial restoration of TTD for slow modes due to pitch-angle randomization, and the partial restoration of gyroresonances for Alfvén modes due to the finite lifetime of the modes.

We also observe that our predictions for nonresonant acceleration, i.e., acceleration in a turbulence whose spectrum is composed of structures rather than linear waves, also match the above resonant values, at least in the range of momenta considered, at the upper end of the inertial range. 




An interesting outcome of our test-particle Monte Carlo simulations is to reveal that, as the particle gains energy through its stochastic interactions, it undergoes superdiffusive spatial transport along the parallel direction. This has consequences for the maximal energy of acceleration, in particular whether a given particle can reach the confinement energy r_g∼ R_s, with R_s the size of the source, or not. Assuming R_s∼ L_max, we plot in Fig. <ref> the evolution of the particle energy in units of the confinement energy as a function of the distance traveled in the parallel direction. 



 §.§ Comparison to recent numerical results

Recently, several groups have provided the first PIC numerical simulations of relativistic turbulence <cit.>. The merit of such simulations is to provide a self-consistent description of the nonlinear interaction between the thermal plasma, the nonthermal population and the electromagnetic fields. The main disadvantages of such simulations are, of course, their high numerical cost and their need to resolve (finely enough) the skin depth scale of the plasma. The numerical cost has to be balanced against the number of particles involved, fewer particles per cell meaning a higher numerical noise. Resolving the skin depth scale allows to properly model the small-scale dissipative physics, which is missing in MHD simulations, but it also means that the maximal energy at acceleration is limited by the dynamic range of the simulation, just as the acceleration mechanism may be affected by small-scale physical effects which would be absent on the large scales of astrophysical sources.

Nevertheless, the dynamical range of the recent simulations shown in Refs. <cit.> is so large that one may expect such effects to be absent. Such simulations thus provide a perfect experimental benchmark for our results. These simulations find that acceleration proceeds in two stages. The first is rapid and likely associated to particle energization in current sheets, through, e.g., reconnection-type acceleration, as discussed in these studies. The second is slower and presumably proceeds through stochastic acceleration in a turbulence that mostly obeys the ideal Ohm's law. In particular, Ref. <cit.> derives a diffusion coefficient D_pp∝ p^2/3 in the first stage, but D_pp∝ p^2 in the second; see also Ref. <cit.>. Our discussion of the previous sections agrees well with such a picture. If the small-scale physics can be described in terms of acceleration in parallel electric fields whose power lies on scale k_max, we have found that D_pp∝ t_scatt for t_scatt≪ k_max^-1, while  D_pp∝ t_scatt^q_E for t_scatt≫ k_max^-1 if 0<q_E<1. Consider, for instance, the latter possibility; the scaling of Ref. <cit.> is then  recovered if t_scatt∝ p^2/(3q_E). Alternatively, if t_scatt∝ p^2, describing the limit in which the particle at such energies is sensitive to the electric field structures on scales k_max but not to the large-scale MHD turbulence, then q_E=1/3 would explain the observed scaling.  A detailed study of the statistics of the small-scale parallel electric fields and a careful follow-up of the particle momentum in the corresponding PIC simulations would allow to test such hypotheses.

Regarding the acceleration on larger scales, the main result of Refs. <cit.> is t_acc∼ L_max/⟨δ u^2⟩. This agrees well with the scalings given previously, since it corresponds to D_pp∼ 0.16 ηβ_m^2(r_gk_min)^2 in units of m^2Ω_0^2/k_min, with β_m the relevant mode velocity, for direct comparison to Eqs. (<ref>) and (<ref>). In particular, it broadly agrees with acceleration in Alfvén and/or slow mode turbulence, accounting for wave decay, or with nonresonant acceleration in a generic turbulence. In order to better understand which of these description better applies, one would need to perform a mode decomposition of the turbulence as in Ref. <cit.>, and to track the scaling of t_scatt with energy in these PIC simulations.



§ CONCLUSIONS

This work has studied the physics of particle acceleration in relativistic turbulence, where relativistic means typical 3-velocity fluctuations at the outer scale ⟨δβ^2⟩^1/2∼ 1. For an MHD turbulence described in terms of eigenmodes, fast and slow magnetosonic, and Alfvén, this corresponds to a magnetized, relativistically hot plasma for which β_A is close to unity (and β_s=1/√(3)); hence, σ≳0.1 in terms of magnetization σ=β_A^2/(1-β_A^2).

We have provided detailed analytical estimates of the pitch-angle scattering (D_μμ) and momentum (D_pp) diffusion coefficient in quasilinear, and extended quasilinear theory for acceleration in isotropic fast wave turbulence, as well as in Goldreich-Sridhar-like anisotropic Alfvén and slow wave turbulence. In turn, this allows to provide analytical scalings for the scattering (t_scatt) and acceleration (t_acc) timescales. We have notably included the effect of a finite mode lifetime, as expected for strong turbulence, and discussed in some detail the resonance broadening associated to the partial randomization of the pitch angle of the particle. We have compared these predictions to dedicated Monte Carlo simulations of test-particle acceleration in synthetic wave turbulence. We have paid attention to the notion of local anisotropy, which is inherent to the Goldreich-Sridhar phenomenology of anisotropic turbulence. To do so, we conduct these anisotropic simulations with a low enough turbulence amplitude to guarantee that the eddies are effectively aligned with respect to the mean field on the scale of the particle gyroradius.
Our Monte Carlo simulations naturally account for the partial randomization of the pitch angle and they include, where necessary, the finite lifetime of turbulent modes. The satisfactory agreement that we reach between those analytical scalings and the numerical simulations suggests that our analytical predictions capture the salient effects of particle acceleration in a prescribed wave turbulence.

We notably observe the following effects in a relativistic turbulent setting:

    
  * For an isotropic cascade of fast magnetosonic modes the transit-time damping contribution is progressively reduced as β_F increases, because the resonant value of the longitudinal phase velocity of the wave is then forced into the superluminal regime. Even though resonance broadening effects preserve part of the transit-time damping contribution, most of the acceleration occurs through gyroresonant interactions. The momentum diffusion coefficient accordingly scales as D_pp^F∼β_F^2η  p^2 k_min(r_gk_min)^q-2 for r_gk_min≲1, with q the 1D index of the turbulence spectrum and k_min^-1=L_max/(2) in terms of the outer scale L_max.
    
    
  * For an anisotropic, Goldreich-Sridhar-like cascade of slow mode waves, the transit-time damping contribution provides the dominant  contribution to scattering and acceleration, thanks to the resonance broadening implied by the partial randomization of the pitch-angle. Gyroresonant-like interactions are negligible. The momentum diffusion coefficient accordingly scales as D_pp^S∼β_S^2η p^2 k_min for r_gk_min≲1, up to a logarithmic correction. 
    
    
  * For an anisotropic, Goldreich-Sridhar-like cascade of Alfvén waves, the resonance broadening associated to the finite lifetime of the modes restores gyroresonant-like interactions, that would otherwise be inefficient. The partial randomization of the pitch angle does not play any significant role. The momentum diffusion coefficient scales according to D_pp^A∼β_A^3η p^2 k_min for r_gk_min≲1, up to a logarithmic correction, as for slow modes. The extra power of β_A results from resonance broadening.
    
    
  * At large rigidity, r_gk_min≳1, the particle effectively sees small-scale turbulence, hence D_pp∼β^2η  m^2Ω_0^2k_min^-1∝ p^0 for all three types of modes, with β the relevant phase velocity.
    
    
  * Overall, we find that, in a relativistic setting, if all three types of modes share a similar fraction of the turbulent energy, they give roughly comparable contributions to the acceleration of particles, for momenta not far below the top of the inertial range where r_g∼ k_min^-1. At lower rigidities, r_g≪ k_min^-1, the contribution of fast modes becomes dominant, because of its softer scaling with r_g.

    
  * We have also provided general arguments concerning the acceleration of particles in a turbulence that cannot be described as a bath of linear waves, but rather as a combination of compressive, shearing, vortical and accelerating fluid motions obeying the ideal Ohm's law. We have shown, notably, that if the scattering timescale of particles is such that t_scatt≲ k_min^-1 (as for all three types of modes above at r_g≪ k_min^-1), then the diffusion coefficient D_pp∼⟨δ u^2⟩ p^2 k_min up to a correction factor (t_scattk_min)^q_u-2 that depends weakly on energy. Here, q_u represents the index of the 4-velocity 1D turbulent spectrum, and ⟨δ u^2⟩ represents its amplitude. 
    
    
  * We have also discussed the possible contributions of violations of ideal Ohm's law, showing that they peak on the small scales ∼ k_max^-1 of the turbulent cascade, and we have characterized their magnitude. We notably find that, if the 1D power spectrum index q_E of the parallel electric field component verifies 0<q_E<1 and t_scatt≫ k_max^-1, then such small-scale effects provide a contribution D_pp∼ q^2⟨ E_∥^2⟩ k_max^-1( t_scattk_max)^q_E which may dominate the acceleration at very small rigidities.

    
  * We have compared our results to recent ab initio simulations of turbulence using kinetic particle-in-cell simulations and shown that the above generally agree with the observed results. 

    
Finally, we provide ready-to-use analytical scalings for applications to high-energy astrophysical phenomenology in Sec. <ref>. 



§ ACKNOWLEDGMENTS

We thank the anonymous referees for useful suggestions.

This work has been financially supported by the ANR-14-CE33-0019 MACH project. C.D. and F.C. acknowledge the financial support from the UnivEarthS Labex program  of Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02). C.D. also gratefully acknowledges the support of the Japan Society for the Promotion of Science (JSPS International Research Fellowship at the Astrophysical Big Bang Laboratory, RIKEN). M.L. acknowledges support by the Emergence-2019 SU program and by the National Science Foundation under Grant No. NSF PHY-1748958. This work was granted access to HPC resources of Institut du développement et des ressources en informatique scientifique (IDRIS) and Centre Informatique National de l13Enseignement Supérieur (CINES) under the allocation  A0040410126 made by Grand Equipement National de Calcul Intensif (GENCI). 






§ RESONANCE BROADENING
 IN QUASILINEAR THEORY

The theoretical estimates of the diffusion coefficients, Eqs. (<ref>) and (<ref>) for pitch-angle diffusion, Eqs. (<ref>) and
(<ref>) for momentum diffusion, involve resonance functions  ℛ_k that characterize the interaction between particles and waves. At a formal level, this resonance function derives from the time-integrated Fourier transform of the propagator that connects the position of the particle at different times in the turbulent bath, e.g., Ref. <cit.>. Accounting for gyromotion around the background magnetic field, this resonance function is expressed as

    ℛ_k = 1/2∫_-∞^+∞dτ ^ı(k_∥μ-ω+nΩ)τ.

In standard quasilinear theory, with ω∈ℝ and μ, the initial  pitch-angle cosine of the particle, this resonance becomes a Dirac-function generating the transit-time damping and the infinite harmonic series of gyrosynchrotron resonances, as explained in the text. 

It has long been appreciated, however, that these idealized resonances are actually broadened to some degree by various physical effects. In our case, the two major causes of broadening are the finite lifetime of the turbulent modes and the partial randomization of the pitch angle of the particle. For the sake of the argument, we treat each case separately.

If the linear eigenmodes are assigned a finite lifetime, the mode frequency can be written  ω=ω-ıγ_d|ω|, with γ_d>0. The resonance function then takes on a Breit-Wigner form,

    ℛ_k = 1/γ_d|ω|/(kμ_kμ - ω + nΩ)^2 + γ_d^2ω^2,

whose finite width is directly governed by γ_d|ω|.

The presence of net magnetic field fluctuations implies that the pitch angle of the particle is modulated at all times by a random quantity, which scales as some power of the fluctuation. This consequently broadens the resonance of particles with waves and permits, in particular, transit-time damping to occur in slow mode turbulence over a broad range of particle pitch angles; see Sec. <ref>. This effect has been introduced in the subrelativistic regime in a number of studies, starting with Ref. <cit.>; see, e.g., Refs. <cit.> for recent implementations. Here, we provide a detailed discussion of this effect in the relativistic limit, and emphasize the differences with respect to these previous studies.

If the pitch-angle cosine μ of the particle becomes a random quantity, the ballistic propagator exp[ı(k_∥μ- ω+nΩ)Δ t] can be approximated, to second-order in the cumulant expansion of μ, by

    ⟨^ı(k_∥μ-ω+nΩ)Δ t⟩≃   exp[ı(k_∥⟨μ⟩-ω+nΩ)Δ t
            -1/2k_∥^2⟨Δμ^2⟩Δ t^2] .

Consequently, the resonance function becomes

    ℛ_k = 1/(2 k_∥^2⟨Δμ^2⟩)^1/2exp[-(k_∥⟨μ⟩-ω+nΩ)^2/2k_∥^2⟨Δμ^2⟩],

and it is entirely characterized by ⟨μ⟩ and ⟨Δμ^2⟩.

For strongly magnetized particles, meaning r_gk_min≪1, these moments can be estimated using the conservation laws of the first adiabatic invariant and of the energy, which guarantee that (1-μ^2)/B is a conserved quantity along the trajectory. In this expression, μ represents the pitch-angle cosine as measured relative to the direction of the total magnetic field B. It differs from the standard pitch-angle cosine μ, which is defined relative to the direction of the unperturbed field, by an amount B_⊥^2/2B_0^2∼η/2(1-η).

Let us first consider the case η≪1, as used in Secs. <ref> and <ref> to build an effective model of the anisotropic Goldreich-Sridhar phenomenology, so that μ≃μ. Then, as the particle travels from one coherence cell of the turbulence, on scale k_min^-1, to another, its pitch-angle cosine evolves from μ' to μ”, and to order η,

    μ”^2   ≃μ'^2-(1-μ'^2)(Δ B_∥/B_0+1/2Δ B_⊥^2/B_0^2),

where Δ B_∥= B”_∥- B'_∥ represents the change in the parallel component of the random contribution of the magnetic field, and Δ B_⊥= B”_⊥- B'_⊥ represents the change in its perpendicular component. One important observation is that for slow modes, Δ B_∥/B_0∼[η/(1-η)]^1/2√(2/3), while for Alfvén modes, Δ B_∥/B_0= B_∥=0. For simplicity, we rewrite the second term in the rhs of Eq. (<ref>) as -(1-μ^2)δ in the following, with δ∼η^1/2 for slow modes and δ∼η for Alfvén modes. To go further, it proves convenient to split the pitch-angle domain according to the sign, using here the symmetry μ↔-μ of our estimates. We thus focus here on μ>0. Moreover, we distinguish escaping particles with initial pitch angle cosine larger than δ^1/2 and trapped ones.

The latter case corresponds to particles undergoing mirror reflections. The pitch-angle cosine of such particles becomes randomized with ⟨μ⟩∼δ^1/2 and rms value ⟨Δμ^2⟩^1/2∼δ^1/2, as indicated by Eq. (<ref>). At values of δ not far below unity, this scaling holds for most particles.

For particles inside the loss cone, meaning μ_0≳δ^1/2, the pitch-angle cosine remains confined around the initial value, ⟨μ⟩≃|μ_0|, to within Δμ≃δ/|μ_0|, since Eq. (<ref>) then gives μ”≃μ' - (1-μ'^2)Δ B_∥/(2μ' B_0) to lowest order. 

As a consequence, the resonance is significantly narrower for particles with μ_0≳δ^1/2 than for those with μ_0≲δ^1/2; hence, the pitch-angle averaged diffusion coefficients are strongly dominated by particles in the former range. We have checked this numerically as follows. We first performed a numerical study of a stochastic system that follows the evolution of the pitch-angle cosine through Eq. (<ref>) step by step, accounting for mirroring whenever it occurs, in order to derive accurate estimates of ⟨μ⟩ and ⟨Δμ^2⟩. We then incorporated these estimates in the resonance function and computed the diffusion coefficients in various regimes of interest and compared the obtained values.


In light of these studies, we find that the broadening of the resonance can be modeled, to the level of accuracy that suits the lowest-order expansion Eq. (<ref>), by assuming ⟨μ⟩≃δ^1/2 and ⟨Δμ^2⟩^1/2∼δ^1/2. For comparison, Ref. <cit.> uses ⟨μ⟩=μ_0 and ⟨Δμ^2⟩∼η^1/4∼δ^1/2 for all μ_0, while Ref. <cit.> rather finds Δμ∼η^1/2∼δ. Which regime applies depends on η and μ_0, as explained above. Recall also that, for magnetosonic modes, δ^1/2∼η^1/4; hence, the amount of resonance broadening is significant, even in low-amplitude turbulence.




§ NUMERICAL SIMULATIONS



 §.§ Field prescription


The magnetic field is described as the sum of a static and uniform magnetic field  B_0 and a turbulent component B, expressed as the superposition of N_ k MHD eigenmodes (pure fast, slow or Alfvén modes) with wave vector moduli between
k_min and k_max. k_max is chosen so that k_max^-1<r_g. Given that most of the energy of the turbulence is concentrated at the largest scales, the exact value of k_max is not important and we typically take k_max∼ 100  r_g^-1, while we set the mode density to 128 per decade (which is rather conservative). 



  §.§.§ Isotropic fast mode turbulence


For pure isotropic fast mode turbulence simulations, we draw N_ k waves vectors with equally log-spaced norms and random directions and build the magnetic perturbations at the coordinates (t,r) in the plasma rest frame as

    B(t,r)=∑_i=1^N_ kB_k_i^Fcos(k_i·r-ω_i t+ϕ_i),

where ω_i and B_k_i^F are obtained by injecting the wave vector k_i in the dispersion and polarization relations (<ref>) and (<ref>), and ϕ_i is the random phase. The vectors are normalized according to the chosen spectral scaling and level of turbulence, namely for a Kolmogorov spectrum,

    b_k_i^2=2η/1-ηB_0^2 (k_i/k_min)^-2/3[∑_j=1^N_ k(k_j/k_min)^-2/3]^-1.

In the same fashion, we compute from Eq. (<ref>) the 4-velocity perturbation[For linear modes, there is no distinction between the 3-velocity and the spatial part of the 4-velocity but for finite amplitude modes, we treat δu as a 4-velocity to ensure that |δE| defined by Eq. (<ref>) remains smaller than |δB|.] as 



    u(t,r)=∑_i=1^N_ ku_k_i^Fcos(k_i·r-ω_i t+ϕ_i),

so that the electric field can be built according to the ideal MHD Ohm's law

    E(t,r)=-u(t,r)/√(1+u^2(t,r))×[B_0 +δB(t,r)].




  §.§.§ Anisotropic slow and Alfvén mode turbulence


For slow and Alfvén waves, we enforce approximately the Goldreich-Sridhar scaling. The difficulty resides in that we need k_i to compute B_k_i^S (B_k_i^A), but the direction of k_i depends on the direction of the local magnetic field that we are trying to build, so a recursive procedure would be required. For simplicity, we assume that the direction of the local field can be approximated to that of B_0. This approximation is reasonable if the perturbation is small enough at the scales of interest, that is, B/B≲ k_∥/k∼(k_∥/k_min)^-1/2, which is increasingly constraining as we move to smaller scales. We therefore limit the dynamic range of these simulations and only consider r_gk_min≳ 10^-2. For this limiting case and at scales close to the resonance condition k_∥∼ r_g^-1, we thus should have B/B≲0.1 or η≲ 10^-2, which is the value that we adopt.

Depending on whether we consider undamped or damped waves, we use one of the following random field construction method.

For undamped waves, the procedure to construct the perturbed fields is essentially the same as the one described above for fast mode turbulence. N_ k wave vectors are drawn with evenly log-spaced perpendicular components between √(2/3)k_min and k_max along random directions in the plane perpendicular to B_0 (in line with the remarks of the precedent paragraph), while the parallel components are defined as in Ref. <cit.>,

    k_i∥=±√(2)/2k_i⊥^2/3k_min⊥^1/3.

For slow modes, we then take in combination with Eqs. (<ref>), (<ref>), and (<ref>),

    B(t,r)    =∑_i=1^N_ kB_k_i^Scos(k_i·r-ω_i t+ϕ_i),
    u(t,r)    =∑_i=1^N_ ku_k_i^Scos(k_i·r-ω_i t+ϕ_i),

with amplitudes scaling as

    b_k_i^2=2η/1-ηB_0^2 (k_i⊥/k_min⊥)^-2/3[∑_j=1^N_ k(k_j⊥/k_min⊥)^-2/3]^-1,

and construct the electric field as previously explained. A similar method is employed for Alfvén modes. 

For simulations with modes of finite time correlation, along the same lines as Ref. <cit.>, we add, compared to the previously described procedure, some pulsation spreading around the solution  ω_k_i of the linear dispersion relation (<ref>),

    B(t,r)    =∑_i=1^N_ k∑_j=1^N_ωB_ij^Scos(k_i·r-ω_ij t+ϕ_ij),
    u(t,r)    =∑_i=1^N_ k∑_j=1^N_ωu_ij^Scos(k_i·r-ω_ij t+ϕ_ij),

using the polarization given by Eqs. (<ref>) and (<ref>) for k_i, the corresponding ω_k_i and the magnitude

    b_ij^2=   2η/1-ηB_0^2 Δ_i 𝒳_ij(k_i⊥/k_min⊥)^-2/3
       ×[∑_mnΔ_m 𝒳_mn(k_m⊥/k_min⊥)^-2/3]^-1,

where Δ_i is the constant logarithmic spacing of ω_ij at fixed i and 𝒳_ij is the function describing the decorrelation of modes away from ω_k_i and whose exact form is model dependent. 
To make contact with our theoretical calculations, we choose

    𝒳_ij=γ_d |ω_k_i|/|ω_ij|/γ_d^2 ω_k_i^2+(ω_ij-ω_k_i)^2,

which is the Fourier transform of t↦^ıω_k_it-γ_d|ω_k_i|t, and corresponds to the "Nonlinear Anisotropic Dynamical Turbulence" model of <cit.>. The pulsation cutoff is chosen to be much larger than the linear pulsation (at least 1000  |ω_k_i|) and the mode density to at least eight modes per decade. 



 §.§ Measurements


For a given physical setup and particle rigidity, a thousand of ultrarelativistic particles are tracked, each injected with a random initial pitch angle in a different turbulence realization. The Bulirsch-Stoer integrator <cit.> is used to evolve the positions r and 4-velocities Γβ in the plasma rest frame, according to the equations of motion 

    dr/dt=β,
    
    
         dΓβ/dt=q/m(E + β×B),

where q and m are the electric charge and the rest mass of the particle and B=B_0+B is the total magnetic field, constructed according to the procedure presented above.

For most simulations, the initial momentum distribution of the particles is a Dirac-delta function in the simulation frame so that all particles have the same initial r_gk_min in that frame. For fast mode turbulence simulations with β_F≥ 0.7 (and η=0.3), we also carried out simulations such that the gyroradius distribution is a delta function in the local rest frame of the fluid (in which the motional electric field is initially zero) in order to be closer to the situation described by our analytical calculations. If r'_g,0=p'_0/(qB') is the value of the gyroradius to be probed in that frame, then in the simulation frame, we initialize the particles with a momentum

    p_0=p'_0/γ_u(1-β_u μ_u)=q (B_⊥^2/γ_u^2+B_∥^2)^1/2/γ_u(1-β_u μ_u)r'_g,0,

where β_u (γ_u) denotes the velocity (the Lorentz factor) of the fluid at the initial position of the particle, μ_u is the cosine of the angle its initial momentum makes with the flow and B_⊥ (B_∥) is the perpendicular (parallel) component of the local magnetic field. For these simulations, we report the measurements of the diffusion coefficients (Figs. <ref> and <ref>) at an effective gyroradius that is defined as the mean momentum, after any transient initial energy kick, divided by qB_0.



  §.§.§ Pitch-angle diffusion coefficient

The pitch-angle diffusion coefficient is estimated by computing the mean square pitch-angle displacement ⟨[Δμ(t)]^2⟩/2t, where ⟨…⟩ denotes the average over the particles, i.e., over the pitch angles and the turbulence realizations. We fit ⟨[Δμ(t)]^2⟩ using least squares to a linear function at early times (but still larger than r_g to omit the ballistic regime) and identify D_μμ to one-half of the slope of the fitting function. One limitation is that irrespective of the energy dependence of the diffusion coefficients, ⟨[Δμ(t)]^2⟩/2t eventually displays a subdiffusive behaviour as t increases due to the bounded nature of Δμ(t) and a well-defined plateau cannot always be observed. 

These results were cross-checked against another estimation relying on the pitch-angle correlation function ⟨μ(t) μ(0) ⟩, where the data are fitted against a decaying exponential model ∝exp(-t/t_scatt) and D_μμ is estimated as 1/t_scatt. Most results were found to differ by a factor of a few[In details, when derived from the correlation function, D_μμ is found to be ≈ 3 times larger for fast modes,  ≈ 4 times larger for Alfvén modes, ≈ equal for damped slow modes, and ≈ 1.7 times smaller for undamped slow modes.] (see, for instance, Fig. <ref>), although for a few simulations (namely that of fast modes with β_F≳ 0.7, r_gk_min≳ 1), the data are not well fitted by an exponential model. We also note that this method requires longer computation times than when using the running displacement. 





  §.§.§ Momentum diffusion coefficient and acceleration timescale

In a similar fashion, the pitch-angle averaged momentum diffusion coefficient is evaluated from the dispersion of the energy distribution with D_pp set to the slope of the linear fit of σ_p^2(t)/2 where σ_p is the standard deviation of the distribution of p at early times. For subrelativistic simulations, this is equivalent to using the mean square displacement around the initial value ⟨[Δ p(t)]^2⟩/2t. For relativistic setups, particles are subject to a quick initial increase of energy due to some first order Fermi acceleration; hence, using ⟨[Δ p(t)]^2⟩/2t would lead to overestimating D_pp. This initial acceleration moreover implies that we do not measure the diffusion coefficient for the initial r_gk_min but rather a slightly larger value (approximately, 1.4 larger for β_A=0.7 and 1.6 larger for β_A=0.9, when the particle momentum distribution is initialized in the simulation frame).
Finally, we note that the energy distribution is not always well described by a Gaussian, e.g., Refs. <cit.>. In particular, for undamped slow modes turbulence in which there are two different populations of particles (those sensitive to resonance broadening and those which are not), the distributions display large wings (see, for instance, Fig. <ref>), and the above described procedure leads to overestimating the acceleration efficiency.

In light of this last remark, the pitch-angle cosine averaged acceleration timescale is not derived from D_pp but defined as the time when the average Lorentz factor has doubled, ⟨Γ (t_acc) ⟩=2 Γ_0. Both approaches lead to similar results for most simulations anyway[In details, for r_gk_min<1, 2t_accD_pp/p^2 takes the following values. For fast modes, ≈ 0.4 at β_A=10^-2, ≈ 0.5 at β_A=10^-1, and ≈ 2 at β_A≥ 0.7.  For undamped and damped Alfvén modes, ≈ 0.7 and ≈ 0.5, respectively. For undamped and damped slow modes, ≈ 4 and ≈ 0.4, respectively. For r_gk_min≈ 1, 1 ≲ t_accD_pp/2p^2 ≲ 1.5, regardless of the turbulence type.]. 








§ POWER SPECTRUM OF ELECTRIC FIELD FLUCTUATIONS
 IN WAVE TURBULENCE

In a pair plasma, inertia effects appear in the Ohm's law and generate the following extra electric field component <cit.>,

    E^×≃∂/∂ t(w/2e^2n^2 j),

emphasizing that E^× is here measured in the plasma rest frame; w represents the enthalpy density, n is the number density, and j is the total current density. In Fourier space, relating j to B to lowest order in ω/k, we obtain for this extra electric field component,

    E^×_k = κω/4ω_p^2 k×B_k,

where κ = [1+T/mΓ̂/Γ̂-1],  Γ̂ denotes the adiabatic index and, T is the temperature. The scaling E^×∝ω k  B_k confirms that the parallel electric field power is maximum on the smallest length scales of the turbulent cascade.

Although each component of this small-scale electric field is orthogonal to the corresponding magnetic field component in Fourier space, this is no longer the case once the sum over k is performed. More specifically, ⟨E^×·B_0⟩ = 0 and ⟨E^×·B⟩ = 0 for all modes, but ⟨(E^×·B_0)^2⟩≠ 0 and ⟨(E^×· B)^2⟩≠ 0 for Alfvén polarization, while ⟨(E^×·B_0)^2⟩ = 0 but ⟨(E^×· B)^2⟩≠ 0 for magnetosonic modes.

Using the polarization of each mode and the corresponding spectrum described earlier, one can calculate the rms parallel electric field E_∥ = ⟨ ( E^×·B)^2⟩^1/2/B, recalling that B = B_0+B. 

For an isotropic cascade of fast modes, we derive

    E^F_∥≃κ/4√(30) β_F η k_max^5/3k_min^1/3/ω_p^2 B.


For a Goldreich-Sridhar cascade of slow modes,

    E^S_∥≃κ/16 β_S η k_max^4/3k_min^2/3/ω_p^2 B.


Finally, for a Goldreich-Sridhar cascade of Alfvén modes, we obtain

    E^A_∥≃   κ/8 β_A k_maxk_min/ω_p^2
       × B [√(2/3)η+√(η(1-η))(k_max/k_min)^1/3].

The second term in brackets corresponds to the component parallel to B_0. 

In the case of electron-ion plasmas, the above corrections are completed by the Biermann battery term, which provides the additional electric field

    E^×_k=ık p_e_k/n

in terms of the electron pressure fluctuation p_e_k. Using p_e_k ∼Γ̂T_e/m_pk·β_k/ω, as applies to linear magnetosonic modes, we find here as well that the power of E_∥ peaks at small scales. Its contribution is dominated by the product E^×_k·B_0, which gives

    E^F_∥≃Γ̂/√(6)T_e/m_p k_max^2/3k_min^1/3/ω_ci √(η(1-η)) B

for fast modes and 

    E^S_∥≃   Γ̂ β_S^2/2√(2)β_s^2T_e/m_p k_max^1/3k_min^2/3/ω_ci √(η(1-η)) B
       ×[2√(2)+√(η/1-η)(k_max/k_min)^1/3]

for slow modes, in terms of ω_ci = q B/(m_p c), the proton cyclotron frequency.

The various dependencies on k_max of the above estimates provide the scale dependence of the corresponding electric field; i.e., if the spectral index of the (one-dimensional) power spectrum of E_∥ is q_E, meaning k^2 ⟨| E_∥_k|^2⟩∝ k^-q_E, then E_∥∝  k_max^(1-q_E)/2, assuming q_E<1. In detail, Eq. (<ref>) for fast modes leads to q_E=-7/3, 
Eq. (<ref>) for slow modes leads to q_E=-5/3, and
Eq. (<ref>) for Alfvén modes leads to q_E=-1. Equivalently, the parallel electric field rms strength on a given scale k scales as . E_∥|_k ∝ k^(1-q_E)/2, to be contrasted with a typical . E|_k ∝ k^-1/3 for the MHD component. Consequently, the parallel electric field energy may represent a non-negligible contribution to the total electromagnetic energy density, but in all cases studied above, its power is strongly suppressed on the largest scales.







