System-size Convergence of Nonthermal Particle Acceleration in Relativistic Plasma Turbulence

For many decades, turbulence has been recognized as a conceivable source of nonthermal energetic particles in collisionless plasmas. Theoretical works have proposed a variety of routes toward nonthermal particle acceleration (NTPA), including diffusive (second-order) acceleration from turbulent fluctuations (Alfvénic modes, Jokipii 1966; Schlickeiser 1989; Chandran 2000; Cho & Lazarian 2006; compressive modes, Schlickeiser & Miller 1998; Yan & Lazarian 2002; Chandran 2003; kinetic modes, e.g., Dermer et al. 1996; Fonseca et al. 2003; Petrosian & Liu 2004; Riquelme et al. 2017) and secular (first-order) acceleration via intermittent structures (shocks, Bykov & Toptygin 1982; Blandford & Eichler 1987; current sheets undergoing magnetic reconnection, Vlahos et al. 2004; Lazarian et al. 2012; Isliker et al. 2017; see also Beresnyak & Li 2016). These mechanisms of NTPA are tantalizing theoretical possibilities, but rely on various assumptions about the nature of turbulence and nonlinear plasma physics. Due the complexity and analytic intractability of the problem, the only practical way to prove the reality of turbulent NTPA (apart from direct experimental confirmation) is with self-consistent, large-scale numerical simulations.

Turbulent NTPA has important implications for space systems such as the solar corona, the solar wind, and planetary magnetospheres, as well as high-energy astrophysical systems such as pulsar wind nebulae, X-ray binaries, supernovae remnants, jets from active galactic nuclei, radio lobes, and gamma-ray bursts. Observations of broadband radiation spectra and cosmic rays imply that nonthermal particles are a significant component of the universe. In this work, we focus on relativistically hot plasmas with modestly relativistic bulk velocities.

In Zhdankin et al. (2017), we applied particle-in-cell (PIC) simulations to demonstrate that driven turbulence can produce a substantial population of nonthermal particles in relativistic pair plasmas, with power-law energy distributions that become harder with increasing magnetization σ (ratio of magnetic enthalpy to relativistic plasma enthalpy). However, these simulations also revealed that the distributions become softer with increasing system size, and were therefore unable to probe distributions in the asymptotic large-system limit. In principle, a lack of convergence can arise from inadequate scale separation, or from the adverse role of physical effects such as scale-dependent anisotropy, intermittency, damping of relevant (e.g., compressive) modes, or the inherent inefficiency of NTPA at magnetohydrodynamic (MHD) scales. Since supercomputers will be unable to simulate systems with sizes comparable to real astrophysical systems in the foreseeable future, it is necessary to understand the scaling of NTPA with system size before applying such simulations to model astrophysical phenomena.

In the present work, we address the system-size dependence of turbulent NTPA. Using an ensemble of PIC simulations with sizes extending beyond those in Zhdankin et al. (2017), we confirm a weak system-size dependence for nonthermal energy distributions when measured at a fixed number of large-scale dynamical times. More importantly, we present evidence that the distributions converge when compared at different times increasing with system size (approximately logarithmically or as a weak power law). Physically, this time dependence arises from the fact that the distributions do not fully develop until particles reach the highest accessible energies via Fermi acceleration. The converged index for the power-law energy distribution (α ≈ 3.0 for the given numerical setup and σ = 3/8) confirms turbulence as an efficient, viable astrophysical particle accelerator.

We perform simulations with the explicit electromagnetic PIC code Zeltron (Cerutti et al. 2013) using charge-conserving current deposition (Esirkepov 2001). The simulation setup is described in detail in Zhdankin et al. (2018). The domain is a periodic cubic box of size L³ (consisting of N³ cells) with uniform mean magnetic field ${{\boldsymbol{B}}}_{0}={B}_{0}\hat{{\boldsymbol{z}}}$. We initialize electrons and positrons from a uniform Maxwell–Jüttner distribution with combined particle density n₀ and temperature T₀ = θ₀mc², where m is the electron rest mass, and we choose θ₀ = 100 (giving an ultrarelativistic initial mean Lorentz factor γ₀ ≈ 300). We then drive strong (δB_rms ∼ B₀) turbulence at low wavenumbers (k = 2π/L) by applying a randomly fluctuating external current density (TenBarge et al. 2014). For these simulations, we fix the initial magnetization to ${\sigma }_{0}\equiv {B}_{0}^{2}/4\pi {h}_{0}=3/8$, where h₀ = 4n₀θ₀mc² is the initial relativistic enthalpy density. The initial Alfvén velocity is given by ${v}_{{\rm{A}}0}\equiv c{[{\sigma }_{0}/({\sigma }_{0}+1)]}^{1/2}\approx 0.52c;$ we perform each simulation for a duration of at least 7.5L/v_A0. As optimized by convergence studies, we set the initial Larmor radius to ρ_e0 = 1.5Δx (where Δx is the cell size), where the characteristic Larmor radius ${\rho }_{e}=\langle \gamma \rangle {{mc}}^{2}/{{eB}}_{\mathrm{rms}}$ is based on the instantaneous mean particle Lorentz factor $\langle \gamma \rangle$ and total rms magnetic field B_rms. We choose 64 particles per cell for the main simulations.

We perform a scan over system size L/2πρ_e0 by varying the number of cells in each simulation, taking N ∈ {256, 384, 512, 768, 1024, 1536}, so L/2πρ_e0 ∈ {27.2, 40.7, 54.3, 81.5, 109, 163}. To check reproducibility, we rerun all cases having N ≤ 1024 with a different random seed (for particle initialization and driving phases); for robustness, we analyze the particle distributions averaged for each simulation pair at a given size. We also perform statistical ensembles of sixteen 384³ cases and eight 768³ cases with 32 particles per cell to investigate statistical variation.

The simulations evolve as discussed in our previous papers (Zhdankin et al. 2017, 2018): external driving disrupts the thermal equilibrium and establishes turbulent fluctuations across a broad range of scales. Turbulence is fully developed by a few Alfvén times, after which internal energy increases at a constant rate due to turbulent dissipation. For reference, in Figure 1, we show the magnetic energy spectrum compensated by ${k}_{\perp }^{5/3}$, where k_⊥ is the wavenumber perpendicular to  B₀, for simulations of varying size, averaged over five snapshots from 3.1 ≤ tv_A0/L ≤ 5.2. Whereas the 512³ case has a spectrum steeper than ${k}_{\perp }^{-5/3}$, larger cases (768³ and above) have spectra close to ${k}_{\perp }^{-5/3}$, consistent with classical MHD turbulence theories (Goldreich & Sridhar 1995; Thompson & Blaes 1998). The 1536³ case exhibits an inertial range from k_⊥ρ_e ∼ 0.06 to k_⊥ρ_e ∼ 0.4.

Zoom In Zoom Out Reset image size

We now turn to the particle energy distribution f(γ), where the Lorentz factor γ = E/mc² is used interchangeably with energy E. The evolution of f(γ) for the 1536³ simulation is shown in the top panel of Figure 2. The distribution develops a power-law tail, f(γ) ∼ γ^−α, over several dynamical times, attaining an index α ≈ 3.0 at tv_A0/L ∼ 7. This power law extends from energies comparable to the instantaneous mean ($\langle \gamma \rangle \sim 1.5\times {10}^{3}$) up to energies limited by the system size (${\gamma }_{\max }\equiv {{LeB}}_{0}/2{{mc}}^{2}\sim 1.5\times {10}^{5}$), a factor of ∼50 in energy. At later times (not shown), particles accumulate at energies near γ_max, causing a high-energy pileup in the distribution.

Zoom In Zoom Out Reset image size

In the middle panel of Figure 2, we show the energy distributions at a fixed number of large-scale dynamical times, taken to be tv_A0/L = 7.0, for simulations of varying size (512³–1536³). For clarity, we compensate the distributions by γ³. The nonthermal tail steepens with increasing system size, ranging from an estimated index of α ≈ 2.7 for the 512³ case to α ≈ 3.0 for the 1536³ case. Thus, when compared at fixed times, there is no clear evidence for convergence of f(γ) with system size; although the scaling of α with size is weak (δα ∼ 0.3 for a factor of 3 increase in size), it can undermine the viability of turbulent NTPA in astrophysical systems if it persists to larger sizes.

The interpretation of the data changes, however, when the energy distributions are compared at different times, chosen to scale with system size. In the bottom panel of Figure 2, we show distributions at tv_A0/L ∈ {5.2, 5.8, 6.4, 7.0} for the simulations with {512³, 768³, 1024³, 1536³} cells (which is an approximately logarithmic increase of time with size). When compared at these times, the 768³, 1024³, and 1536³ simulations all exhibit converged distributions with an index near α ≈ 3.0, to within ±0.1 accuracy. Notably, these times approximately coincide with the initial formation of the pileup at γ_max. This leads to our main proposal, that turbulent NTPA produces a power-law particle energy distribution that converges with increasing system size, but the time required to fully form this distribution slowly increases with system size. We suggest that the apparent size dependence of the (fixed-time) index α in our simulations is due to the power laws becoming contaminated by the pileup at γ_max, which develops earlier for smaller systems. We next provide a physical motivation for this proposal.

To understand the acceleration process in detail, we tracked a random sample of 8 × 10⁵ particles in each simulation. In Figure 3, we show the energy evolution for the four tracked particles that attain the highest energies at the end of the 1536³ simulation. These four particles are the only tracked particles with final energies γ > γ_max. At early times, particle energies exhibit rapid oscillations on a timescale comparable to their Larmor period. The particles occasionally undergo acceleration episodes in which their energy rapidly increases by a factor of 2 or more, but the overall acceleration occurs gradually over several large-scale dynamical times.

Zoom In Zoom Out Reset image size

In the same figure, we show the average energy evolution for all tracked particles with final energies γ > 10⁵ (yielding 44 particles), denoted $\langle \gamma {\rangle }_{\mathrm{he}}$. After turbulence is fully developed, $\langle \gamma {\rangle }_{\mathrm{he}}$ increases at a slightly subexponential rate, until it approaches γ_max. As we now show, this evolution is consistent with Fermi acceleration with a slowly evolving acceleration timescale. We suppose that the scattering process causes the typical energy of the most energetic particles to increase as

$$\begin{eqnarray}&&\displaystyle \frac{d\gamma }{{dt}}\sim \displaystyle \frac{\gamma }{{\tau }_{\mathrm{acc}}},\end{eqnarray} \tag{ 1 }$$

where τ_acc(t) is the acceleration timescale. For second-order Fermi acceleration, assuming scattering by large-scale Alfvénic fluctuations, the acceleration timescale is

$$\begin{eqnarray}&&{\tau }_{\mathrm{acc}}\sim \displaystyle \frac{3}{4}\displaystyle \frac{{\lambda }_{\mathrm{mfp}}c}{{u}_{{\rm{A}}}^{2}},\end{eqnarray} \tag{ 2 }$$

where ${u}_{{\rm{A}}}={v}_{{\rm{A}}}/{(1-{v}_{{\rm{A}}}^{2}/{c}^{2})}^{1/2}={\sigma }^{1/2}c$ is the Alfvén four-velocity and λ_mfp is the scattering mean free path (Longair 2011). In particular, for scattering by classical MHD turbulence (i.e., a Goldreich–Sridhar spectrum of Alfvén waves), λ_mfp is expected to be effectively independent of energy and set by the largest scale of turbulence (e.g., Schlickeiser 1989; Miller et al. 1990; Chandran 2000). For time-independent τ_acc, Equation (1) leads to an exponential increase in particle energy. The time required for particles to reach γ_max from an initial energy of γ_i ∼ γ₀ is then $t/{\tau }_{\mathrm{acc}}\sim \mathrm{log}({\gamma }_{\max }/{\gamma }_{i})\,\sim \mathrm{log}(L/{\rho }_{e0})$. In relativistic plasmas with no energy sink, however, the acceleration timescale evolves with time since the Alfvén velocity decreases due to turbulent energy dissipation increasing the relativistic plasma inertia. As discussed in Zhdankin et al. (2018), for a constant energy injection rate, $\langle \gamma \rangle \sim {\gamma }_{0}(1+\eta {\sigma }_{0}{v}_{{\rm{A}}0}t/L)$ (where η ≈ 1 is the measured injection efficiency for the given simulations), v_A(t) is given by

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{v}_{{\rm{A}}}}{c} & = & \sqrt{\displaystyle \frac{\sigma }{\sigma +1}}={\left(1+\displaystyle \frac{\langle \gamma \rangle }{{\sigma }_{0}{\gamma }_{0}}\right)}^{-1/2}\\ & \sim & {\left(\displaystyle \frac{{c}^{2}}{{v}_{{\rm{A}}0}^{2}}+\eta \displaystyle \frac{{{tv}}_{{\rm{A}}0}}{L}\right)}^{-1/2}.\end{array}\end{eqnarray} \tag{ 3 }$$

The energy growth from second-order Fermi acceleration is then a power law in time (solving Equations (1)–(3)),

$$\begin{eqnarray}&&\gamma \sim {\gamma }_{i}{\left(1+\eta {\sigma }_{0}\displaystyle \frac{{{tv}}_{{\rm{A}}0}}{L}\right)}^{4{Lc}/3\eta {\lambda }_{\mathrm{mfp}}{v}_{{\rm{A}}0}}.\end{eqnarray} \tag{ 4 }$$

Note that, using ${(1+x/n)}^{n}\to \exp x$ as $n\to \infty$, this equation approaches an exponential in the limit ησ₀ ≪ 1, consistent with time-independent Fermi acceleration. We find that Equation (4) provides a good fit to $\langle \gamma {\rangle }_{\mathrm{he}}$, as shown in Figure 3, if we take λ_mfp/L = 1/2 and γ_i = 0.7γ₀ (giving γ ∝ t^5.1 at late times). These results suggest that second-order Fermi acceleration can account for NTPA observed in the simulations; measuring the relative importance of first-order and second-order mechanisms is left to future work.

Inverting Equation (4) gives the time required for the particle to reach a given energy γ,

$$\begin{eqnarray}&&\displaystyle \frac{{{tv}}_{{\rm{A}}0}}{L}\sim \displaystyle \frac{1}{\eta {\sigma }_{0}}\left[{\left(\displaystyle \frac{\gamma }{{\gamma }_{i}}\right)}^{3\eta {\lambda }_{\mathrm{mfp}}{v}_{{\rm{A}}0}/4{Lc}}-1\right].\end{eqnarray} \tag{ 5 }$$

This is used to estimate the time taken for Fermi-accelerated particles to reach the system size limit (γ ∼ γ_max). Note that Equation (5) approaches a logarithmic function for ησ₀ ≪ 1.

We now relate Fermi acceleration to the late-time evolution of the energy distributions. As previously discussed, the distributions form a power law and subsequently develop a broad pileup near γ_max. It is natural to focus on the distribution just prior to the pileup formation, when the power law has its maximum extent. After the power law is fully formed, a pair of inflection points appear due to the pileup (which makes the distribution no longer concave down). We define the inflection time, t_inf, as the latest time at which the difference between local power-law indices $[\alpha (\gamma )\equiv -\partial \mathrm{log}f/\partial \mathrm{log}\gamma ]$ at the two inflection points [local extrema of α(γ)] is less than 0.1. For t > t_inf, the distributions become influenced by the pileup, making a power-law index difficult to define precisely.

We show the normalized inflection time t_infv_A0/L versus system size L/2πρ_e0 in the top panel of Figure 4. We find that t_infv_A0/L increases with size, consistent with particles requiring a longer number of larger-scale dynamical times to reach γ_max. In fact, t_inf is consistent with the second-order Fermi acceleration timescale calculated in Equation (5) (with L_mfp/L = 1/2 and γ_i = γ₀, giving t_inf ∝ (γ_max/γ₀)^0.2), also shown in the top panel of Figure 4. This scaling is close to logarithmic over the given sizes. We note that the time taken for the distribution to reach energies slightly beyond γ_max exhibits a similar scaling as for t_inf (not shown). The system-size dependence of the inflection time, and related pileup, motivates comparing f(γ) at times that increase according to Equation (5).

Zoom In Zoom Out Reset image size

In the bottom panel of Figure 4, we show the index α (measured at the logarithmic center of the power-law segment) versus system size L/2πρ_e0 at various times: the inflection time t = t_inf, logarithmic times $t\propto \mathrm{log}(L/2\pi {\rho }_{e0})$, arbitrary fixed time t = 7L/v_A0, and times with fixed mean particle energy $\langle \gamma \rangle \sim 4.2{\gamma }_{0}$ (nominally identical to fixed time, but sensitive to statistical variations). We find that measuring the distribution at fixed time or fixed injected energy shows a clear system-size dependence, although the dependence weakens with size; in particular, α exhibits an approximately logarithmic dependence on L/2πρ_e0 (weaker than suggested in Zhdankin et al. 2017). In contrast, the distribution at the inflection time or logarithmic times shows no systematic variation for L/2πρ_e0 ≳ 80. To sum up, distributions attain the same power-law index, independent of system size, just prior to pileup formation.

We conclude with a comment about the statistical significance of our results. In our simulations, the amount of energy injected from external driving fluctuates randomly in time, since driven mode phases are evolved randomly. While the mean energy injection rate approaches a universal value for sufficiently long simulations, the distributions presented in this paper were measured after a limited duration (≲7L/v_A0), and thus the amount of injected energy at that point can vary significantly between different runs (by up to ∼30%). In principle, a larger energy injection may supply a harder nonthermal population, bringing up an important question: do the measured nonthermal distributions exhibit significant statistical variability (from run to run) due to random driving? To build confidence that our largest simulations are not statistical outliers, we analyzed ensembles of sixteen 384³ (L/2πρ_e0 = 40.7) and eight 768³ (L/2πρ_e0 = 81.5) simulations. For the 384³ ensemble, we obtain an average fixed-time index of $\langle \alpha \rangle \approx 2.53$ and rms spread of δα_rms ≈ 0.06 (at t = 7L/v_A0), with α weakly correlated with injected energy. For the 768³ ensemble, we find $\langle \alpha \rangle \approx 2.79$ and δα_rms ≈ 0.07, indicating that the statistical spread is similar at both sizes and is less than the size dependence. When measured at logarithmic times, we find $\langle \alpha \rangle \approx 2.86$, δα_rms ≈ 0.09 for the 384³ ensemble (at tv_A0/L ≈ 4.2) and $\langle \alpha \rangle \approx 2.94$, δα_rms ≈ 0.06 for the 768³ ensemble (at tv_A0/L ≈ 5.2), in agreement with Figure 4. In addition, we find that the variations are modest in each pair of equal-size simulations from our main system-size scan. Hence, we believe that the trends in our simulations are robust.

Based on the simulations presented in this paper, we are prepared to declare the system-size independence of power-law indices of particle energy distributions produced by driven, large-scale turbulence in relativistic collisionless plasmas. We empirically observe that convergence occurs at times that depend on system size (approximately logarithmically or as a weak power law). We propose a physical interpretation for this dependence: as the size is increased, an increasing number of scatterings is required for particles to acquire the highest accessible energies via Fermi acceleration. The energy distribution becomes fully developed only once particles reach the energy limit γ_max due to finite system size, and subsequently the distributions experience a high-energy pileup that may complicate measurements of the power law.

The pileup forms in our simulations due to continuous energy injection into a closed system, which is a convenient benchmark scenario but unrealistic in nature. In general, particle energies may be limited by other physics, such as open boundaries or radiative cooling, which may cause a similar pileup (Schlickeiser 1984; Stawarz & Petrosian 2008). At present, we can only postulate that, regardless of the physical mechanism that sets the high-energy cutoff, distributions will attain the converged index measured in this Letter during transient evolution.

This work demonstrates that accurately measuring the converged power-law index is feasible with modestly large PIC simulations of turbulence (with L/2πρ_e ≳ 80). We find a converged index α ≈ 3.0 at magnetization σ = 3/8 and turbulence amplitude δB_rms ∼ B₀; as discussed in Zhdankin et al. (2017), NTPA becomes more efficient with increasing σ, leading to smaller α and faster formation of the nonthermal population (as implied by Equation (5)). Whether the converged index is universal with respect to other system parameters remains to be seen. A similar converged value of the index, for comparable plasma parameters, is measured in PIC simulations of relativistic magnetic reconnection (e.g., Werner & Uzdensky 2017); steeper indices are measured for relaxation of magnetostatic equlibria (Nalewajko et al. 2016). Likewise, our measured index is consistent with to the value ≈3.2 inferred from the continuum synchrotron spectrum in the Crab Nebula, where a magnetization below unity is expected (Meyer et al. 2010).

We thank referee Jonathan Zrake for valuable comments. The authors acknowledge support from NSF grant AST-1411879 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. This work also used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. This work used the XSEDE supercomputer Stampede2 at the Texas Advanced Computer Center (TACC) through allocation TG-PHY160032 (Towns et al. 2014).

Software: Zeltron (Cerutti et al. 2013).