Stochastic particle acceleration in magnetized turbulence has emerged as a key process in high-energy astrophysics, as it is likely at play in a wide variety of sources and in a diverse set of physical conditions, from the solar atmosphere [1][2][3] up to more exotic objects, e.g., blazars, gamma-ray bursts, and other relativistic outflows [4][5][6][7][8][9][10]. In phenomenological applications, particle acceleration is commonly characterized by a diffusion coefficient in momentum space, whose magnitude and scaling control the evolution of the particle distribution function in a Fokker-Planck description. Most often, this diffusion coefficient is calculated in the framework of quasilinear theory, which ascribes particle energization to resonant interactions with linear eigenmodes of the plasma, e.g., [11,12]. Quasilinear theory, however, is a perturbative description restricted to the regime of small-amplitude turbulence, i.e., δB ≪ B, with δB the turbulent rms magnetic field fluctuation on the integral scale l c of the turbulence cascade, and B the mean field strength. To the contrary, the turbulence can be regarded as of large amplitude in most astrophysical settings. Furthermore, stochastic acceleration is expected to be rather slow in the small-amplitude limit, since the acceleration timescale, i.e., the time needed on average to double a particle energy, scales in proportion to ðδB=BÞ -2 .
The regime of strong turbulence is thus of broader applicability and interest, if not for practical purposes. Recently, kinetic particle-in-cell (PIC) simulations have offered ab initio numerical probes of that regime, in the relativistic limit where the Alfvén velocity v A ∼ c [13][14][15][16][17][18][19][20][21][22][23]. Quite interestingly, those experiments have measured a momentum diffusion coefficient D pp ∼ 0.1hδufoot_1 ip 2 c=l c , here written in terms of the variance of four-velocity fluctuations hδu 2 i and particle momentum p, whose p 2 dependence is suggestive of a nonresonant form of acceleration.
Nonresonant acceleration can be described in various ways. In the original Fermi picture [24], acceleration proceeds through repeated encounters with discrete, pointlike moving magnetized structures. The particle then draws energy from the motional electric field carried by those highly conducting plasma elements; accordingly, in the reference frame of the scattering center, the interaction is assumed to be purely elastic and described as pitch-angle scattering. If the mean free path to interaction t mfp does not depend on energy, and if the particle distribution function is assumed isotropic before each interaction, then the energization process can be indeed described as a random walk in momentum space with diffusion coefficient D pp ∝ p 2 =t mfp [24][25][26].
Generalizing this picture to a turbulent flow, the acceleration of a particle can now be related in general terms to the time variation of the velocity flow that the particle encounters along its journey, i.e., to its interaction with the time-dependent sheared and compressive parts of the drift velocity field v E ¼ E × B=B 2 (in units of c), which characterizes the velocity of magnetic field lines in the magnetohydrodynamics (MHD) approximation [27][28][29][30][31][32][33][34]. More specifically, for particles whose gyroradius r g is much smaller than the integral scale l c , nonresonant acceleration derives from three main force terms expressed in terms of u E , the four-velocity generalizing v E : an inertial contribution characterizing the time dependence of u E , its shear as projected along the mean magnetic field direction, and its compression in the plane transverse to the mean magnetic field [28]. Those force terms and their overall combination will be detailed further below.
The objective of the present paper is to test this nonresonant acceleration model through a direct comparison to numerical simulations of strong turbulence. To do so, we rely on dedicated large-scale PIC simulations of forced turbulence in 2D and 3D, as well as of decaying turbulence in 2D; we also make use of an incompressible MHD simulation, borrowed from the Johns Hopkins University turbulence database. We frame the discussion as follows. In Sec. II, we detail the theoretical model and explain the metric that we use to test it on numerical simulations. In Sec. III, we then present the numerical simulations and carry out our test. This comparison, which proves quite satisfactory, is summarized and discussed in Sec. IV. Finally, we present some similar tests of our model against simulations of test-particle acceleration in synthetic wave turbulence in the Appendix, which are known to be well described by quasilinear theory, even at relatively large turbulence amplitude δB=B ∼ 1.
Let us first emphasize some differences between the more conventional approach based on quasilinear theory and the nonresonant acceleration scenario that we study here; see also [28] for a detailed discussion of those issues. In its simplest formulation, quasilinear theory describes the turbulence as a linear superposition of eigenmodes (fast and slow MHD magnetosonic waves, and Alfvén waves in the MHD approximation) of infinite spatial extent and well defined frequency [11]. Nonlinear extensions include resonance-broadening effects or stochastic corrections to the particle orbits, see [35] for a recent general discussion, and references therein. In the framework of modern, anisotropic MHD turbulence theories [36,37], resonance broadening effects associated with wave damping are actually necessary to restore a part of wave-particle interactions, which would otherwise disappear at small r g =l c [38], with the exception of fast MHD modes, which preserve resonant interactions [39]. Such resonant interactions would lead to a scaling D pp ∝ p q , with q the (absolute value of the) index of the one-dimensional power spectrum of magnetic fluctuations contained in fast modes. This disagrees with the scaling observed in recent PIC simulations, if q < 2 as expected and observed [40,41]. 1By contrast, resonance-broadened interactions yield D pp ∝ p 2 up to logarithmic corrections [35,38]. This scaling emerges because, once the resonance is broadened, all scales above the gyroradius give equal contributions to the diffusion coefficient. In the context of strong turbulence, this contradicts the quasilinear assumption, because the wave-particle interaction timescale (∼l c =c) becomes larger than the timescale over which the magnetic field strength and direction have changed by the order of unity. Moreover, in such turbulence, the eigenmodes that contribute on all scales up to l c are aligned with respect to different mean field directions, since this orientation changes with scale [43]. The polarization of those misaligned eddies, as seen by the particle, thus does not correspond to that of linear eigenmodes, further quelling any resonance. From its point of view, meaning as seen on a scale r g , the particle rather experiences those various modes as random velocity structures distributed on all scales up to l c , and it interacts with them in a nonresonant way.
In that sense, nonresonant acceleration as we discuss it is not antagonistic to the presence of waves or wave packets in a turbulent bath; it assumes, however, that sharp waveparticle resonances are absent, or at least subdominant. The results that we present further below will support that point of view.
Although the nonresonant acceleration process that we test here can be regarded as a generalization of the original Fermi scenario to a continuous turbulent flow, there are noteworthy differences. In particular, the particle does not gain energy depending on whether its interaction with the velocity structure is head-on or tail-on in the present case. It rather gains or loses energy depending on the sign of the variation of the velocity flow while inside the structure [28,34]. For instance, in a region that undergoes compression, energy gain is positive, while it is negative if the flow undergoes dilation, as discussed further below. That difference can be traced back to the underlying assumptions: the original model of E. Fermi depicts discrete, pointlike interactions in a fixed laboratory frame, while our nonresonant model rather discusses how particles interact with modes of extent ≳r g in a comoving frame. The two pictures do not contradict each other, however. To see this, consider a moving magnetic mirror (type-A interaction in the Fermi model): in the comoving frame, the particle will experience in (comoving) time a compression or a dilation depending on whether the mirror is moving towards (head-on) or away from (tail-on) the particle, leading, respectively, to energy gain or loss.
Nonresonant acceleration is best described by making use of the formalism developed in [28,34], which follows the momentum of the particle in a sequence of frames in which the electric field vanishes. One assumption of that model is that this reference frame exists, meaning that the Lorentz invariant quantity E 2 -B 2 is everywhere negative (or at least, over most of space). This is in particular satisfied in the ideal MHD approximation, in which case E ¼ -v p × B (henceforth all velocities are written in units of c), with v p the plasma bulk velocity. We assume here that this MHD regime applies. In the following, the frame in which E vanishes is denoted R = E ; with respect to the laboratory frame, it moves with the velocity v E defined earlier. Quantities expressed in R = E will be annotated with a prime; while this distinction does not matter in a nonrelativistic setting, it matters here, because we will compare our results to PIC simulations of relativistic turbulence. We emphasize, however, that the present treatment is general, and that it can be applied equally well to the subrelativistic or to the relativistic regime.
The advantage of this approach is that it allows one to substitute the electric field for the gradients of v E in the expression of forces acting on the particles. Energy gains and losses are thus directly connected to the statistics of the velocity structures. At this stage, this can be seen as plain rewriting, no information has been lost. For practical purposes, though, this model can be simplified further by noting that turbulent modes of wavelengths significantly smaller than the gyroradius of the particle contribute little, if at all, to the process of energization. Retaining only the contribution of modes of scales ≳r g , one can approximate the gradients as their average over a gyrating orbit of the particle around the field line. Note that we do not follow unperturbed trajectories around a mean field line, as in quasilinear theory, but the exact trajectory around the perturbed field line.
In that approximation, the energization of the particle can be written in terms of three force terms only, which are proportional to the quantities Θ k , Θ ⊥ and a E • b. These have the following definitions and interpretations:
denotes the projection of the shear of the (field line) four-velocity field
along the direction of the mean field direction b (b a unit vector). Here α; β ∈ f0; …; 3g are spacetime indices, b α ¼ f0; bg and u E α ¼ fγ E ; u E g. The mean-field direction is interpreted as the sum of the coherent magnetic field B 0 and of all modes on scales larger than the gyroradius, i.e., b ¼ B= B with
where δB l represents the perturbation coarse grained on scale l ¼ r g . It should be understood as a fluctuation seen at that position, after filtering out scales smaller than l. Consequently, the mean field is here a function of both particle position and momentum.
The force term associated with Θ k leads to energy gain if negative, to energy loss if positive. It can be related to the projection on E of the drift velocity (in a guiding center formulation) due to the curvature of the magnetic field line [28]. This term thus describes acceleration through the curvature term, or in the phrasing of the original Fermi mechanism, a Fermi type-B interaction which energizes the particle as in a slingshot.
The perpendicular shear term Θ ⊥ is defined as
with η αβ the Minkowski metric. It thus corresponds to the shear of the field line velocity field in the plane transverse to the magnetic field line. This term leads to energy gain or loss in much the same way as a particle gains or loses energy through the compression or expansion of a collisional plasma. Here, the particle is tied to the mean field through its orbit, thus transverse compression energizes the particle, while transverse expansion draws energy from it. This term also depicts the effect of a magnetic mirror, or Fermi type-A interaction, and as such it can be seen as a form of betatron acceleration [28].
Finally, the acceleration term a E • b depicts the influence of the effective gravity associated with the noninertial nature of R = E as the particle travels along the field line. More explicitly, a E is the (Lagrangian) three-acceleration of u E , viz.
With these definitions, the theoretical acceleration rate of a particle of mass m, as expressed in R = E , can be written as
with γ 0 ¼ ϵ 0 =mc 2 (respectively, ϵ 0 ) represents the Lorentz factor (respectively, the energy) of the particle in R = E ; u 0 k ¼ p 0 • b=mc (respectively, p 0 ) its four-velocity projected along the mean magnetic field, in units of c (respectively, its three-momentum); correspondingly,
denotes the perpendicular component, with u 0 ¼ p 0 =mc. Finally, dτ ¼ dt 0 =γ 0 represents an element of proper time.
We test the above nonresonant acceleration model against numerical simulations of particle acceleration in turbulence in the following way. For a given particle in a given simulation, we measure the history of its Lorentz factor, which we write γ 0 obs ðtÞ, as a function of t the time measured in the simulation frame. In parallel, we reconstruct a theoretical history γ 0 th ðtÞ defined as follows in terms of the quantities that appear in Eq. ( 5):
where t 0 represents some initial time. To do so, we extract from the numerical simulation, at each point of the trajectory, the various quantities that enter this equation, namely u E β and its gradients ∂ α u E β . We then reconstruct the fields a E • b, Θ k and Θ ⊥ and use them to predict, at each time step of the trajectory, the change in γ 0 . This reconstruction is affected by several effects. For one, the above model is an approximation obtained in the limit r g ≪ l c , where l c denotes the coherence scale of the turbulent power spectrum, i.e., the length scale on which most of the turbulent power lies. We therefore expect some effects of order r g =l c to alter the model predictions. Because numerical simulations are restricted in their dynamic range (an effective rigidity ρ ¼ 2πr g =l c ∼ 0.03-0.1 is close to what that can be currently achieved at best), such effects can be significant, especially in regions of low magnetic field strength, as r g can then take large values relative to its average at a given energy. Likewise, particles can experience substantial acceleration over a period of time, which also leads to an increase in r g , and hence to a loss of accuracy of the model predictions.
Furthermore, the model assumes that the fields u E , B 0 and their gradients are coarse-grained quantities, meaning that sub-Larmor scales have been filtered out. Such a procedure is too costly to be put in place in PIC simulations, as the Larmor scale changes from particle to particle, and even from time step to time step, since the energy of a particle is not constant. We thus use the actual fields and gradients, as measured in the simulation on the scale of the numerical grid, and discard any filtering. This introduces high-frequency noise in the reconstruction of γ 0 th , associated with small-scale effects. To test how this may affect our comparison, we have also performed a reconstruction of the trajectory including time filtering, which allows us to smooth the time profile of the history of γ 0 obs ðtÞ on timescales ≃r g =c. More explicitly, we smooth the fields u E , B 0 and their gradients that the particle encounters on its trajectory before performing the reconstruction. While the reconstructed trajectory differs from that obtained in the absence of this time filtering, the overall result remains similar.
Finally, the comparison with the model is made further complicated by nonideal MHD effects. Kinetic simulations have demonstrated that particles of the thermal pool initially gain energy through nonideal parallel electric field components in reconnection layers, then start to probe the large-scale turbulence once their gyroradius exceeds the typical scale of those layers [14,17]. The contribution of nonideal parallel electric fields is observed to weaken with increasing particle energy, in general agreement with the idea that on large scales, the physics tends toward the ideal MHD regime, as assumed by the nonresonant acceleration model. Those nonideal electric fields, however, can affect the energy gain process on the simulation scales, and thus perturb the comparison. For this reason, on the one hand, we have evaluated in kinetic simulations the fraction of points along the particle trajectories where the ideal MHD condition is violated (E > B), which turns out not to exceed 0.05% in 3D geometry. On the other hand, we have compared the model with the dynamics of test particles propagated in the fields extracted from a 3D ideal MHD simulation.
From the above considerations, we do not expect an exact match between γ 0 obs and γ 0 th . To test the model, we thus calculate, for each particle i, a Pearson correlation coefficient r i , r i ≡ cov½γ 0 obs ðiÞ; γ 0 th ðiÞ cov½γ 0 obs ðiÞ; γ 0 obs ðiÞ 1=2 cov½γ 0 th ðiÞ; γ 0 th ðiÞ 1=2 ; ð7Þ
where cov½AðiÞ; BðiÞ represents the covariance of the histories of the quantities A and B over the trajectory of particle i.
We then collect, for a large number of particles, the sample of correlation coefficients and establish a probability density. We thus seek to see what fraction of those correlation coefficients lies within the vicinity of þ1, that value denoting perfect correlation, hence perfect reconstruction.
As will be detailed in the following, sub-Larmor effects (or nonideal electric fields on small length scales) can lead to a sharp departure in the history of γ 0 obs , and to a different departure in γ 0 th . This departure is sharp, because smallscale effects are associated with short timescales. For some particles, the two histories can thus reveal strong correlation before and after this sudden event, but the global trajectory itself will show a lesser degree of correlation. To bypass such effects, we have performed several reconstructions, which differ in the duration over which the trajectories are examined. In the following, we present reconstructions for intermediate timescales and for the entire duration of the simulation. The duration of the interval over which we follow the trajectories is written Δt. Additionally, we perform the correlation test for a sample of particles for which the energy varies by a significant factor, in order to test Eq. ( 5), treating on an equal footing energy gains and losses. To select the particles, we adopt a threshold g min and consider those trajectories, or chunks of trajectories, that satisfy Δγ 0 =γ 0 ≥ g min with Δγ 0 ¼ maxðγ 0 Þminðγ 0 Þ over the interval of duration Δt.
In the following section, we describe the numerical simulations, the results of the reconstruction, including some details specific to each. In the Appendix, we conduct a similar experiment on simulations that follow test particles in a synthetic turbulence, meaning a turbulence that is constructed from a sum of noninteracting linear eigenmodes (Alfvén, fast or slow magnetosonic modes) of the plasma, following the study of Ref. [35]. The interest of that experiment is that the physics of particle acceleration in such turbulence is relatively well understood, as it follows the predictions of quasilinear theory, and that part of it (transit-time damping acceleration related to magnetic mirroring effects) can be captured by the above model. It can therefore be used to gauge the amount of information contained in the probability density of correlation coefficients that we reconstruct.
In this section, we report on the comparison between the histories predicted by the model and those observed in several numerical experiments: (i) a 2D decaying turbulence 10 000foot_2 PIC simulation; (ii) a 2D forced turbulence 10 000 2 PIC simulation; (iii) a 3D forced turbulence 1080 3 PIC simulation; and (iv) a 3D forced turbulence 1024 3 MHD simulation.
All PIC simulations assume a pair plasma composition. They have been conducted using the finite-difference time domain, relativistic PIC CALDER code [44], to which a turbulence stirring module has recently been added. The MHD simulation is that made available for public use on the Johns Hopkins Turbulence database 2 [45,46]. The purpose of using different physical parameters and simulation frameworks is to test the capability of our model to describe the acceleration process under various turbulent regimes.
We initialize a 2D decaying turbulence PIC simulation with the following characteristics: domain size N x × N y ¼ 10 000 2 cells, corresponding to physical size L x × L y ¼ 1000 2 c 2 =ω 2 p , integrated over a time T ¼ 5000 ω -1 p . Here, ω p ¼ ð4πn AE e 2 =mÞ 1=2 , with n AE the initial (uniform) proper density of positrons/electrons and e the elementary charge, represents the nonrelativistic plasma frequency of one species, so 1= ffiffi ffi 2 p of the total plasma frequency. The cell size is δx ¼ δy ¼ 0.1 c=ω p , and the time step is δt ¼ 0.099 ω -1 p . The plasma is initialized with 10 particles per species per cell, sampled from a Maxwell-Jüttner distribution function with a temperature T 0 ¼ 1 MeV. The particle count ensures satisfactory statistics when making use of additional filtering as discussed below. Note also that, due to their large gyroradius, the high-energy particles studied here are less sensitive than thermal particles to small-scale field fluctuations. Periodic boundary conditions are used for both fields and particles in all directions.
Turbulence is excited as in Refs. [14,17,18], i.e., a decaying turbulence initialized as a sum of plane waves. Here, we use 24 wave numbers, with mean wave number hki ¼ ð2π=L x Þ × 2.9, corresponding to a stirring scale l c ≃ L x =2.9 ≃ 350 c=ω p . The square root of the average of the squared wave numbers give a similar estimate,
Here, we will not distinguish the stirring scale from the coherence (or integral) scale; various definitions exist for the latter, which give values within a factor of the order of unity of l c . The modes initialized at time 0 excite δB x and δB y fluctuations which are left to evolve freely with the plasma at time t > 0. The mean field lies in the out-of-plane direction (along z).
We define the magnetization parameter as
where hB 2 i is the mean-squared (coherent or turbulent) magnetic field and w the plasma enthalpy density. For reference, w ≃ 8ðn þ þ n -Þmc 2 for an electron-positron pair plasma of 1 MeV temperature. The magnetization associated with the mean-field component is σ 0 ¼ 1.6. Since δB=B 0 ≃ 2.8, the turbulent magnetization is σ δB ≃ 13. In Fig. 1, we show the power spectrum of magnetic fluctuations as measured in this 2D decaying turbulence PIC simulation. It reveals a general scaling close to k -5=3 at large scales, followed by a steeper spectrum characteristic of the dissipation range. This shape generally matches that observed in previous PIC simulations of decaying turbulence [14,17,18].
In that figure, we indicate by a dashed line the scale corresponding to the inverse gyroradius of particles with the initial-meaning, at the time t ¼ 1500 ω -1 p ∼ 4l c =c at which we initiate the test-Lorentz factor γ ∼ 50, which we follow in order to compare the model to the data using the method described earlier. 3 We recall that this model assumes r g ≪ l c , hence r g cannot be made arbitrarily larger. However, it cannot be made arbitrarily small either, otherwise the particle gyroradius will lie out of the range of the inertial (nondissipative) spectrum. On small spatial scales, corresponding to gyroradii of particles with energies in the thermal part of the spectrum-γ ∼ 10-particle energization is furthermore mostly controlled by parallel electric fields, as recalled above [14,17]. We thus conclude that Lorentz factors in the range ∼20-60 provide a reasonable compromise to test the theoretical model of nonresonant acceleration. In the present case, the effective rigidity 2πr g =l c of particles with Lorentz factor γ ∼ 50 is of the order of 0.1, in the range anticipated earlier.
In Fig. 2, we show the energy spectrum of particles in this 2D decaying turbulence simulation, at time t ≃ 5l c =c. The power-law tail, with index s ≃ -2, extends from Lorentz factors γ ∼ 20 up to γ ∼ 10 3 , at which point the gyroradius of particles becomes comparable to the maximal scale of the turbulent cascade, implying less efficient acceleration at larger energies.
We now turn to the comparison between the model predictions and the PIC simulation. The trajectories of a sample/subset of particles are recorded from t ¼ 1500 ω -1 p ∼ 4l c =c up to 5000 ω -1 p ∼ 14l c =c. The initial time ensures that the turbulence has had time to cascade down to small scales by the time the test starts. In the PIC simulations, both time and space derivatives are calculated using simple first-order differences, i.e., from one cell to the next (or one step to the next for time). The time derivatives are smoothed through 16 repeated applications of binomial filtering. The spatial derivatives are computed from fields that also underwent 16 successive applications of binomial filtering. This helps eliminate shot noise on the scale of the mesh (here, ∼0.1 c=ω p ) that would otherwise pollute the reconstruction of derivatives which, as discussed before, are meant to be calculated on scales significantly larger than the grid size.
The test has been carried out for two typical durations: Δt ≃ 1l c =c, see Fig. 3 and Δt ≃ 10l c =c, see Fig. 4. Those figures present the probability density function (PDF) of the correlation coefficients r i , as defined in Eq. (7).
To construct the histogram shown in Fig. 3, we have selected at random, for each test particle, chunks of trajectories in which the energy of the particle changes by an amount at least equal to unity, i.e., Δγ 0 =γ 0 ≥ 1 with Δγ 0 ¼ maxðγ 0 Þminðγ 0 Þ. We typically use 10 4 test particles to construct such a histogram; each test particle history, extending over ≃10l c , is sampled at most 10 times to obtain a chunk of extent 1l c =c. In Fig. 4, we integrate over the whole trajectory of each test particle, provided jΔγ 0 =γ 0 j ≥ 2.
Figure 3 indicates a genuinely positive degree of correlation for the contribution of the Θ k force terms, and similarly when all contributions are summed together as in Eq. ( 5). More specifically, to plot the probability density of the correlation coefficients for one force contribution, we use Eq. ( 5) but set the contributions of the other two terms to zero. This figure suggests that neither the force term Θ ⊥ nor a E • b appear to contribute strongly to the evolution of the particle energy. The dominance of Θ k is a common trait to our PIC simulations, which will also hold in 3D as discussed further on.
The trend observed in Fig. 4 is similar. The level of noise is larger in that figure, because we can select only one trajectory for each test particle instead of a number of distinct time intervals, and because our stronger constraint on the amount of energy variation within the interval limits further the number of test particles that are selected for the test.
We note that the above figures and results are relatively insensitive to the choice of the threshold of energy variation jΔγ 0 =γ 0 j, as we have verified. It is also somewhat insensitive to the duration of the interval that we consider. The latter must be large enough, obviously, to accommodate a large number of gyroperiods, since the model considers only contributions from scales larger than r g .
We now analyze a 2D driven turbulence PIC simulation with characteristics similar to that for the decaying turbulence scenario: domain size N x × N y ¼ 10 000 2 cells, corresponding to physical size L x × L y ¼ 1000 2 c 2 =ω 2 p , integrated over time T ¼ 5000 ω -1 p ; the cell and step size are, as before, δx ¼ δy ¼ 0.1 c=ω p and δt ¼ 0.099 ω -1 p . The initial magnetizations are the same as for the decaying turbulence scenario, σ 0 ≃ 1.6 and σ δB ≃ 13.
Turbulence is excited using a Langevin antenna scheme [47], in a way similar to the implementation of Refs. [13,15,19,20]. For the present 2D simulations, we excite external current fluctuations along the mean magnetic field (z axis) only, with wave modes oriented in the ðx; yÞ plane. We use 24 modes, with mean wave number hki ¼ ð2π=L x Þ × 2.9 (and similar hk 2 i 1=2 ), implying l c ∼ 350 c=ω p as before. Those external currents generate δB x and δB y . We then tune the amplitude of the antenna to FIG. 3. Histogram of correlation coefficients between the expected and observed evolution of γ 0 along chunks of test particle trajectories with initial 25 < γ < 50, for the 2D decaying turbulence PIC simulation. The chunks are selected at random among the whole of test particle histories, provided they fulfill the following criteria: the duration Δt ≃ 1l c =c and the energy change within that time interval verifies jΔγ 0 =γ 0 j > 1. This histogram shows that the parallel shear Θ k contribution, and more generally the nonresonant model as described by Eq. ( 5), match relatively well the observed variations. FIG. 4. Same as Fig. 3 (2D decaying turbulence PIC simulation), now considering the whole trajectory of each test particle with jΔγ 0 =γ 0 j > 2.
reproduce the chosen initial turbulent magnetization. The Langevin antenna is also characterized by a real frequency ω 0 and a damping term Γ 0 . We found it useful to set the real frequency to low values, in practice ω 0 ≈ 0, in order to avoid excessive heating of the plasma at early times, caused by the rapid generation of non-MHD electric fields on large scales. Regarding the damping term, we tune it in order to ensure that the autocorrelation time of the turbulent magnetic field matches roughly l c =c; in practice, we set Γ 0 ≃ 0.6hkic.
The power spectrum of magnetic fluctuations shown in Fig. 5 reveals a shape similar to that seen in the decaying turbulence case, with a (roughly) k -5=3 generic scaling over the inertial domain, followed by the steeper dissipative range at kinetic scales. The spectrum amplitude is more pronounced at the stirring scale, a trend which is characteristic of forced turbulence PIC simulations, if one compares Refs. [17,16]. This is expected insofar turbulence is continuously injected at the stirring scale in forced turbulence, while the spectrum moves in time to larger k in decaying turbulence.
As in the decaying turbulence scenario, we indicate with a dashed line the inverse gyroradius of particles with Lorentz factor γ ∼ 50. Again, their effective rigidity is of the order of 0.1, which falls in the right range to test the nonresonant acceleration model.
In Fig. 6, we plot the particle energy distribution, which reveals a power-law tail extending from γ ∼ 10 up to γ ∼ 10 3 , as for the decaying turbulence scenario. The best-fitting spectral index, s ≃ -2.2, is also close to that found previously. We note that in forced turbulence simulations, the spectrum evolves slowly in time, as the energy that is continuously injected in the simulation maintains δB=B (and, to a lesser degree, the overall magnetization) at values not far from its initial state, thereby guaranteeing that acceleration can proceed at all times. In decaying turbulence simulations, the drop in magnetization associated with magnetic dissipation within ∼5-10l c =c implies that acceleration becomes much slower, so that the spectrum essentially freezes on those timescales [14,17,18].
To reconstruct the probability density distributions of the correlation coefficients between observed and reconstructed histories, we follow test particles from t ¼ 1500 ω -1 p ∼ 4l c =c up to 5000 ω -1 p ∼ 14l c =c, as for decaying turbulence. We present those probability densities in Figs. 7 and8. The parameters (duration, amount of p ∼ 5l c =c. A power-law tail is clearly seen, extending from γ ∼ 10 up to γ ∼ 10 3 , with spectral index s ≃ -2, defining s through dN=dγ ∝ γ s . The test particles that we study, with 25 < γ < 50, are located in the power-law tail. FIG. 7. Histogram of correlation coefficients between the expected and observed evolution of γ 0 along chunks of test particle trajectories with initial 25 < γ < 50, for the 2D forced turbulence PIC simulation. This histogram shows that the nonresonant model provides a satisfactory match to the observed variations, and that the parallel shear term Θ k provides the dominant contribution to the force terms. variation of the energy) are the same as in the decaying turbulence scenario. We observe a similar trend, namely the nonresonant model captures fairly well the observed energy histories, and Θ k provides the dominant contribution among the three force terms. The perpendicular contribution Θ ⊥ and the inertial term do not show such significant degrees of correlation, although that of Θ ⊥ is skewed towards positive values, at least for short Δt ≃ 1l c =c timescales.
Generally speaking, the degree of agreement between model and simulations appears more satisfactory for the present forced turbulence scenario than for the decaying one. In that respect, we note that the shape of the spectrum can impact this comparison in the following way. The power spectrum of the forced simulation shows a larger amplitude on the smallest k modes, meaning on the largest length scales, rather than the decaying turbulence one, all things being considered equal. This difference can be read off Figs. 1 and5, but it is actually more pronounced at later times, since the power spectrum of the decaying turbulence shifts to larger k as time progresses. This implies that, on the whole, particles experience a turbulence on larger scales in the forced turbulence case than in the decaying turbulence one, as measured relatively to their gyroradius. Since the model works to order r g =l c , this larger degree of agreement is therefore not unexpected, at least at a qualitative level.
We now turn to 3D simulations. We have first performed a 3D forced turbulence PIC simulation with domain size N x × N y × N z ¼ 1080 3 cells, corresponding to physical size L x × L y × L z ¼ 540 3 c 3 =ω 3 p , integrated over 5000 time steps, corresponding to a time t ≃ 1500 ω -1 p ; the mesh size is now δx ¼ δy ¼ δz ¼ 0.5 c=ω p , δt ¼ 0.495 ω -1 p and we use 15 particles per species per cell. This choice of parameters is motivated by the need to optimize the execution time, while avoiding excessive shot noise associated with the number of macroparticles per skin depth volume. When measured in terms of the total relativistic plasma frequency, Ω p ¼ ½4πðn þ þ n -Þe 2 =w 1=2 , the mesh size reads δx ≃ 0.25 c=Ω p ; given that the plasma further heats with time in the turbulence, this provides a relatively fair sampling of the skin depth volume. Figure 9 offers a general view on the simulation at time t ≃ 600 ω -1 p : magnetic energy density (top panel), current density component along the mean field direction (middle panel) and plasma bulk velocity (bottom panel).
The initial mean field magnetization is σ 0 ¼ 1.6 as before, while σ δB ≃ 8. The forced turbulence is excited using the same Langevin antenna scheme as in 2D, with the following parameters: in 3D, we generate 24 modes of external current density fluctuations along x, along y, and along z separately, with mean wave numbers hki ¼ 2π=L max × 2.; for reference, hk 2 i 1=2 ≃ 2π=L max × 2 as well, guaranteeing a stirring scale l c ≃ L max =2 ≃ 270 c=ω p . We use a real frequency ω 0 ¼ 0 to avoid the generation of external electric fields in the forcing scheme, and a damping term Γ 0 ¼ 0.4hkic.
As before, we present in Fig. 10 the power spectrum of magnetic fluctuations and in Fig. 11 the energy distribution of particles, at t ≃ 600 ω -1 p ∼ 2.2l c =c. To compute the 3D power spectrum (and preserve memory usage), the field values have been rebinned by ten, so that the minimum length scale plotted is 10δ x ¼ 5 c=ω p . Consequently, the power spectrum shown in Fig. 10 lacks data at large wave numbers (in the dissipative range); it covers about two decades, even though the grid size contains 1080 cells along each its axis.
That timescale t ∼ 2.2l c =c is shorter than that used in 2D PIC simulations for plotting purposes, because of the shorter duration of that 3D simulation. Consequently, the peak amplitude associated with the externally injected energy appears more prominent in the 3D simulation, and the power-law tail of the energy distribution has not yet reached the maximum energy fixed by the coherence length, of the order of several hundreds here.
To compute the probability density histograms of the correlation coefficients r i , shown in Figs. 12 and 13, we have followed the test particle trajectories from t ≃ 500 ω -1 p ∼ 2l c =c up to 1500 ω -1 p ∼ 6l c =c, which marks the duration of the simulation. As anticipated, the PDF is sharply peaked around þ1 for this 3D simulation, indicating a nice match between the energy variations predicted by the model and those observed in the simulation. The parallel compression term Θ k provides as before the leading contribution; the PDF of the perpendicular force term is slightly biased toward positive values, as for the 2D forced turbulence simulation, while the inertial term does not show any clear signature, as in 2D. Interestingly, the correlation appears slightly enhanced when all force terms are taken FIG. 8. Same as Fig. 7 (2D forced turbulence PIC simulation), now considering the whole trajectory of each test particle. together as in Eq. ( 5), than when they are taken one by one, at least for the case in which intervals of duration 1l c =c are examined.
In Fig. 14, we present the energy evolution of two test particles, which are fair representatives of their parent population. The dotted blue line shows the evolution of γðtÞ, i.e., the Lorentz factor of the particle as measured in the simulation frame. It displays characteristic oscillations associated with the gyromotion of the particle around magnetic field lines that move at velocity v E : depending on the phase of that gyromotion, the particle motion is aligned or antialigned with v E , leading to a larger or smaller apparent energy in the simulation frame, see also Refs. [19,35]. The period of those oscillations thus provides an estimate of 2πr g =c, which takes different values at different times, depending on the strength of the magnetic field and of the particle energy. The solid purple line shows the evolution of γ 0 obs ðtÞ, in the frame R = E in which the motional electric field vanishes. The oscillations have disappeared and γ 0 obs ðtÞ evolves as the particle gains or loses energy through Fermi processes. Finally, the dashed red line shows the reconstructed particle history γ 0 th ðtÞ, using Eq. ( 5) with initial condition γ 0 th ðt 0 Þ ¼ γ 0 obs ðt 0 Þ at the initial time t 0 ¼ 1.8l c =c.
In the upper panel, we observe that the match between the reconstructed and the observed trajectories is rather tight in regions where the frequency of the oscillations increase, e.g., 3l c =c ≲ t ≲ 4l c =c. This is not unexpected, insofar as an increase in the frequency of oscillations corresponds to a decrease in the particle gyroradius, and the model works to order r g =l c . On the contrary, at later times t ≳ 5l c =c, the particle has achieved a larger energy, and it seemingly propagates in a region of lower-than-average magnetic strength, hence the ratio r g =l c is no longer small compared to unity, as evidenced by the timescale of the oscillations. Deviations from the observed trajectory can thus be expected at that stage, although they remain rather mild.
In the lower panel, the energy history is well reconstructed at early times t ≲ 3l c =c. We observe an offset in the vertical direction between the predicted and observed trajectories at later times, although those two histories maintain a rather strong degree of correlation. Had we chosen as initial time t 0 ≃ 3-3.5l c =c, we would thus have obtained a nice match to the observed history at late times. In effect, the departure between the model and the simulation is limited to the interval ∼2.7-3.3l c =c, and likely related to some small scale effect. As mentioned before, this observation has motivated our choice to adopt two timescales for the comparison of the model to the simulations: one reduced timescale of the order of 1l c =c, and one integrating over the whole history.
Finally, we compare the model to trajectories of test particles that were tracked in the 3D forced MHD simulation of the JHU turbulence database [45,46]. This 3D direct numerical simulation solves the incompressible MHD equations on a 1024 3 periodic grid with a time resolution δt ≃ 0.04δx. The database output provides 1024 time snapshots, with sampling interval 10δt. The simulation is visco resistive, with magnetic Prandtl number unity, and magnetic Reynolds number R λ ∼ 140 at the Taylor FIG. 12. Histogram of correlation coefficients between the expected and observed evolution of γ 0 along chunks of test particle trajectories with initial 25 < γ < 50, for the 3D forced turbulence PIC simulation. This histogram shows that the combination of all force terms provides a good match to the observed variations. Among the three force terms, the parallel shear Θ k provides the dominant contribution. FIG. 13. Same as Fig. 12 (3D forced turbulence PIC simulation), now integrating over the whole trajectory for each test particle. FIG. 14. Example of the time evolution of the energy for two test particles, in the 3D forced turbulence simulation. In dotted blue: the energy of the particle as measured in the simulation frame; in solid purple: the particle Lorentz factor γ 0 obs in the R = E frame; in solid red: the Lorentz factor γ 0 th , as reconstructed using Eq. ( 5). scale λ ∼ 1.05 × 10 -2 L; here, L represents the size of one side of the simulation cube and the Taylor scale is defined as
, where S k denotes the onedimensional power spectrum of magnetic fluctuations. The Alfvén velocity is v A ¼ 0.41c and the rms velocity hδu 2 i 1=2 ≃ 0.4c. The turbulence is excited through an external force acting on the velocity field at a stirring wave number k f ≃ 12.6L -1 . At the reference time t ¼ 0, the simulation, as made available on the database, has already achieved a steady state.
The integral scale of the turbulence, as defined in the database, is L w ∼ 0.1 in units of the cube size. The simulation volume thus comprises many coherence cells of the turbulence, hence the effective dynamic range is restricted to L w =δx ∼ 100. The power spectrum of magnetic fluctuations is shown in Fig. 15. It correspondingly reveals a lack of power at wave numbers k ≲ 3l -1 c followed by the k -5=3 scaling in the inertial range. We adopt here l c ¼ 0.1L.
We follow test particles with a gyroradius r g ≃ 2δx, in order to maintain r g =l c as small as possible while preserving a reasonable reconstruction of the particle trajectory. As can be seen from Fig. 15, the inverse gyroscale r -1 g lies at the transition between the inertial and the dissipative range. The effective rigidity is 2πr g =l c ∼ 0.1. Experiments conducted with a gyroradius twice as large provide similar results. We propagate 24000 particles over 4.2l c =c ∼ 200r g ; those particles were initialized with a common Lorentz factor (in the simulation frame), corresponding to the desired gyroradius, at random positions and velocity orientations. For each test particle, we integrate its trajectory over the duration of the simulation, using a numerical Monte Carlo code which, at each time step, queries the database to retrieve the values of the magnetic and velocity field at the particle location. The field values are determined at the particle spatial location using high-order (fourth or sixth) Lagrangian interpolation. Although, we sample the particle trajectory with a time step of 0.1r g =c, we do not seek to interpolate the field values at the corresponding intermediate times, and rather use the values calculated from the nearest snapshot. Given that the typical velocity on the grid size is of the order of ðδx=l c Þ 1=3 hδu 2 i 1=2 ∼ 0.050c-assuming a standard Kolmogorov scaling-this represents a reasonable approximation. This also allows us to maintain the computational time within reasonable limits, since computational time is here dominated by the queries to the database, which are performed online.
The code computes the electric field at the particle location using ideal Ohm's law, then advances the particle using a Boris pusher. All along the trajectories, we record the time and space derivatives of the magnetic and electric fields; the latter is computed from the magnetic and velocity derivatives. Those derivatives, provided by the database as fourth-order centered finite differencing, are used to calculate the quantities that enter the force terms in Eq. ( 5), as for the PIC simulation. The database does not directly provide time derivatives; those are thus calculated using first-order finite differencing from values obtained at consecutive times. We note that the force terms that enter Eq. ( 5) are dominated by the spatial derivatives in the sub-or mildly relativistic conditions of the present MHD simulation.
In Fig. 16, we plot the resulting energy distribution after a time t ≃ 4l c =c. It reveals a power-law tail at large momenta, as in the PIC simulation. To our knowledge, such a behavior had not been observed in time-evolving MHD simulations before. The spectral index s ≃ -4 is somewhat larger (in absolute value) than that observed in FIG. 15. Magnetic power spectrum of the 3D MHD simulation, rescaled by k 5=3 , vs wave numbers in units l -1 c . The inverse gyroradius of the test particles is indicated by a dashed red line. It lies at the transition between the inertial and the dissipative range. FIG. 16. Energy distribution of test particles propagated in the 3D MHD simulation, at a time t ≃ 4l c =c. At the initial time t ¼ 0, all particles shared a common Lorentz factor γ 0 , corresponding to a gyroradius r g ≃ 0.02l c . The spectrum takes a power-law shape at large energies, with spectral index s ≃ -4, assuming dN=dγ ∝ γ s . the PIC simulation, as expected for particle acceleration in a turbulence of smaller magnetization level [14,28].
In Figs. 17 and 18, we present the histograms of the correlation coefficients. As a threshold of energy variation and duration of integration, we have adopted jΔγ 0 j=γ 0 > 0.5 and Δt ≃ 1.7l c =c in a first case (Fig. 17), jΔγ 0 j=γ 0 > 1 and Δt ≃ 4.2l c in a second one (Fig. 18). The lesser threshold in energy variation and longer duration of the interval, comparatively to the PIC simulations, are meant to compensate for slower acceleration in the present simulation.
We recover here a high degree of correlation, as observed in the 3D PIC simulation. Interestingly, that degree of correlation is, in the present case, substantially higher when all force terms are combined together using Eq. ( 5) to reconstruct γ 0 th ðtÞ, than when they are taken individually. We also note that both Θ k and Θ ⊥ seem to provide contributions with a net positive degree of correlation, when taken individually, while the influence of a E • b is not visible here. Interestingly, the degrees of correlation of Θ ⊥ and Θ k appear to be of the same or order of magnitude, while the PIC simulations showed a clear dominance of Θ k . We cannot identify here the exact reason why that is so, but we speculate that this difference may indicate that their relative contribution depends on how the turbulence is driven: incompressible turbulence driven by external velocity fluctuations in the MHD case vs compressible turbulence driven by magnetic perturbations in the kinetic regime. Recent PIC simulations have similarly demonstrated that the efficiency of acceleration depends on the stirring procedure [23]. This difference may also be affected by the different velocity regimes (subrelativistic for MHD, relativistic for PIC).
With respect to the test of our model, we stress here the significance of observing such a significant degree of correlation for both 3D PIC and MHD simulations, up to the above difference in individual contributions: on the "large" length scales that we are interested in (comparatively to the kinetic scales), both should in principle reproduce the same physics of acceleration; however, both rely on different schemes of approximations. In particular, the MHD case neglects all kinetic effects and all deviations of Ohm's law that are inherently included in PIC simulations.
In the present work, we report on a statistical comparison of a recent model of nonresonant particle acceleration in magnetized turbulence vs 2D and 3D kinetic numerical simulations as well as a 3D (incompressible) MHD simulation. This model describes energization as the continuous interaction of the particle with the random velocity flow of the turbulence, in the frame of ideal MHD; it can thus be regarded as the direct generalization to a continuous turbulent flow of the original Fermi picture of discrete, pointlike interactions [28,34]. It does so by following the evolution of the particle momentum in the frame R = E that moves with the magnetic field lines at velocity v E ¼ E × B=B 2 , and where the electric field vanishes. This allows us to relate the sources of energy gains and losses to the gradients of the velocity field v E , and more particularly to three main contributions: an inertial term a E • b, a longitudinal shear term Θ k and a perpendicular compressive mode Θ ⊥ , the notions of longitudinal/perpendicular being defined relative to the mean magnetic field direction at that location. To lowest order in the ratio of particle gyroradius to coherence scale of the turbulence, r g =l c , the evolution of the particle energy is captured by Eq. (5).
To test this theoretical model, we have conducted PIC simulations of 2D decaying turbulence, of 2D and 3D driven turbulence in the relativistic regime v A ∼ c, and we have made use of the 3D forced MHD simulation of the FIG. 17. Histogram of correlation coefficients between the expected and observed evolution of γ 0 along blocks of the energy histories of test particles that have been propagated in the 3D MHD simulation. The initial Lorentz factor for all particles is γðt ¼ 0Þ ¼ 10 (simulation frame). FIG. 18. Same as Fig. 17 (3D MHD simulation), now considering the whole trajectory of each test particle.
JHU database. We have then followed the time histories of the energy for a large sample of particles and compared the observed time histories to those reconstructed by the model. For what regards the MHD simulation, we have propagated test particles through the simulation, properly taking into account the time evolution of the fields. In all simulations, we have selected particles whose inverse gyroradius corresponds to wave numbers at or below the transition between the inertial and the dissipative range of the turbulence, in order to test the model in conditions in which it applies, namely a ratio r g =l c as small as possible and near-MHD conditions. To obtain the reconstructed particle histories, we have extracted from the simulations the quantities a E • b, Θ k and Θ ⊥ then used Eq. ( 5) at each point along the particle trajectory to integrate in time the particle energy using Eq. ( 5). We have then computed for each particle trajectory a Pearson correlation test between the two histories (observed vs reconstructed), then derived from the sample of particles a probability density of the correlation coefficients r. A perfect adequation of the model to the date would translate in a probability density sharply peaked around þ1, while an complete inadequacy would rather yield a featureless, roughly uniform histogram over the interval ½-1; þ1. We have verified this using Monte Carlo simulations of test-particle transport in a synthetic turbulence composed of a sum of linear eigenmodes of the plasma, see the Appendix.
Our main result is that we observe a clearcut correlation between the model predictions and the numerical experiments, with histograms of the Pearson correlation coefficients distinctly peaked around þ1, for all numerical simulations. This indicates that the nonresonant model can successfully account for the bulk of particle energization through stochastic Fermi processes. Let us recall here that particle acceleration in a magnetized turbulence appears to proceed in two distinct stages: an injection into the nonthermal population through nonideal electric fields, then acceleration à la Fermi up to much higher energies [14,17]. The model and the test that we study here thus apply to the second stage, where the influence of nonideal electric fields can be neglected.
In our PIC numerical simulations, we observe that the longitudinal shear term Θ k appears to provide the dominant contribution to particle energization, because the correlation histogram when neglecting the other two force terms in the theoretical reconstruction of the energy histories lies close to that obtained when considering all force terms. This longitudinal shear term can be depicted as a form of slingshot acceleration in a moving, curved magnetic field, as in the Fermi type-B interaction of the original Fermi model [24]. Contrariwise, the MHD simulation reveals about similar degrees of correlation of Θ ⊥ and Θ k , with a slight preference for the former, which characterizes magnetic mirroring effects, or type-A Fermi interactions. This MHD simulation also shows a significantly higher degree of correlation when all contributions are summed together as in the model Eq. ( 5) than when only one force term is considered individually, and the other two discarded.
This difference in contributions between the MHD and the PIC simulations suggests that the physics of acceleration, in particular the dominant energization process, depends on the stirring process, on the nature of the turbulence and/or the velocity regime: while the (subrelativistic) turbulence of the MHD simulation is by construction incompressible and forced through solenoidal velocity motions, the (relativistic) turbulence in the PIC simulations is compressible and driven through external magnetic perturbations. A dependence of the energy distribution of accelerated particles on the stirring process (solenoidal vs compressible) has been noted before in Ref. [23].
More work is clearly needed to shed light on those issues and to develop further the above nonresonant model. It would be useful, in particular, to extend the microphysical random walk introduced in [28] in extract predictions for the particle energy spectrum, and to derive a kinetic equation for the distribution function. More work is also needed to understand how the present picture extends into the subrelativistic regime, to other domains of the turbulence cascade (e.g., dissipative) and to electron-ion plasmas. This will be the subject of future work.
In this appendix, we conduct comparisons similar to those discussed in the main text for PIC and MHD simulations, although using a synthetic turbulence constructed as a sum of noninteracting linear MHD eigenmodes. Quasilinear theory predicts that particle energization takes place through two types of interactions: gyroresonant interactions of the form
Gyroresonant interactions can take place for all modes, at least in the absence of local anisotropy à la Goldreich Sridhar [38,39], while the Landau synchrotron (transit-time damping) are specific to fast magnetosonic modes. We recall that those transit-time damping interactions are related to magnetic mirroring effects, which are captured by the theoretical model of nonresonant interactions (Θ ⊥ contribution). However, that model does not predict any gyroresonant interaction. Conversely, Alfvén modes predict, in the linear limit, Θ k ¼ 0 as well as Θ ⊥ ¼ 0, while fast magnetosonic modes lead to Θ k ¼ 0 but Θ ⊥ ≠ 0 [28].
The numerical code used to build the synthetic turbulence and track particles therein is presented in [35]. In brief, particle trajectories are integrated using a Bulirsch-Stoer algorithm. At each timestep, electromagnetic and velocity fields at the location of particles are constructed as the sum of a background magnetic field and the superposition of the fluctuations carried by a collection of waves with dispersion relation and polarizations of (special relativistic) MHD eigenmodes. The electric field is derived from the total magnetic field and total velocity field through ideal Ohm's law. The wave vectors and amplitudes of the waves are initialized so as to achieve the desired power spectrum of turbulence over a range of scales ½L min ; L max . Particles are injected along random directions in different turbulence realizations with the energy corresponding to the gyroradius of interest. To reconstruct the energy histories using Eq. ( 5), we calculate the spatial and temporal derivatives of the magnetic field and the velocity field, then derive those of the electric field through ideal Ohm's law. In this synthetic turbulence, the derivatives can be expressed analytically in terms of the plane wave expansion.
We conduct two experiments on such synthetic turbulence comprised of 256 modes, with wavelengths extending from L max ¼ l c down to L min ¼ L max =100. In experiment (A), we simulate a turbulence of isotropic fast magnetosonic modes with δB=B 0 ¼ 1, Alfvén velocity v A ¼ 0.6c, sound velocity v s ≪ v A , which implies a phase velocity for each wave v F ≃ v A . Isotropic means here that the turbulent magnetic power spectrum does not depend on the direction of the wave number; its scaling is assumed to follow Kolmogorov S k ∝ k -5=3 . We inject particles with a gyroradius outside the range of scales of the turbulence, r g ¼ 0.1L min . In that configuration, gyroresonant interactions are suppressed because restricted to high harmonics (large n) so that acceleration is dominated by transit time damping acceleration [48]. We thus expect the theoretical model to provide a fair reconstruction of the trajectories, at least its Θ ⊥ part. In a second experiment, (B), we simulate an opposite situation, namely a turbulence of Alfvén modes with v A ¼ 0.6c, δB=B 0 ¼ 1 and set the gyroradius of the particles to fall in the range of wavelengths of the turbulence, r g ¼ 0.1L max ¼ 10L min , which thus permits gyroresonant interactions at the first harmonic n ¼ 1. We simulate here simple Alfvén waves, meaning that we neglect any wave damping term and that we assume an isotropic Kolmogorov cascade. Our aim indeed is to bring to light the effect of gyroresonances, or rather, the lack of correlation between observed and reconstructed trajectories in a situation in which most of energy gain is known to result from gyroresonant interactions; we thus deliberately render those resonances sharp. We choose Alfvén waves in order to erase any magnetic mirroring effect. Consequently, we expect the theoretical model to behave poorly in that limit, given that it ignores such gyroresonances, by construction.
Figure 19 shows the power spectrum (normalized by k 5=3 ) of magnetic fluctuations in this synthetic wave turbulence, with the locations of r -1 g indicated as dashed lines for both models: model (A) with r g below the minimum scale, and model (B) with r g in the inertial range.
The histogram of the probability density of the Pearson correlation coefficients between the observed and reconstructed trajectories for model (A) is shown in Fig. 20. The concentration of the probability density of r around þ1 indicates that, as anticipated, the model is highly successful in reproducing the trajectories, most notably so for the magnetic mirror (Θ ⊥ ) contribution. A closer inspection reveals that the nonresonant model, summing over the contributions of all three force terms, provides a better match to the observed energy histories than the contribution of magnetic mirrors alone, as its probability density is more sharply peaked around þ1. The difference comes from the inertial term a E • b, which provides a net contribution, as evidenced by its overall positive degree of correlation. Within the frame of the model that we are testing, this is not altogether surprising, insofar as this inertial term characterizes the influence of accelerations/decelerations of the frame R = E in which the notion of a magnetic mirror can be properly defined. In that sense, the inertial term should not be left aside when considering the influence of Θ ⊥ (or Θ k , for similar reasons).
Figure 21 shows the corresponding histograms in the case of model (B), which appear relatively structureless and uniformly distributed over the interval ½-1; þ1. There appears to be a slight bias toward positive values of the correlation coefficients, but nothing of the sort discussed previously in Sec. III for PIC or MHD simulations. This suggests that most of the energy gains/losses are indeed not captured by the nonresonant model in the present case and that gyroresonant wave-particle interactions do not leave a strong signature in this histogram of correlation coefficients. The positive correlations observed in Sec. III can thus be interpreted as genuine evidence in favor of nonresonant acceleration. In that case, the model is unable to reproduce the energy gains, hence the PDF of correlation coefficients appears devoid of structure, revealing no particular preference for values close to þ1.
Reference[42] reports a spectrum characterized by q ≃
for fast modes in isothermal compressible, relativistic turbulence; however, this isothermal configuration, which is obtained by adding an ad hoc cooling term to the plasma evolution, does not match the conditions of the PIC simulations of relativistic turbulence which have reported D pp ∝ p 2 .
Available from: http://turbulence.pha.jhu.edu/Forced_MHD_ turbulence.aspx.
We select here a range in γ, not γ 0 , but this does not significantly influence our results, as we have explicitly checked for the 2D PIC driven turbulence and the MHD simulations discussed further below.