First-principles demonstration of diffusive-advective particle acceleration in kinetic simulations of relativistic plasma turbulence

TL;DR

This study uses kinetic simulations to validate the Fokker-Planck framework for nonthermal particle acceleration in relativistic plasma turbulence, revealing diffusive and advective behaviors consistent with astrophysical observations.

Contribution

It provides the first ab-initio kinetic simulation evidence supporting the Fokker-Planck model of particle acceleration in relativistic plasma turbulence.

Findings

Diffusive nature of NTPA confirmed through particle tracking.

Energy diffusion coefficient scales as γ^2 at high energies.

Advection coefficient influences particle energy distribution.

Abstract

Nonthermal relativistic plasmas are ubiquitous in astrophysical systems like pulsar wind nebulae and active galactic nuclei, as inferred from their emission spectra. The underlying nonthermal particle acceleration (NTPA) processes have traditionally been modeled with a Fokker-Planck (FP) diffusion-advection equation in momentum space. In this paper, we directly test the FP framework in ab-initio kinetic simulations of driven magnetized turbulence in relativistic pair plasma. By statistically analyzing the motion of tracked particles, we demonstrate the diffusive nature of NTPA and measure the FP energy diffusion ($D$) and advection ($A$) coefficients as functions of particle energy $\gamma m_e c^2$. We find that $D(\gamma)$ scales as $\gamma^2$ in the high-energy nonthermal tail, in line with 2nd-order Fermi acceleration theory, but has a much weaker scaling at lower energies. We also find that $A$ is not negligible and reduces NTPA by tending to pull particles towards the peak of the particle energy distribution. This study provides strong support for the FP picture of turbulent NTPA, thereby enhancing our understanding of space and astrophysical plasmas.

Equations

$$\gamma (t)={\gamma }_D{\gamma }^{\prime }\left(1+{\beta }_D\frac{v_{\perp }^{\prime }}{c}\cos {\omega }^{\prime }{t}^{\prime }\right),$$

$${\partial }_{t}f={\partial }_{\gamma }(D{\partial }_{\gamma }f)-{\partial }_{\gamma }(Af).$$

$$\Delta \overline{\gamma }({\gamma }_0,\Delta t)$$

$$\delta {\gamma }_{\text{rms}}({\gamma }_0,\Delta t)$$

Full text

acceleration of particles, turbulence

\setabbreviationstyle

[acronym]long-short

1 Introduction

Relativistic plasmas with nonthermal power-law energy distributions are ubiquitous in astrophysical systems such as \acPWN Meyer et al. (2010); Bühler & Blandford (2014), jets from \acAGN Begelman et al. (1984); Hartman et al. (1992) and their radio lobes Hardcastle et al. (2009), and black-hole accretion-disk coronae Yuan et al. (2003). The underlying \acNTPA processes have been studied theoretically for decades; proposed mechanisms include collisionless shocks Blandford & Eichler (1987), turbulence Kulsrud & Ferrari (1971), and magnetic reconnection Hoshino & Lyubarsky (2012). The most common turbulent \acNTPA models posit that particles gain energy in a stochastic process (e.g., scattering off magnetic fluctuations) that can be modeled using \pglsFP advection-diffusion equation in momentum space Fermi (1949); Kulsrud & Ferrari (1971); Melrose (1974); Skilling (1975); Blandford & Eichler (1987); Schlickeiser (1989); Miller et al. (1990); Chandran (2000); Cho & Lazarian (2006).

Numerical tests of the \acFP framework for \acNTPA were originally performed by injecting test particles into \aclMHD simulations Dmitruk et al. (2003, 2004); Kowal et al. (2012); Lynn et al. (2014); Kimura et al. (2016); Isliker et al. (2017a) or artificially prescribed fields Arzner et al. (2006); O’Sullivan et al. (2009). These test-particle simulations are relatively inexpensive, but have physical limitations such as ad-hoc particle injection and the absence of particle feedback on the fields, which can only be resolved by considering more physically complete simulations.

Recently, first-principles kinetic (and hybrid) \acPIC simulations have confirmed that turbulence Kunz et al. (2016); Makwana et al. (2017); Zhdankin et al. (2017, 2018b, 2019); Comisso & Sironi (2018); Arzamasskiy et al. (2019), shocks Hoshino et al. (1992); Amato & Arons (2006); Spitkovsky (2008); Sironi & Spitkovsky (2011); Marcowith et al. (2016), and relativistic reconnection Zenitani & Hoshino (2001); Jaroschek et al. (2004); Lyubarsky & Liverts (2008); Sironi & Spitkovsky (2014); Guo et al. (2014, 2016); Werner et al. (2016, 2018); Werner & Uzdensky (2017) can generate efficient \acNTPA in collisionless plasma. \acPIC simulations contain complete microphysical information including the self-consistent trajectories and energy histories of individual particles. However, this wealth of data has not yet been employed directly to test stochastic acceleration models (e.g., \acFP) or to measure the energy diffusion and advection coefficients.

In this paper, we use tracked particles to demonstrate stochastic acceleration and directly measure the \acFP coefficients in \ac3D \acPIC simulations of driven turbulence in collisionless relativistic plasma. Stochastic particle acceleration in relativistic plasma turbulence has important applications to astrophysical systems such as \acPWN Bucciantini et al. (2011); Tanaka & Asano (2017), \acAGN accretion flows Dermer et al. (1996); Kimura et al. (2015), \acAGN jets Rieger et al. (2007); Asano et al. (2014), and \aclGRB Dermer & Humi (2001). We consider pair plasma both for theoretical and computational simplicity, and for its relevance to high-energy astrophysical systems like \acPWN and \acAGN jets. However, our methods also apply to future investigations of \acNTPA in turbulent non-relativistic and electron-ion plasmas, as well as to other processes, e.g., magnetic reconnection.

2 Method

We analyze 3D simulations [performed with our \acPIC code Zeltron Cerutti et al. (2013)] of externally driven turbulence in relativistic pair plasma Zhdankin et al. (2018b). We focus on the largest simulation with $1563^{3}$ grid cells and 64 particles per cell (electrons and positrons combined), totalling ${\sim}$2.4\text{\times}{10}^{11}$$ particles; smaller simulations give similar results. The simulation domain is a periodic cube of size $L$, with an initially uniform magnetic guide field $B_{0}\mathbf{\hat{z}}$. The plasma is initially uniform and isotropic, with total charged particle density $n_{0}$, and a Maxwell-Jüttner thermal distribution with a relativistically hot temperature of $T_{0}=\nobreak 100m_{e}c^{2}$, corresponding to the average Lorentz factor $\overline{\gamma}_{\textrm{init}}\approx\nobreak 3T_{0}/m_{e}c^{2}=\nobreak 300$. The initial magnetization is $\sigma_{0}=B_{0}^{2}/16\pi n_{0}T_{0}=\nobreak 3/8$. In the fiducial simulation, the normalized system size is $L/2\pi\rho_{e_{0}}=\nobreak 163$, where $\rho_{e_{0}}\equiv\nobreak\overline{\gamma}_{\textrm{init}}m_{e}c^{2}/eB_{0}$ is the initial characteristic Larmor radius. Turbulence is electromagnetically driven TenBarge et al. (2014) and becomes fully developed after several light crossing times Zhdankin et al. (2018a), with rms turbulent magnetic fluctuations $\delta B_{\textrm{rms}}\sim\nobreak B_{0}$. The turbulence is essentially Alfvénic Zhdankin et al. (2018a), with initial Alfvén velocity $v_{A0}/c\equiv[\sigma_{0}/(\sigma_{0}+1)]^{1/2}\simeq 0.52$.

Our previous studies Zhdankin et al. (2017, 2018b) have shown that such turbulence reliably produces nonthermal power-law particle spectra. In this section, we describe our procedure to investigate \acNTPA in relation to the \acFP framework. First, we examine gyro-scale oscillations in particle energy and explain their physical origin. Then, we present our methodology for averaging out these oscillations, which is critical for accurately measuring energy diffusion. Finally, we detail our tests of diffusive \acNTPA in our turbulence simulations, and our procedure for measuring the energy diffusion and advection coefficients as functions of particle energy.

Our analysis tracks the positions, momenta, and local electromagnetic field vectors for a statistical ensemble of $8\text{\times}{10}^{5}$ randomly-chosen particles. We observe order-unity oscillations in particle energy, $\gamma{}m_{e}c^{2}$, at the gyrofrequency, as shown for a representative particle in Fig. 1a. The energy oscillates once per gyro-orbit (Fig. 1b) because of the large-scale electric field accelerating and decelerating the particle as it gyrates. To describe this analytically, we consider a charged particle moving in constant, uniform electromagnetic fields. We use primed variables for the frame moving with the $\mathbf{E\times B}$ drift velocity, $\mathbf{v}_{D}$, given by $c\mathbf{v}_{D}/(c^{2}+v_{D}^{2})=\mathbf{E\times B}/(E^{2}+\nobreak B^{2})$.

In the primed frame, where $\mathbf{B}^{\prime}$ and $\mathbf{E}^{\prime}$ are parallel, the particle gyrates about $\mathbf{B}^{\prime}$ while being accelerated along $\mathbf{B}^{\prime}$ by $\mathbf{E}^{\prime}$. Typically, $E^{\prime}\ll\nobreak B^{\prime}$, and so $\gamma{{}^{\prime}}m_{e}c^{2}$ is slowly-varying on the oscillation timescale. Then, the motion in the primed frame is approximately a simple gyration with $\mathbf{E}^{\prime}\approx 0$, and, applying the inverse Lorentz transformation, the lab-frame energy can be found:

[TABLE]

with $t=t_{0}+\gamma_{D}(t^{\prime}+\beta_{D}v_{\perp}^{\prime}\sin{\omega^{\prime}t^{\prime}}/\omega^{\prime}c)$. Here, $t$ is the coordinate time, $v_{\perp}^{\prime}$ is the particle’s primed-frame velocity perpendicular to $\mathbf{B}^{\prime}$, $\beta_{D}=\nobreak v_{D}/c$, $\gamma_{D}=(1-\nobreak v_{D}^{2}/c^{2})^{-1/2}$, $\omega^{\prime}=\nobreak eB^{\prime}/\gamma^{\prime}m_{e}c$ is the cyclotron frequency with $B^{\prime}=\nobreak B/\gamma_{D}$, and $t_{0}$ is a phase. Since we are considering relativistic particles ($v_{\perp}^{\prime}\sim\nobreak c$) and relativistic turbulence ($E_{\textrm{rms}}\sim\nobreak B_{0}$ and $\beta_{D}\sim\nobreak 1$), the predicted oscillation magnitude is comparable to $\gamma$ (see Fig. 1a), and cannot be ignored.

Without further processing, these energy oscillations are incompatible with a random-walk-based energy-diffusion model. However, we aim to measure the statistical properties of \acNTPA on the Alfvénic timescale ${\sim}L/v_{A}$ relevant to the formation of the nonthermal power law. Since this is generally much longer than the gyroperiod, we analyze just the secular component of the energy histories.

Our oscillation removal procedure is informed by Eq. (1), which indicates that in the idealized case of uniform constant fields, the secular component of the lab-frame energy is $\gamma_{D}\gamma^{\prime}$. However, $\mathbf{v}_{D}$ fluctuates as the particle traverses small-scale fields, so we average $\mathbf{v}_{D}$ over the smoothed gyroperiod $2\pi{\gamma^{\prime}}m_{e}c/e{B}$, where $\gamma^{\prime}=\gamma_{D}\gamma(1-\mathbf{v}_{D}\cdot\mathbf{v}/c^{2})$ is obtained from boosting the lab-frame four-velocity $\gamma\mathbf{v}$ by $\mathbf{v}_{D}$. We denote by $\langle\mathbf{v}_{D}\rangle$ the average of $\mathbf{v}_{D}$ over this smoothed period. We then define the smoothed particle energy to be $\langle\gamma_{D}\rangle\langle\gamma^{\prime}\rangle$, where $\langle\gamma^{\prime}\rangle=\langle\gamma_{D}\rangle\gamma(1-\langle\mathbf{v}_{D}\rangle\cdot\nobreak\mathbf{v}/c^{2})$ and $\langle\gamma_{D}\rangle\equiv(1-\langle\mathbf{v}_{D}\rangle^{2}/c^{2})^{-1/2}$. This transformed energy (Fig. 1a) has greatly reduced oscillations. Thus, this procedure extracts the secular component of particle energy, allowing us to test the \acFP picture of \acNTPA. Hereafter, $\gamma$ and “energy” refer to $\langle\gamma_{D}\rangle\langle\gamma^{\prime}\rangle$, except in the case of the overall particle energy distributions (Fig. 4a) and magnetic energy spectra (Fig. 4b).

We now describe our tests of diffusive acceleration using tracked particles. We use $t_{0}=\nobreak 10.0L/c$ as a fiducial initial time for our measurements, by which point a power-law particle energy spectrum has fully formed. We then bin tracked particles by their energy at $t_{0}$ in logarithmic intervals spaced by 10%. For each bin, we measure the standard deviation, $\delta\gamma_{\textrm{rms}}$, and the mean, $\overline{\gamma}$, of the particle energy distribution as a function of subsequent times $\Delta t\equiv\nobreak t-t_{0}$. For a classical diffusion process, one expects $\delta\gamma_{\textrm{rms}}(\Delta t)\propto\nobreak\sqrt{\Delta t}$ if the diffusion coefficient $D(\gamma,t)$ varies slowly compared with $\Delta t$ and $\delta\gamma_{\textrm{rms}}$.

We then measure the diffusion and advection coefficients, $D(\gamma)$ and $A(\gamma)$, respectively, for the simplest Fokker-Planck (FP) equation for the energy distribution $f(\gamma,t)$, ignoring pitch angle (with respect to $\mathbf{B}$):

[TABLE]

Limiting our measurements to times where $\Delta t\lesssim\nobreak{}L/v_{A}$, $\delta\gamma_{\textrm{rms}}\ll{\gamma_{0}}$, and $\Delta\overline{\gamma}\equiv\nobreak\overline{\gamma}(t)-\nobreak\overline{\gamma}(t_{0})\ll{\gamma_{0}}$, we approximate the bin distribution as narrow and the coefficients as constant in time. Applying Eq. (2), we find:

[TABLE]

We first measure $D(\gamma)$ and $M(\gamma)$ by applying Eqs. (3) and (4) to each energy bin and then calculate $A(\gamma)\equiv\nobreak M-\partial_{\gamma}D$.

3 Results

We first describe the evolution of the overall lab-frame distribution, $f(\gamma)$. Starting from a thermal distribution, $f(\gamma)$ acquires a power-law tail extending to the system-size limit, $\gamma_{\rm max}\equiv\nobreak LeB_{0}/2m_{e}c^{2}\simeq\nobreak$1.5\text{\times}{10}^{5}$$, and gradually hardening over time—its index converges to approximately -3 by $12.3L/c$. At the start of our measurements, $t_{0}=10.0L/c$, the index is approximately $-3.2$, the peak of $f(\gamma)$ is at $\gamma_{\textrm{peak}}\simeq 520$, and the mean at $\gamma_{\textrm{avg}}\simeq 1170$. Because the system lacks an energy sink, $\gamma_{\textrm{peak}}$ and $\gamma_{\textrm{avg}}$ increase at a rate of about $40c/L$ and $100c/L$, respectively.

We now present tests of energy diffusion. For illustration, Fig. 2a shows the evolution of the energy distribution of a single bin of particles with ${\gamma_{0}}=\nobreak$5\text{\times}{10}^{3}$$, deep in the power-law section. We find that $\overline{\gamma}\propto\Delta t$ for small $\Delta t$ (Fig. 2b), while $\delta\gamma_{\textrm{rms}}(\Delta t)\propto\nobreak\sqrt{\Delta t}$ (Fig. 2c), consistent with simple diffusion.

Fig. 3 shows $\delta\gamma_{\textrm{rms}}(\Delta t)$ for several bins along with corresponding $\sqrt{\Delta t}$ fits. To avoid artifacts of the smoothing procedure, each fit begins after one gyroperiod $T_{\rm L}(\gamma_{0},B_{\textrm{rms}})\equiv\nobreak 2\pi{\gamma_{0}}m_{e}c/e{B_{\textrm{rms}}}$. To ensure Eqs. (3) and (4) are valid, each fit ends when $\delta\gamma_{\textrm{rms}}/{\gamma_{0}}$ reaches 0.3, $\Delta\overline{\gamma}/{\gamma_{0}}$ reaches 0.1, or $\Delta t=\nobreak 2L/c$, whichever is earliest. Under these criteria, almost all fits end before $\Delta t=1L/c$. The fits generally agree well with the data over the fitted intervals. While some of the plotted $\delta\gamma_{\textrm{rms}}(\Delta t)$ have intervals of weakly anomalous energy diffusion these regimes are not our current focus, and we refer to studies of anomalous energy diffusion in plasma turbulence Isliker et al. (2017a, b). In summary, Figs. 2 and 3 confirm our expectations of a standard diffusive process, supporting the \acFP model of turbulent \acNTPA.

We now report on our measurements of $D(\gamma)$ and $A(\gamma)$. Fig. 4c shows $D(\gamma)$, extracted from the fits of $\delta\gamma_{\textrm{rms}}(\Delta t)$ using Eq. (4). In the high-energy nonthermal tail (2\text{\times}{10}^{3}$\lesssim\nobreak\gamma\lesssim\nobreak$3\text{\times}{10}^{4}, see Fig. 4a), $D=\nobreak 0.06(c/L)\gamma^{2}$ is an excellent fit, while for lower energies $\gamma\lesssim\gamma_{\textrm{peak}}$, there is a much shallower scaling roughly consistent with $D\propto\gamma^{2/3}$. We observe that $D\propto\gamma^{2}$ for particles gyroresonant with fluctuations in the inertial range of the magnetic energy spectrum (Fig. 4b), while the lower-energy scaling corresponds to the sub-inertial range of turbulence. In simulations with different magnetization, system size, and number of particles per cell (not shown), the high-energy scaling of $D\propto\gamma^{2}$ is maintained while the low-energy behavior varies slightly.

The high-energy scaling of $D\propto\gamma^{2}$ is commonly predicted by \acNTPA theories Kulsrud & Ferrari (1971); Skilling (1975); Blandford & Eichler (1987); Schlickeiser (1989); Chandran (2000); Cho & Lazarian (2006). We compare our high-energy fit, $D/\gamma^{2}=\nobreak 0.06c/L$, to the theoretical prediction from 2nd-order Fermi acceleration, $D/\gamma^{2}=\nobreak u_{A}^{2}/3c\lambda_{\textrm{mfp}}$, for ultra-relativistic particles interacting with isotropic scatterers moving at the Alfvén velocity, where $u_{A}=v_{A}/(1-v_{A}^{2}/c^{2})^{1/2}$ and $\lambda_{\textrm{mfp}}$ is the mean free path between scattering events Blandford & Eichler (1987); Longair (2011). At $t_{0}=\nobreak 10.0L/c$, $v_{A}=\nobreak 0.51c$, giving a theoretical scaling of $D/\gamma^{2}=\nobreak 0.12c/\lambda_{\textrm{mfp}}$, which agrees with our fit if $\lambda_{\textrm{mfp}}\sim\nobreak 2L$. We also find, from simulations with reduced driving wavelength (not shown), that the turbulence driving scale, and not the box size, sets $D/\gamma^{2}$ (and hence $\lambda_{\textrm{mfp}}$). While the value of $D/\gamma^{2}$ is consistent with 2nd-order Fermi acceleration from gyroresonant scattering of particles by Alfvén modes, further work is needed to test this view. This includes studying the effect of varying magnetization, as well as the spatial transport characteristics.

To compare the effects of the first- and second-order derivative terms in Eq. (2), we separate the contributions of $A$ and $\partial_{\gamma}D$ to the average acceleration rate $M(\gamma)=\partial_{t}\overline{\gamma}$. We extract $M(\gamma)$ from linear fits of $\overline{\gamma}(\Delta t)$ [see Eq. (3)], using the same time intervals as those used for fitting $\delta\gamma_{\textrm{rms}}(\Delta t)$ to measure $D(\gamma)$. As shown in Fig. 4d, $M(\gamma)$ is positive, as expected with external energy injection, and has a minimum near $\gamma_{\textrm{avg}}\approx\nobreak 1200$.

We then compute the contribution of energy diffusion to acceleration, $\partial_{\gamma}D$, (Fig. 4d) and subsequently the advection coefficient $A(\gamma)\equiv M-\partial_{\gamma}D$ (Fig. 4e). In the high-energy power-law section ($\gamma>\gamma_{\textrm{peak}}$), $A$ is negative, while for low energies ($\gamma<\gamma_{\textrm{peak}}$), $A$ is positive. Overall, $A$ tends to pull particle energies towards $\gamma_{\textrm{peak}}$, narrowing $f(\gamma)$. We find that $A$ is reasonably approximated by a logarithmic scaling with energy. For a momentum-space \acFP equation containing only diffusion [see, e.g., Ramaty (1979)], the expected energy-space advection coefficient for Eq. (2) is $A=2D/\gamma$. The measured negative logarithmic scaling thus implies that the momentum-space \acFP equation contains a significant advection term. In the nonthermal section, the magnitude of $A$ is generally smaller than, but still comparable to that of $\partial_{\gamma}D$. Hence, the evolution of the nonthermal population in this simulation cannot be interpreted as being due to $D$ alone. This also complicates estimates of the acceleration time based on $D/\gamma^{2}$. The rate of overall energy increase from advection is $\int d\gamma Af\approx-23c/L$ while that from diffusion is $\int d\gamma(\partial_{\gamma}D)f\approx 140c/L$. The general effects of systematic and stochastic acceleration (due to $A$ and $D$ respectively) are discussed in, e.g., Pisokas et al. (2018) and Vlahos & Isliker (2019). We note that our measurement of $A$ has considerable scatter outside of the central range, $10^{2}\lesssim\nobreak\gamma\lesssim\nobreak 10^{4}$, as it depends on the difference between two noisy quantities.

Finally, we test whether the \acFP equation with coefficients measured by our methodology can reproduce the evolution of $f(\gamma,t)$ from the \acPIC simulations. We first repeat the measurements of $D$ and $A$ at almost 200 different $t_{0}$ evenly spaced over the entire simulation, thus obtaining time-dependent coefficients. We then insert the measured coefficients into the \acFP equation [Eq. (2)] and solve it numerically using a finite-volume method. We use linear interpolation in $\gamma$ and nearest-neighbor interpolation in $t$. The \acFP coefficients are extrapolated as constant in $\gamma$ for energies without enough particles to measure them. Comparing with the particle energy distributions obtained from the \acPIC simulation, we find that the \acFP equation gives excellent agreement at all subsequent times when initialized with the corresponding \acPIC distribution as early as $t=\nobreak{}4.5L/c$ (Fig. 5a), with moderate errors if initialized at $t=\nobreak{}0$ (Fig. 5b). This may be due to the turbulence not being fully developed at early times. Finally, evolving from $t=\nobreak{}4.5L/c$ with $A$ set to zero shows a clear mismatch with the \acPIC distributions (Fig. 5c), giving us confidence that the measurements of $A$ are reasonably accurate and that $A$ significantly affects the evolution of $f(\gamma,t)$. Thus, in this regime of fully-developed turbulence, we find that the \acFP equation with our coefficient measurement methodology is appropriate for modeling \acNTPA. We leave further investigation of the early-time behavior to future work.

4 Conclusions

In this study we rigorously demonstrate, for the first time, diffusive \acNTPA in first-principles \acPIC simulations of driven relativistic plasma turbulence, through direct statistical measurements using large numbers of tracked particles. We introduce a procedure to suppress large-amplitude gyro-oscillations of particle energy, which is critical for revealing the diffusive nature of \acNTPA and measuring the \acFP coefficients. We find that the energy diffusion coefficient $D$ scales with particle energy $\gamma{}m_{e}c^{2}$ as $D\simeq\nobreak 0.06(c/L)\gamma^{2}$ in the high-energy nonthermal power-law region, in line with theoretical expectations Kulsrud & Ferrari (1971); Skilling (1975); Blandford & Eichler (1987); Schlickeiser (1989); Chandran (2000); Cho & Lazarian (2006), while there is a much shallower scaling at energies below the peak of the energy distribution. We also measure the energy advection coefficient $A(\gamma)$, though with more uncertainty. We find that $A$ is not negligible, and tends to narrow the distribution by accelerating low-energy particles and decelerating high-energy particles. Furthermore, a numerical solution of the \acFP equation with the measured coefficients reproduces the evolution of the particle energy spectrum from the \acPIC simulation over a significantly longer time interval than was used for measuring the coefficients. This suggests that this simple \acFP model can fully account for \acNTPA in our simulations. These results thus lend strong first-principles numerical support to a broad class of turbulent \acNTPA theories.

Our new methodology can also be applied to future tracked-particle studies of \acNTPA in other contexts such as shocks or magnetic reconnection, and over broader ranges of physical regimes. Future work may investigate the effects of system parameters such as magnetization, plasma beta, and guide field strength; radiative cooling; relativistic vs non-relativistic regimes; and plasma composition (e.g., pair vs electron-ion plasma). In addition, the analysis can be extended to include pitch-angle dependence and scattering. Characterizing stochastic \acNTPA in various regimes has important implications in a broad range of contexts such as solar flares, the solar wind, \aclPWN, \aclAGN, \aclGRB, and cosmic ray acceleration in supernova remnants. Thus, this study will facilitate more detailed tests of \acNTPA theories against \acPIC simulations exploring various physical situations, thereby advancing our understanding of space, solar, and high-energy astrophysical phenomena.

The authors acknowledge support from NSF grants AST-1411879 and AST-1806084, and NASA ATP grants NNX16AB28G and NNX17AK57G. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357.