The Interplay of Magnetically Dominated Turbulence and Magnetic Reconnection in Producing Nonthermal Particles

Turbulence structures from our reference 2D simulation are illustrated in Figure 1. We plot the magnetic field squared fluctuations δB², the out-of-plane electric current density J_z, the bulk dimensionless four-velocity Γβ, and the particle density ratio n/n₀. Here ${\rm{\Gamma }}=1/\sqrt{1-{({\boldsymbol{V}}/c)}^{2}}$ indicates the plasma bulk Lorentz factor and $\beta =| {\boldsymbol{V}}| /c$ is the dimensionless plasma bulk speed obtained by averaging the velocities of individual particles. We can see that the fluctuations δB² are generally stronger in large-scale flux tubes (see the circular structures of size comparable to the energy-carrying scale l), but high values of δB² are also obtained in small-scale structures identified with reconnection plasmoids (see the circular structures with size ≪l). These are "secondary" magnetic islands (flux ropes in 3D) that are produced by magnetic reconnection (Biskamp 2000). In such plasmoids, the particle number density n exhibits strong enhancements in excess of n ∼ 15 n₀. High values of particle number density occur in large-scale flux tubes as well. In general, the particle density displays strong compressions in coherent quasi-circular structures spanning a range of scales.

Zoom In Zoom Out Reset image size

In between flux tubes, reconnection layers reveal the formation of plasmoids within narrow current sheets. Indeed, current sheets with high aspect ratio tend to fragment into plasmoids and secondary current sheets as a result of magnetic reconnection. Smaller-size current sheets are also ubiquitous, spanning a wide range of scales. We will see in the following sections that reconnecting current sheets, which are a natural by-product of turbulent cascades in magnetized plasmas (e.g., Servidio et al. 2009; Wan et al. 2013; Cerri & Califano 2017; Franci et al. 2017; Haggerty et al. 2017; Comisso & Sironi 2018; Dong et al. 2018; Papini et al. 2019), play an important role for particle injection into the acceleration process (Comisso & Sironi 2018). Finally, we also point out that in the strongly magnetized regime of plasma turbulence investigated here, the plasma bulk speed can reach very high values. In particular, we observe ultrarelativistic flows with bulk Lorentz factor as high as Γ ∼ 5. Such high speeds develop predominantly in between the large-scale flux tubes, although high-velocity fluctuations occur all over the spatial domain.

We now consider the fluid structures that develop in 3D plasma turbulence. Our reference 3D simulation has L/d_e0 = 820, which is half the size of the reference 2D simulation. However, since in 3D we adopt perturbation numbers up to N_max = 4 (as compared to N_max = 8 in 2D), we still have a well-extended turbulence inertial range. In fact, the ratio of the initial energy-carrying scale l = 2π/k_N to the plasma skin depth d_e0 remains the same between our reference 2D and 3D simulations, leading to the same high-energy cutoff of the particle energy spectrum (see Comisso & Sironi 2018, as well as Equation (9) in the following subsection).

The turbulence structures from our reference 3D simulation are displayed in Figure 2. The magnetic field squared fluctuations δB² present both large-scale and small-scale structures. However, in this case, the large-scale fluctuations are not organized in coherent flux tubes (as it was in 2D, where they were a result of the constrained 2D dynamics). Despite differences in the large-scale structure of the magnetic field, there is still a copious presence of current sheets (current ribbons when considering the third direction). Due to the presence of the mean magnetic field ${{\boldsymbol{B}}}_{0}={B}_{0}\hat{{\boldsymbol{z}}}$, current ribbons are mostly elongated along $\hat{{\boldsymbol{z}}}$. We can see that the z component of the electric current density, J_z, displays a variety of current sheets of different sizes. Some of these current layers break into plasmoids (see also Section 4), as highly elongated layers cannot be stable against the plasmoid instability, also in 3D geometry (e.g., Daughton et al. 2011; Sironi & Spitkovsky 2014; Guo et al. 2015; Huang & Bhattacharjee 2016; Ebrahimi 2017; Werner & Uzdensky 2017; Baalrud et al. 2018; Stanier et al. 2019). Here we show that plasmoids/flux ropes are self-consistently created in fully 3D plasma turbulence (see Section 4), where current sheets are self-consistently generated by the turbulence itself. As for 2D plasma turbulence, we will see that these current sheets play an important role in the initial stages of particle acceleration (Sections 4–6).

Zoom In Zoom Out Reset image size

Locations characterized by strong electric current densities are typically accompanied by strong gradients in particle density. In localized regions, the particle density can exceed n ∼ 12 n₀, similar to the 2D case. On the other hand, large-scale structures like the overdense regions at the core of 2D large-scale flux tubes are missing in 3D. Finally, we observe that also in 3D, due to the high magnetization of the system, the plasma flow speed is generally very high. We can see regions with ultrarelativistic flow speeds having bulk Lorentz factor as high as Γ ∼ 4.

We now present the time evolution of the magnetic power spectrum from the reference 2D and 3D simulations. In our simulations, turbulence develops from the initialized magnetic fluctuations. The magnetic energy decays in time, as no continuous driving is imposed, and a well-developed inertial range and kinetic range of the turbulence cascade develop within the outer-scale nonlinear timescale. In Figure 3, we show the time evolution of the magnetic power spectrum P_B(k) for the reference 2D simulation, where

$$\begin{eqnarray}&&{P}_{B}(k){dk}=\displaystyle \sum _{{\boldsymbol{k}}\in {dk}}\displaystyle \frac{\delta {{\boldsymbol{B}}}_{{\boldsymbol{k}}}\cdot \delta {{\boldsymbol{B}}}_{{\boldsymbol{k}}}^{* }}{{B}_{0}^{2}}\end{eqnarray} \tag{ 4 }$$

is computed from the discrete Fourier transform $\delta {{\boldsymbol{B}}}_{{\boldsymbol{k}}}$ of the fluctuating magnetic field. Each curve refers to a different time (from brown to orange), as indicated by the corresponding vertical dashed lines in the inset, where we present the temporal decay of the magnetic field fluctuations δB_rms/B₀. We can see that at MHD scales (kd_e0 ≲ 0.5) the magnetic power spectrum is consistent with a Kolmogorov scaling P_B(k) ∝ k^−5/3 (Biskamp 2003) (compare with the dotted–dashed line), while the Iroshnikov–Kraichnan scaling P_B(k) ∝ k^−3/2 (Iroshnikov 1963; Kraichnan 1965) (triple-dotted–dashed line) is possibly approached at late times. At kinetic scales (kd_e0 ≳ 0.5), the spectrum steepens and approaches a power-law slope P_B(k) ∝ k^−4.3 (compare with the dashed line). A similar slope was proposed for magnetized turbulence at subinertial scales in a cold plasma (Abdelhamid et al. 2016; Passot et al. 2017). We finally observe that the turbulence integral scale

$$\begin{eqnarray}&&{\ell }(t)=\displaystyle \frac{2\pi }{{k}_{I}(t)}=2\pi \displaystyle \frac{{\int }_{0}^{\infty }{k}^{-1}{P}_{B}(k,t){dk}}{{\int }_{0}^{\infty }{P}_{B}(k,t){dk}},\end{eqnarray} \tag{ 5 }$$

which is close to the energy-carrying scale associated with the wavenumber where P_B(k, t) peaks, increases as the magnetic energy decays in time. This is due to the merging of the large-scale magnetic flux tubes, which drives an inverse energy transfer to scales larger than the initial integral scale (e.g., Biskamp & Schwarz 2001).

Zoom In Zoom Out Reset image size

We now consider the time evolution of the magnetic power spectrum in 3D. Due to the presence of the large-scale mean magnetic field ${{\boldsymbol{B}}}_{0}$, turbulence becomes increasingly anisotropic toward small scales, within the inertial range. To account for this global anisotropy with respect to ${{\boldsymbol{B}}}_{0}$, we consider the magnetic power spectrum with respect to the wavenumber ${k}_{\perp }={({k}_{x}^{2}+{k}_{y}^{2})}^{1/2}$ perpendicular to the mean field, obtained from the discrete Fourier transform of the fluctuating magnetic field as

$$\begin{eqnarray}&&{P}_{B}({k}_{\perp }){{dk}}_{\perp }=\displaystyle \sum _{{\boldsymbol{k}}\in {{dk}}_{\perp }}\displaystyle \frac{\delta {{\boldsymbol{B}}}_{{\boldsymbol{k}}}\cdot \delta {{\boldsymbol{B}}}_{{\boldsymbol{k}}}^{* }}{{B}_{0}^{2}}.\end{eqnarray} \tag{ 6 }$$

Figure 4 shows the time evolution of the magnetic power spectrum ${P}_{B}({k}_{\perp })$, which does exhibit inertial and kinetic ranges of the turbulence cascade at times ct/l ≳ 1. As the magnetic field fluctuations decay (see inset), the inertial range (${k}_{\perp }{d}_{e0}\lesssim 0.5$) of the magnetic power spectrum tends to flatten from ${P}_{B}({k}_{\perp })\propto {k}_{\perp }^{-5/3}$ (Goldreich & Sridhar 1995; Thompson & Blaes 1998; dotted–dashed line) to ${P}_{B}({k}_{\perp })\propto {k}_{\perp }^{-3/2}$ (Boldyrev 2006; triple-dotted–dashed line). At kinetic scales (${k}_{\perp }{d}_{e0}\,\gtrsim 0.5$), the spectrum steepens to a power law ${P}_{B}({k}_{\perp })\propto {k}_{\perp }^{-4.3}$ (dashed line), similar to the 2D result and in agreement with theoretical predictions for magnetized turbulence at subinertial scales in cold plasmas (Abdelhamid et al. 2016; Passot et al. 2017). Note also that in the 3D case the magnetic energy decays faster than in the 2D case (compare insets of Figures 3 and 4). We will show that this leads to a reduced particle acceleration rate at late times.

Zoom In Zoom Out Reset image size

The most interesting outcome of the turbulent cascade is the generation of a large population of nonthermal particles. This is shown in Figure 5 (for the 2D setup), where the time evolution of the particle energy spectrum ${dN}/d\mathrm{ln}(\gamma -1)$ is presented ($\gamma -1={E}_{k}/{{mc}}^{2}$ is the normalized particle kinetic energy). As a result of turbulent field dissipation, the spectrum shifts to energies much larger than the initial Maxwellian, which is shown by the blue line peaking at $\gamma -1\sim {\gamma }_{{th}0}-1\simeq 0.6$. At late times, when most of the turbulent energy has decayed, the spectrum stops evolving (orange and red lines): it peaks at γ − 1 ∼ 5 and extends well beyond the peak into a nonthermal tail of ultrarelativistic particles that can be described by a power law

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{d\gamma }={N}_{0}\,{\left(\displaystyle \frac{\gamma -1}{{\gamma }_{{st}}-1}\right)}^{-p},\quad \mathrm{for}\ {\gamma }_{{st}}\lt \gamma \lt {\gamma }_{c},\end{eqnarray} \tag{ 7 }$$

and a sharp cutoff for γ ≥ γ_c. Here N₀ is the normalization of the power law and p is the power-law index, which is about 2.8 for the simulation results presented in the main panel of Figure 5 (note that in our figures we plot dN/dln(γ − 1) to emphasize the particle content, proportional to ${(\gamma -1)}^{-p+1}d\gamma$ for the distribution in Equation (7)). The percentage of the particles in the nonthermal tail (measured as the number of particles with Lorentz factor exceeding twice the thermal peak) is high, ζ_nt ∼ 17%, and it corresponds to a high value of the normalization N₀, which is close to the thermal peak. The starting point of the power law, γ_st, is roughly only a factor of two larger than the peak of the particle energy spectrum at late times. Therefore, dropping O(1) factors, the starting point of the power law can be estimated as

$$\begin{eqnarray}&&{\gamma }_{{st}}\sim {\gamma }_{\sigma }=\left(1+\displaystyle \frac{{\sigma }_{0}}{2}\right){\gamma }_{{th}0},\end{eqnarray} \tag{ 8 }$$

since most of the magnetic energy is converted to particle energy by the time the particle energy spectrum has saturated (see inset of Figure 3). On the other hand, the high-energy cutoff γ_c depends on the system size. As discussed in the following sections, stochastic acceleration by turbulent fluctuations dominates the energy gain of the most energetic particles. High-energy particles cease to be efficiently scattered by turbulent fluctuations when their Larmor radius ${\rho }_{L}=\left(\gamma {mc}/{eB}\right){v}_{\perp }$ exceeds the integral length scale ${\ell }=2\pi /{k}_{I}$, implying an upper limit to their Lorentz factor of

$$\begin{eqnarray}&&{\gamma }_{c}\sim e\sqrt{\langle {B}^{2}\rangle }\displaystyle \frac{{\ell }}{{{mc}}^{2}}\sim \displaystyle \frac{2\pi }{{k}_{I}{d}_{e0}}\sqrt{{\sigma }_{z}}{\gamma }_{{th}0},\end{eqnarray} \tag{ 9 }$$

where $\langle {B}^{2}\rangle$ is the space-averaged mean-square value of the magnetic field and

$$\begin{eqnarray}&&{\sigma }_{z}=\displaystyle \frac{{B}_{0}^{2}}{4\pi {n}_{0}{w}_{0}{{mc}}^{2}}={\sigma }_{0}{\left(\displaystyle \frac{{B}_{0}}{\delta {B}_{\mathrm{rms}0}}\right)}^{2}.\end{eqnarray} \tag{ 10 }$$

This argument assumes that the turbulence survives long enough to allow the particles to reach this upper limit (we also assumed B₀/δB_rms ≳ 1). A numerical confirmation of Equation (9), with k_I ∼ k_N, was presented in Comisso & Sironi (2018) by performing simulations with different domain sizes. We point out that inverse magnetic energy transfer (e.g., Biskamp & Schwarz 2001; Zrake 2014; Brandenburg & Kahniashvili 2015) can possibly drive a substantial decrease in time of k_I, which, in turn, can allow the most energetic particles to reach even higher energies.

Zoom In Zoom Out Reset image size

We observe that the slope of the power law is not universal, but it depends on the magnetization σ₀ and the ratio δB_rms0/B₀ (Comisso & Sironi 2018). The inset of Figure 5 shows how the power-law index changes with the magnetization σ₀ from two series of simulations having δB_rms0/B₀ = 1 and δB_rms0/B₀ = 2. We can see that the slope of the power law becomes harder for larger magnetization, and that for fixed σ₀ it is harder when increasing δB_rms0/B₀ (see also Figure 7). The decrease of the power-law index p for increasing magnetization σ₀ (see also Zhdankin et al. 2017; Comisso & Sironi 2018) is in analogy with the results of PIC simulations of relativistic magnetic reconnection (Guo et al. 2014; Sironi & Spitkovsky 2014; Werner et al. 2016; Lyutikov et al. 2017; Petropoulou & Sironi 2018). We will see that magnetic reconnection plays an important role also in the turbulence scenario considered here. However, as we show below, its role is confined to the initial stages of particle acceleration, while the dominant acceleration process is given by stochastic scattering off turbulent fluctuations, which determines the slope and the cutoff of the high-energy power-law tail.

A similar picture holds in 3D, i.e., a generic by-product of the magnetized turbulence cascade is the production of a large number of nonthermal particles. Figure 6 shows the evolution of the particle energy spectrum ${dN}/d\mathrm{ln}(\gamma -1)$ starting from the initial Maxwellian peaked at $\gamma -1\sim {\gamma }_{{th}0}-1\simeq 0.6$. As time progresses, the particle energy spectrum shifts to higher energies and develops a high-energy tail containing a large fraction of particles. At late times, when most of the turbulent energy has decayed, the particle energy spectrum stops evolving (orange and red lines), and it peaks at γ − 1 ∼ 7. It extends well beyond the peak into a nonthermal tail of ultrarelativistic particles that can be described by a power law with an index p ∼ 2.9 (main panel of Figure 6). As in the 2D case, the normalization of the power law is close to the peak of the spectrum, giving a large fraction of nonthermal particles. At ct/l = 12 we find that about 16% of particles have or exceed twice the energy of the spectral peak, which provides an indication of the percentage of particles in the nonthermal tail ζ_nt.

Zoom In Zoom Out Reset image size

In order to understand the dependence of the high-energy power-law slope on the initial magnetization in 3D, we performed four large-scale 3D simulations with ${\sigma }_{0}\in \left\{5,10,20,40\right\}$ and same δB_rms0/B₀ = 1, L/d_e0 = 820. The power-law index p decreases for increasing σ₀ (see top inset in Figure 6), with values that are close to the ones from the corresponding 2D simulations with δB_rms0/B₀ = 1 (blue curve from the inset in Figure 5). Here we also show the scaling of the high-energy cutoff γ_c (bottom inset in Figure 6), defined as the Lorentz factor where the spectrum drops one order of magnitude below the power-law best fit. The high-energy cutoff γ_c increases as ${\gamma }_{c}\propto {\sigma }_{0}^{1/2}$ (compare with dashed line in the inset), which is consistent with the expectation from Equations (9) and (10) for a σ₀-independent domain size L/d_e0 and fixed δB_rms0/B₀.

Several astrophysical systems are thought to have δB_rms/B₀ larger than unity (e.g., $\delta {B}_{\mathrm{rms}}^{2}/{B}_{0}^{2}\sim 6$ in some regions of the Crab Nebula; Lyutikov et al. 2019). Therefore, we have performed three additional 2D simulations with initial ratios δB_rms0/B₀ = 1, 2, 4, with fixed initial magnetization σ₀ = 40 and a larger domain size L/d_e0 = 3280. Figure 7 shows that the power law becomes harder with increasing δB_rms0/B₀, with p < 2 for large initial fluctuations. In this case, both Equations (8) and (9) should be understood as upper limits that are subject to energy constraints, as we now discuss. The starting point of the power-law tail, γ_st, could be lower than indicated in Equation (8), if only a minor fraction of the available energy goes into thermal particles, while most of the energy goes into the nonthermal tail. Also, while in the case p > 2 one can have from Equation (9) that ${\gamma }_{c}\to \infty$ as k_Id_e0 → 0, the case 1 < p < 2 has a lower attainable high-energy cutoff γ_c, since the mean energy per particle in the power-law tail has to be (Sironi & Spitkovsky 2014)

$$\begin{eqnarray}&&\displaystyle \frac{1-p}{2-p}\,\displaystyle \frac{{({\gamma }_{c}-1)}^{2-p}-{({\gamma }_{{st}}-1)}^{2-p}}{{({\gamma }_{c}-1)}^{1-p}-{({\gamma }_{{st}}-1)}^{1-p}}=\left(1+\chi \displaystyle \frac{{\sigma }_{0}}{2}\right){\gamma }_{{th}0},\end{eqnarray} \tag{ 11 }$$

where χ is the fraction of turbulent magnetic energy converted into particles belonging to the power-law tail.

Zoom In Zoom Out Reset image size

We conclude this section with the results of 2D simulations having different initial plasma temperature θ₀. From Figure 8, we can see that the slope p, the fraction of nonthermal particles, and the extent of the nonthermal tail γ_c/γ_st do not depend on θ₀. Indeed, this plot shows that spectra obtained from simulations with different θ₀ nearly overlap, when shifted by an amount equal to the initial thermal Lorentz factor γ_th0. The role of the initial choice of temperature is only to produce an energy rescaling, since both γ_st and γ_c are proportional to γ_th0, as can be seen from the relations (8) and (9), and the definitions of σ₀ and $\sqrt{{\sigma }_{z}}/{d}_{e0}$ already take into account relativistic thermal effects.

Zoom In Zoom Out Reset image size

Up to this point, we have discussed general features of the particle spectrum generated as a by-product of the plasma turbulence. We have found that despite some differences between 2D and 3D settings, the produced particle spectrum does not depend on the dimensionality of the simulation domain (see also Comisso & Sironi 2018). In both cases, the high-energy power-law range extends from (about) the thermal peak to a maximum energy set by the energy-containing scale of turbulence. These common features, combined with the fact that the slope of the power law is also similar, yield a similar percentage of particles in the power-law tail. In the next sections, we will shed light on the particle acceleration mechanisms that produce the nonthermal particle spectrum.

We begin our analysis from the reference 2D case, and then we extend the analysis to the reference 3D case. For the injection analysis presented in this section, we employed a subsample of ∼10⁶ tracked particles for the 2D case and a subsample of ∼10⁷ tracked particles for the 3D case.

We show in Figure 9(a) the time evolution of the Lorentz factor for 10 representative particles that eventually populate the nonthermal tail at ct/l = 12 (see particle spectrum in Figure 5). These particles have a distinct moment in which they are "extracted" from the thermal pool at γ ∼ γ_th and injected to higher Lorentz factors γ ≫ γ_th. To identify this moment, which we call injection time t_inj, we evaluate when the rate of increase of the particle Lorentz factor (averaged over cΔt/d_e0 = 45) satisfies ${\rm{\Delta }}\gamma /{\rm{\Delta }}t\geqslant {\dot{\gamma }}_{\mathrm{thr}}$, and prior to this time the particle Lorentz factor was γ ≤ 4γ_th0 ∼ 6. We take the threshold ${\dot{\gamma }}_{\mathrm{thr}}\simeq 0.01\sqrt{{\sigma }_{0}}{\gamma }_{{th}0}{\omega }_{p0}$, but we have verified that our identification of t_inj is nearly the same when varying ${\dot{\gamma }}_{\mathrm{thr}}$ around this value by up to a factor of three (the factor 0.01 is much lower than the typical collisionless reconnection rate [∼0.1, in units of the Alfvén speed], which is the appropriate reference scaling here, as shown in Comisso & Sironi 2018 and below).

Zoom In Zoom Out Reset image size

Once t_inj is determined for the population of particles at hand, it is possible to explore the properties of the electromagnetic fields at the injection location. In this case, by analyzing the fields at injection, we find that the out-of-plane current density J_z is particularly revealing. In particular, J_z has, in general, high values at injection locations. To provide a statistical measure of the likelihood of this occurrence, we can construct the probability density function (pdf) of the magnitude of the out-of-plane electric current density experienced by the particles at their injection time, $| {J}_{z,p}|$, normalized by J_z,rms, i.e., the standard deviation of the current density ${J}_{z,\mathrm{rms}}={\left\langle {J}_{z}^{2}\right\rangle }^{1/2}$ in the whole domain at that time. The outcome of this analysis is shown in Figure 9(b) by the red circles. The pdf of the high-energy particles at injection should be contrasted with the pdf of the entire population of particles at a representative time (here, ct/l = 3.5), shown by the blue diamonds in Figure 9(b). The difference between the two pdf's is striking. The main difference is that the pdf of the overall particle population is peaked around zero, while the pdf of the high-energy particles at injection is peaked at much higher values corresponding to $| {J}_{z,p}| \sim 4\,{J}_{z,\mathrm{rms}}$. In particular, ∼95% of the high-energy particles are injected at locations with $| {J}_{z,p}| \geqslant 2\,{J}_{z,\mathrm{rms}}$. On the other hand, by taking all the particles at the representative time ct/l = 3.5, only ∼9% of them happen to be in regions where $| {J}_{z,p}| \geqslant 2\,{J}_{z,\mathrm{rms}}$. Note also that the pdf of the overall particle population does not follow Gaussian statistics owing to the intermittent nature of current sheets in turbulence (e.g., Servidio et al. 2009; Cerri et al. 2017; Haggerty et al. 2017; Dong et al. 2018), which is therefore reflected in the pdf of the particles that sample the entire domain.

To obtain further insight, we look now at the morphology of regions with out-of-plane current density $| {J}_{z}| \geqslant 2\,{J}_{z,\mathrm{rms}}$, and we correlate it with the spatial locations of the particles undergoing injection at t_inj. This is shown in Figure 10(a), where we can see that the vast majority of the structures with $| {J}_{z}| \geqslant 2\,{J}_{z,\mathrm{rms}}$ are sheet-like structures, namely, current sheets, and the overwhelming majority of particles at injection reside in these regions. A large fraction of these current sheets are active reconnection layers, fragmenting into plasmoids. A typical case of such reconnecting current sheets is illustrated in Figure 10(b), where we show a small portion of the domain, corresponding to the area within the rectangular blue contour in Figure 10(a), at different times ct/l = 3.3, 3.4, 3.5. The reconnecting current sheet evolves in time and breaks up in shorter sheets owing to the formation of plasmoids. During this period of time, particles are constantly injected up to nonthermal energies, as shown by the red circles in Figure 10(b).

Zoom In Zoom Out Reset image size

These results are also robust in 3D, for which we have performed the same type of analysis. Figure 11(a) shows the time evolution of the Lorentz factor for 10 representative particles selected to end up in different energy bins of the nonthermal tail at ct/l = 12 (see particle spectrum in Figure 6). As in 2D, we can see a sudden acceleration episode with particles that are extracted from the thermal pool and injected into the acceleration process. We identify the injection time t_inj as for the 2D case, by evaluating when the Lorentz factor increases at a rate exceeding the same threshold ${\dot{\gamma }}_{\mathrm{thr}}$ adopted for 2D, starting from a value that is γ ≤ 5γ_th0 ∼ 8 (this value is slightly higher than the 2D case, since in 3D a larger fraction of magnetic energy is dissipated by the end of the simulation). Note that, as in 2D, after the injection phase the particles continue to gain energy owing to stochastic scattering off turbulent fluctuations. We will discuss in detail this second acceleration stage in Section 7.

Zoom In Zoom Out Reset image size

By constructing the pdf of $| {J}_{z,p}| /{J}_{z,\mathrm{rms}}$ for the high-energy particles at injection and for all particles at a representative time (taken at ct/l = 2.5), we find results that are similar to the ones we have obtained for the 2D case. Figure 11(b) indeed shows that the pdf of the particles at injection (red circles) peaks at $| {J}_{z,p}| /{J}_{z,\mathrm{rms}}\sim 2.5$, as opposed to the pdf of the entire population of particles at the representative time ct/l = 2.5 (blue diamonds), which peaks at $| {J}_{z,p}| /{J}_{z,\mathrm{rms}}\sim 0$. Again, particles at injection feel a substantial electric current density in the direction of the mean magnetic field. The peak of the pdf for the particles at injection is at a lower value of $| {J}_{z,p}| /{J}_{z,\mathrm{rms}}$ than in 2D, and in general there are weaker $| {J}_{z,p}| /{J}_{z,\mathrm{rms}}$ wings for both the pdf of all particles and the pdf of particles experiencing injection. This can be attributed to the lower levels of intermittency that characterize 3D magnetized turbulence with respect to its 2D counterpart (e.g., Biskamp 2003). Nevertheless, about 80% of the particles are injected in regions with $| {J}_{z,p}| \geqslant 2\,{J}_{z,\mathrm{rms}}$. On the other hand, only approximately 11% of the entire population of particles (at the representative time ct/l = 2.5) reside at $| {J}_{z,p}| \geqslant 2\,{J}_{z,\mathrm{rms}}$. Therefore, also in 3D, special locations of high electric current density are associated with particle injection.

The spatial locations with $| {J}_{z}| \geqslant 2\,{J}_{z,\mathrm{rms}}$ are associated with current ribbons that are predominantly elongated along the mean magnetic field ${{\boldsymbol{B}}}_{0}$. In Figure 12, we show the morphology of these regions for two representative planes perpendicular to ${{\boldsymbol{B}}}_{0}$ (taken at ct/l = 2.5). These regions are sheet-like structures with a variety of length scales. We can see that the majority of the particles undergoing injection, whose location is shown by the red circles, resides at these current sheets. A large fraction of these current sheets are active reconnection layers, fragmenting into plasmoids. A typical example of such reconnecting current sheets is shown in Figure 13. We can see four flux ropes (3D plasmoids) that are formed within the current sheet (and elongated in the direction of the mean magnetic field), which is the typical signature of fast plasmoid-mediated reconnection. We will see in the next subsection that current sheets undergoing fast reconnection are important for having efficient particle injection, as they are capable to "process" a significant fraction of particles (from the thermal pool) during their lifetime in the turbulent plasma.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Reconnecting current sheets are a viable source of particle injection in typical astrophysical systems (${\ell } \ggg {d}_{e0}$) only if the injection efficiency (i.e., the fraction of particles going through the injection phase) is large and independent of system size. Here we show that this is indeed expected for our turbulence studies.

The rate at which a reconnecting current sheet can process particles is proportional to the normalized reconnection speed β_R = v_R/c, which essentially quantifies the speed of the reconnection process. This rate would be low for very elongated current sheets, as the large aspect ratio has the effect of throttling the reconnection rate. Indeed, a stable current sheet would be able to reach an asymptotic width determined by the microphysics of the plasma. For a collisionless relativistic pair plasma, the steady-state solution for the half-width of a reconnecting current sheet is (Comisso & Asenjo 2014)

$$\begin{eqnarray}&&{\lambda }_{\infty }\simeq {d}_{w}=\sqrt{\displaystyle \frac{{{mc}}^{2}}{4\pi {{ne}}^{2}}w},\end{eqnarray} \tag{ 12 }$$

where w = K₃(1/θ)/K₂(1/θ) is the enthalpy per particle in units of mc². For a thinning current sheet, ${\lambda }_{\infty }$ is the asymptotic limit of its half-width. For θ = k_BT/mc² ≫ 1, ${d}_{w}=\sqrt{{\gamma }_{{th}}{{mc}}^{2}/3\pi {{ne}}^{2}}\sim {d}_{e}$. Then, for a current sheet of half-length $\xi \gg {\lambda }_{\infty }$ and a compression ratio between inflow and outflow of order unity, the steady-state reconnection rate is

$$\begin{eqnarray}&&{\beta }_{R}\sim \displaystyle \frac{{d}_{w}}{\xi }\ll 1.\end{eqnarray} \tag{ 13 }$$

Since current sheets generated by outer-scale eddies (which, as we discuss below, are the ones that dominate particle injection) have half-length ξ ∼ ℓ larger than d_w by many orders of magnitude, the reconnection rate, as well as the injection efficiency, would be extremely low in this scenario.

However, plasmoids (which form copiously in our simulations) can break the reconnection layer into shorter elements, consequently leading to a regime of fast nonlinear reconnection (Daughton et al. 2006; Daughton & Karimabadi 2007; Bhattacharjee et al. 2009; Daughton et al. 2009; Huang & Bhattacharjee 2010; Uzdensky et al. 2010). This can happen if the plasmoids disrupt the current sheet within its characteristic lifetime, i.e., within one nonlinear eddy turnover time (Carbone et al. 1990; Boldyrev & Loureiro 2017; Loureiro & Boldyrev 2017; Mallet et al. 2017; Comisso et al. 2018; Dong et al. 2018; Walker et al. 2018). Fast magnetic reconnection essentially begins when plasmoids become nonlinear, namely, when the current density fluctuations caused by the growing plasmoids are of the same order as the current density of the reconnection layer (see Figure 4 in Huang et al. 2017). Therefore, understanding the plasmoid formation in the context of a forming current sheet is essential to understand the onset of fast magnetic reconnection and ensuing particle injection.

In order to evaluate the conditions for plasmoid formation and current sheet disruption, we need to analyze the growth rate of tearing (or "reconnecting") modes in such a current sheet. The collisionless tearing mode dispersion relation for a relativistic pair plasma can be obtained from the relativistic pair-plasma fluid equations (e.g., Koide 2009) by applying the standard tearing mode analysis (Furth et al. 1963; Coppi et al. 1976; Ara et al. 1978). In this way, one can obtain, for arbitrary values of the tearing stability parameter Δ' (Furth et al. 1963), the dispersion relation

$$\begin{eqnarray}&&{\gamma }^{1/2}{\tau }_{H}^{1/2}{\left(\displaystyle \frac{\lambda }{{d}_{w}}\right)}^{3/2}\displaystyle \frac{{\rm{\Gamma }}\left[({\rm{\Upsilon }}-1)/4\right]}{{\rm{\Gamma }}\left[({\rm{\Upsilon }}+5)/4\right]}=-\displaystyle \frac{8}{\pi }{\rm{\Delta }}^{\prime} ,\end{eqnarray} \tag{ 14 }$$

where

$$\begin{eqnarray}&&{\tau }_{H}=\displaystyle \frac{1}{{k}_{\xi }{v}_{A\lambda }},\,{\rm{\Upsilon }}=\gamma \,{\tau }_{H}\,\displaystyle \frac{\lambda }{{d}_{w}},\end{eqnarray} \tag{ 15 }$$

γ is the growth rate, k_ξ is the wavenumber in the ξ-direction, λ is the current sheet half-width, v_Aλ is the Alfvén speed based on the reconnecting magnetic field, and Γ(z) indicates the gamma function. This dispersion relation matches the nonrelativistic one (Porcelli 1991) when w → 1, i.e., when the plasma is cold.¹ Equation (14) can be further simplified for short-wavelength modes (small Δ') and long-wavelength modes (large Δ'), which is convenient in order to derive analytically the conditions for current sheet disruption. For ϒ ≪ 1, the small-Δ' regime, the growth rate, and the inner tearing layer half-width (where the ideal MHD approximation breaks down owing to the finite electron and positron inertia) of the instability are

$$\begin{eqnarray}&&{\gamma }_{s}={\left[\displaystyle \frac{{\rm{\Gamma }}\left(\tfrac{1}{4}\right)}{2\pi {\rm{\Gamma }}\left(\tfrac{3}{4}\right)}\right]}^{2}\displaystyle \frac{{\left({\rm{\Delta }}^{\prime} \lambda \right)}^{2}}{{\tau }_{H}}{\left(\displaystyle \frac{{d}_{w}}{\lambda }\right)}^{3},\qquad {\delta }_{\mathrm{in}}={d}_{w}^{2}{\rm{\Delta }}^{\prime} .\end{eqnarray} \tag{ 16 }$$

On the other hand, for ${\rm{\Upsilon }}\to {1}^{-}$, in the large-Δ' regime, the growth rate and the inner tearing layer half-width are

$$\begin{eqnarray}&&{\gamma }_{l}=\displaystyle \frac{1}{{\tau }_{H}}\left(\displaystyle \frac{{d}_{w}}{\lambda }\right),\qquad {\delta }_{\mathrm{in}}={d}_{w}.\end{eqnarray} \tag{ 17 }$$

Using these relations, together with the tearing stability index²

$$\begin{eqnarray}&&{\rm{\Delta }}^{\prime} =\displaystyle \frac{2}{\lambda }\left(\displaystyle \frac{1}{{k}_{\xi }\lambda }-{k}_{\xi }\lambda \right),\end{eqnarray} \tag{ 18 }$$

it can be shown that the dominant tearing mode at current sheet disruption scales like the fastest-growing mode (see Comisso et al. 2016, 2017; Huang et al. 2017), so that the instability wavenumber at current sheet disruption turns out to be simply

$$\begin{eqnarray}&&{k}_{\xi ,d}\sim \displaystyle \frac{{d}_{w}}{{\lambda }_{d}^{2}},\end{eqnarray} \tag{ 19 }$$

where the subscript "d" denotes current sheet disruption. This implies also that the growth rate and the inner tearing layer half-width at current sheet disruption are

$$\begin{eqnarray}&&{\gamma }_{d}\displaystyle \frac{\xi }{{v}_{A\lambda }}\sim \displaystyle \frac{{d}_{w}^{2}}{{\lambda }_{d}^{3}}\xi ,\qquad {\delta }_{\mathrm{in},d}\sim {d}_{w}.\end{eqnarray} \tag{ 20 }$$

From these expressions, one still needs to determine the current sheet half-width at disruption, λ_d, in order to know the wavenumber k_ξ,d and the growth rate γ_d. We calculate the width of the current sheet at disruption by using the principle of least time introduced in Comisso et al. (2016, 2017), substituting the resistive tearing mode dispersion relation with the collisionless dispersion relation discussed above. Then, for a rapid current sheet that forms on the Alfvénic timescale, the mode that becomes nonlinear in the shortest time disrupts the current sheet when

$$\begin{eqnarray}&&\displaystyle \frac{{\lambda }_{d}}{\xi }\,{\left[\mathrm{ln}\left(\displaystyle \frac{1}{2{\hat{\epsilon }}^{1/2}}{\left(\displaystyle \frac{{d}_{w}}{\xi }\right)}^{1+\alpha /2}{\left(\displaystyle \frac{\xi }{{\lambda }_{d}}\right)}^{1/2+\alpha }\right)\right]}^{1/3}\simeq {c}_{\lambda }{\left(\displaystyle \frac{{d}_{w}}{\xi }\right)}^{2/3},\end{eqnarray} \tag{ 21 }$$

where c_λ is an O(1) constant, $\hat{\epsilon }=\epsilon /(\delta {B}_{\lambda }\xi )$ is a normalized amplitude of the noise that seeds the instability (evaluated at the disruption scale), δB_λ is the characteristic magnetic field fluctuation at scale λ, and α is an index that depends on the spectrum of the noise, which is related to the turbulence spectrum as ${P}_{B}({k}_{\xi })\propto {{k}_{\xi }}^{1-2\alpha }$ (Comisso et al. 2018). Equation (21) can be solved exactly in terms of the Lambert W function, but here we prefer to consider an asymptotic solution that yields more transparent results. Therefore, we solve Equation (21) by iteration, obtaining, at the first order, the solution

$$\begin{eqnarray}&&{\lambda }_{d}\sim {d}_{w}^{2/3}{\xi }^{1/3}{\left[\mathrm{ln}\left(\displaystyle \frac{1}{2{\hat{\epsilon }}^{1/2}}{\left(\displaystyle \frac{{d}_{w}}{\xi }\right)}^{\tfrac{4-\alpha }{6}}\right)\right]}^{-1/3},\end{eqnarray} \tag{ 22 }$$

which gives us the critical current sheet width that determines the layer disruption and the onset of fast reconnection. Finally, the growth rate of the instability when the current sheet reaches this ratio is

$$\begin{eqnarray}&&{\gamma }_{d}\sim \displaystyle \frac{{v}_{A\lambda }}{\xi }\mathrm{ln}\left(\displaystyle \frac{1}{2{\hat{\epsilon }}^{1/2}}{\left(\displaystyle \frac{{d}_{w}}{\xi }\right)}^{\tfrac{4-\alpha }{6}}\right),\end{eqnarray} \tag{ 23 }$$

while the wavenumber of the dominant mode becomes

$$\begin{eqnarray}&&{k}_{\xi ,d}\sim {d}_{w}^{-1/3}{\xi }^{-2/3}{\left[\mathrm{ln}\left(\displaystyle \frac{1}{2{\hat{\epsilon }}^{1/2}}{\left(\displaystyle \frac{{d}_{w}}{\xi }\right)}^{\tfrac{4-\alpha }{6}}\right)\right]}^{2/3}.\end{eqnarray} \tag{ 24 }$$

We obtain from Equation (22) that ${\lambda }_{d}\gg {\lambda }_{\infty }\simeq {d}_{w}$ for outer-scale current sheets with ξ ∼ ℓ ≫ d_w, as is expected under typical astrophysical conditions. Therefore, an outer-scale current sheet disrupts in a chain of plasmoids before reaching the kinetic scale d_w, while interplasmoid layers, being shorter, can reach the thickness d_w. Equation (22) tells us also that current sheets disrupt at a larger thickness for larger noise levels. However, the dependence is only logarithmic. From the other two relations, Equation (23) and Equation (24), we have that the growth rate of the dominant mode is ${\gamma }_{d}\xi /{v}_{A\lambda }\gg 1$ at current sheet disruption, as is required for the instability to amplify the perturbation to a significant level within the lifetime of the current sheet (Comisso et al. 2018). Also, the number of plasmoids fragmenting the outer-scale current sheets, which is $\propto {k}_{\xi ,d}{\ell }$, increases as ${\ell }/{d}_{w}$ increases and the noise of the system decreases. As an example, we show in Figures 14 and 15 that a larger number of plasmoids form when the domain is increased by a factor of 4 with respect to the reference 2D simulation. In this simulation, as we argue below in this section, efficient plasmoid formation keeps the reconnection speed and the injection efficiency high when increasing system size. As a result, the fraction of nonthermal particles remains about the same when moving from the reference box size L/d_e0 = 1640 up to L/d_e0 = 6560 (see Figure 2(b) in Comisso & Sironi 2018).

Zoom In Zoom Out Reset image size

When the reconnection layer becomes dominated by the presence of plasmoids, soon after the condition λ ∼ λ_d is met, the complexity of the dynamics gives rise to a strongly time-dependent process. Nevertheless, in a statistical steady-state, we may expect that the reconnection layer containing the main X-point, which is the one that determines the global reconnection rate, has a bounded aspect ratio ξ_X/λ_X. If ξ_X/λ_X ∼ 1, the reconnection process would choke itself off, since this would imply β_R ∼ 0 (Comisso & Bhattacharjee 2016). This means that ξ_X/λ_X ≫ 1 in a steady reconnection process. On the other hand, the reconnection layer at the main X-point cannot be longer than the marginally stable sheet. Indeed, the fractal-like process of current sheet disruption due to the plasmoid instability terminates when the length of the innermost local current layer of the chain is shorter than the critical length ξ_c (Huang & Bhattacharjee 2010; Uzdensky et al. 2010; Comisso et al. 2015; Comisso & Grasso 2016). Therefore, ξ_X ≲ ξ_c is also expected. At present there are no analytical estimates for the aspect ratio ξ_c/λ_X, which might also depend on the noise level (e.g., Ni et al. 2010; Huang et al. 2017; Shi et al. 2018). However, numerical simulations have found ξ_c/λ_X ∼ 50 in the collisionless regime (e.g., Daughton et al. 2006; Daughton & Karimabadi 2007; Ji & Daughton 2011). As a consequence, for a compression ratio between inflow and outflow of order unity, the reconnection rate is bounded from above and below as

$$\begin{eqnarray}&&1/50\lesssim {\beta }_{R}\ll 1,\end{eqnarray} \tag{ 25 }$$

which classify it as a fast reconnection rate. More precisely, numerical simulations consistently indicate that β_R is an O(0.1) quantity (for relativistic pair plasmas, see, e.g., Zenitani et al. 2009; Bessho & Bhattacharjee 2012; Cerutti et al. 2012; Guo et al. 2014; Kagan et al. 2015; Liu et al. 2015, 2017; Sironi et al. 2016; Werner & Uzdensky 2017).

The aforementioned properties of reconnecting current sheets are important in regulating the particle injection efficiency. Here we show that the fraction of particles processed by reconnecting current sheets is independent of the system size and is quite large (despite the small filling fraction of current sheets) as long as the reconnection rate is high. To this purpose, let us consider a generic current sheet of characteristic length 2ξ and thickness 2λ, whose lifetime is approximately given by the local eddy turnover time ${\tau }_{\mathrm{nl}}\sim {\tau }_{A\xi }=\xi /{v}_{A\lambda }$, assuming critical balance (Goldreich & Sridhar 1995; Boldyrev 2006). If fast reconnection occurs for a time close to the eddy turnover time (see, e.g., Figure 15), a single reconnecting current sheet can "process" the upstream plasma up to a distance

$$\begin{eqnarray}&&{\lambda }_{R,j}={\beta }_{R,j}\,c\,{\tau }_{\mathrm{nl},j}\sim {\beta }_{R,j}\,c\displaystyle \frac{{\xi }_{j}}{{v}_{A\lambda ,j}},\end{eqnarray} \tag{ 26 }$$

where j labels the jth current sheet among the population of current sheets present at a given time, and the subscript "R" stands for reconnection. Since the surface processed by the entire population of reconnecting current sheets is in good approximation the one processed by the largest-scale ones, whose length scale corresponds to the turbulence integral length ℓ (e.g., Servidio et al. 2009), we have that magnetic reconnection can process a plasma surface

$$\begin{eqnarray}&&{{ \mathcal A }}_{R}={\sum }_{j}{\lambda }_{R,j}\,{\xi }_{j}\sim {\beta }_{R}{L}^{2}\end{eqnarray} \tag{ 27 }$$

in one large-eddy turnover time. Here we have used n_cs ∼ (L/ℓ)² as an estimate for the number of outer-scale current sheets, and β_R is the average reconnection rate. Furthermore, if we consider that current sheets in 3D are sheet-like structures with $2\lambda \ll 2\xi \lesssim 2{l}_{\parallel }$, with 2l_∥ indicating the direction along the magnetic field, we can obtain that, in one large-eddy turnover time, the reconnecting current sheets process a plasma volume

$$\begin{eqnarray}&&{{ \mathcal V }}_{R}={\sum }_{j}{\lambda }_{R,j}\,{\xi }_{j}\,{l}_{\parallel ,j}\sim {\beta }_{R}{L}^{3}.\end{eqnarray} \tag{ 28 }$$

Therefore, according to Equations (27)/(28), the plasma surface/volume processed by the reconnecting current sheets is a fixed fraction of the domain if β_R is independent of the system size, as discussed above. Moreover, since β_R is an O(0.1) quantity, magnetic reconnection can process large volumes of magnetic energy in few outer-scale eddy turnover times.

Zoom In Zoom Out Reset image size

In the next sections, we will address how particles are energized both in the injection phase and in the subsequent stochastic acceleration phase, and we will analyze the signatures of the acceleration process on the particle distribution.

The time evolution of the overall particle distribution with respect to $\cos \alpha$ is shown in Figure 19(a) for the reference 2D simulation and in Figure 19(b) for the reference 3D simulation. As turbulence evolves, the pitch-angle distribution becomes anisotropic with strong peaks at $\cos \alpha =\pm 1$, i.e., for particles moving along the magnetic field lines. Pronounced peaks of the pitch-angle distribution near $\cos \alpha =\pm 1$ have also been found in nonrelativistic plasma turbulence at low ${\beta }_{p}$ (e.g., Pecora et al. 2018), with ${\beta }_{p}={n}_{0}{k}_{B}T/(\langle {B}^{2}\rangle /8\pi )$ indicating the plasma beta, i.e., the ratio of thermal pressure to magnetic pressure. Indeed, the low-β_p regime is similar to the high-σ regime investigated here, in the sense that in both cases the magnetic energy density dominates over the thermal energy density (in our simulations the initial plasma beta is ${\beta }_{p}=2{\theta }_{0}/[{w}_{0}({\sigma }_{z}+{\sigma }_{0})]$). The pdf's of 2D and 3D simulations are similar; nevertheless, the 3D case exhibits higher probability peaks near $\cos \alpha =\pm 1$. Furthermore, the pitch-angle distribution evolves more rapidly in 3D as a consequence of the faster conversion of magnetic energy into particle energy. The large fraction of particles having velocity strongly aligned/antialigned with the local magnetic field is a natural expectation of injection mediated by magnetic reconnection, which can efficiently energize particles through the work done by ${{\boldsymbol{E}}}_{\parallel }$. As we have shown, reconnecting current sheets can process a large fraction of particles in just a few c/l (see Equations (27) and (28)).

Zoom In Zoom Out Reset image size

The pdf's illustrated in Figure 19 are dominated by low-energy particles (i.e., near the spectral peak), since they control the number census (see Figures 5 and 6). In order to characterize the anisotropy of particles of higher energy, we construct pdf's of $\cos \alpha$ for different populations of particles depending on their Lorentz factor. We collected particle data from a time range ${ct}/l\in [3,12]$ in order to increase statistics. However, we have verified that decreasing this range (up to a single time snapshot taken at late times) does not modify the results. These results are shown in Figures 20(a) and (b) for 2D and 3D turbulence, respectively. At very low energies (γ ∼ 1), the particle distribution remains nearly isotropic. These are the particles of our initial Maxwellian, which have not been energized. At moderate Lorentz factors (γ ∼ 15), the particle distribution displays strong peaks close to $\cos \alpha =\pm 1$, in analogy with the results shown in Figure 19. At higher Lorentz factors, the pitch-angle distribution evolves into a "butterfly distribution" with minima at both $\cos \alpha =\pm 1$ and $\cos \alpha =0$. This phenomenon occurs at Lorentz factors γ ∼ 50 for the simulations with σ₀ = 10 shown in Figure 19. At even higher energies (γ ≫ 50), the pitch-angle distribution becomes eventually peaked at $\cos \alpha =0$, i.e., for particles moving in the plane perpendicular to the local magnetic field. This trend can be displayed using a distribution $f\left(\cos \alpha ,\gamma \right)$ of particles with respect to $\cos \alpha$ and γ. This distribution, shown in Figures 20(c) and (d) for 2D and 3D turbulence, respectively, has been normalized such that

$$\begin{eqnarray}&&{\int }_{-1}^{1}f\left(\cos \alpha ,\gamma \right)d(\cos \alpha )=1.\end{eqnarray} \tag{ 33 }$$

In these plots, the peaks of $f\left(\cos \alpha ,\gamma \right)$ are located at $\cos \alpha =\pm 1$ for low energies, and then they move toward $\cos \alpha =0$ until γ ∼ 80. At higher energies, the peak of the distribution remains located at $\cos \alpha =0$, with particles that lie progressively more perpendicular to the local magnetic field as their energy increases.

Zoom In Zoom Out Reset image size

The energy-dependent anisotropy illustrated in Figure 20 reflects the different acceleration mechanisms that operate at different energies (see Section 5). At low energies, the contribution of the ${{\boldsymbol{v}}}_{\parallel }\cdot {{\boldsymbol{E}}}_{\parallel }$ energization is dominant, so that particles end up being strongly aligned/antialigned with the magnetic field ($\cos \alpha \sim \pm 1$). On the other hand, as the energy increases, the ${{\boldsymbol{v}}}_{\perp }\cdot {{\boldsymbol{E}}}_{\perp }$ energization takes over and propels the particles in the direction perpendicular to the local magnetic field. The timescale of this acceleration is fast compared to the pitch-angle scattering timescale, so that particles retain their orientation $\cos \alpha \sim 0$ for long times.

The results shown in Figure 20 for magnetization σ₀ = 10 also hold for the other magnetizations we investigate. In Figure 21, we present the results from four simulations that differ in magnetization ${\sigma }_{0}\in \left\{5,10,20,40\right\}$. Here we show only the results from 3D simulations, since those from 2D simulations are analogous. The ranges in γ are scaled with σ₀, which provides the typical energy scale (e.g., the starting point of the high-energy nonthermal tail, γ_st, increases linearly with σ₀, as illustrated by Equation (8)). For $\gamma \sim ({\sigma }_{0}/2){\gamma }_{{th}0}$ (solid lines), we have a pitch-angle distribution peaked at $\cos \alpha \sim \pm 1$ (the only difference is that the percentage of particles aligned/antialigned with the local magnetic field slightly increases with σ₀). The butterfly distribution with minima at $\cos \alpha =\pm 1,0$ appears for $\gamma \sim 5({\sigma }_{0}/2){\gamma }_{{th}0}$ (long-dashed lines). Finally, for $\gamma \gg 5({\sigma }_{0}/2){\gamma }_{{th}0}$, well into the nonthermal tail, the particle velocities become mostly perpendicular to the magnetic field, and we can see that all the distributions are peaked at $\cos \alpha =0$ (see dashed lines).

Zoom In Zoom Out Reset image size

The results on the anisotropy of the pitch-angle distributions, computed with respect to the local magnetic field ${\boldsymbol{B}}={{\boldsymbol{B}}}_{0}+\delta {\boldsymbol{B}}$, suggest that the four-velocity distribution function, with respect to the mean magnetic field ${{\boldsymbol{B}}}_{0}={B}_{0}\hat{{\boldsymbol{z}}}$, should also display significant anisotropy. Indeed, as the turbulence fluctuations decay, the local magnetic field becomes progressively more aligned with the direction of the mean magnetic field.

We calculated the domain-averaged four-velocity distributions in the x-y plane, $f\left(\gamma {\beta }_{x},\gamma {\beta }_{y}\right)$, and in the x-z plane, $f\left(\gamma {\beta }_{x},\gamma {\beta }_{z}\right)$, from the same samples of particles used to analyze the local pitch-angle distributions. The results for our reference 3D simulation (the 2D case is analogous) are shown in Figure 22 (and a zoom-in in Figure 23). As for Figure 20, we collected particle data in the time range ${ct}/l\in [3,12]$ in order to increase statistics, but we have also verified that decreasing this range (up to a single time snapshot taken at late times) does not modify the results. As we expected, from Figure 22 we find that the four-velocity distribution is isotropic in the plane perpendicular to ${{\boldsymbol{B}}}_{0}$ (top panel), while it develops more complex features with respect to planes that contain ${{\boldsymbol{B}}}_{0}$, as for the case of $f\left(\gamma {\beta }_{x},\gamma {\beta }_{z}\right)$ (bottom panel). The results are analogous when considering $f\left(\gamma {\beta }_{y},\gamma {\beta }_{z}\right)$ or $f\left({({\gamma }^{2}{\beta }_{x}^{2}+{\gamma }^{2}{\beta }_{y}^{2})}^{1/2},\gamma {\beta }_{z}\right)$. The distribution $f\left(\gamma {\beta }_{x},\gamma {\beta }_{z}\right)$ displays a core region elongated in the $\gamma {\beta }_{z}$ direction, as particle velocities are mostly aligned/antialigned with the magnetic field at low energies. Furthermore, a close inspection shows that there is an intermediate-energy region, with the majority of particles residing in a double cone whose axis is the direction of the mean magnetic field ${{\boldsymbol{B}}}_{0}$ (see Figure 23). This is the intermediate-energy range, in which the peak of pitch-angle distribution moves from $\cos \alpha =\pm 1$ toward $\cos \alpha =0$. At even higher energies, Figure 22 shows that $f\left(\gamma {\beta }_{y},\gamma {\beta }_{z}\right)$ becomes elongated in the direction perpendicular to ${{\boldsymbol{B}}}_{0}$, consistently with the dominance of ${{\boldsymbol{v}}}_{\perp }\cdot {{\boldsymbol{E}}}_{\perp }$ energization at higher energies and the resulting anisotropy of the pitch-angle cosine.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We show that while the qualitative and quantitative features of the pitch-angle distributions are similar in our 2D and 3D simulations, the turbulent mixing (in space) of the energized particles is quite different. Particle mixing in 2D is expected to be less efficient than in 3D, since the translation-invariant symmetry along $\hat{{\boldsymbol{z}}}$ seriously constrains the 2D dynamics. As a consequence, regions of space devoid of high-energy particles can be retained for a larger number of outer-scale eddy turnover times in 2D simulations.

In Figure 24, we show how the energized particles are distributed in the spatial domain in 2D (left column) and 3D (right column). For both cases, we plot the cell-averaged kinetic energy per particle, $\langle \gamma -1{\rangle }_{\mathrm{cell}}$, at three different times (from top to bottom). The mean kinetic energy is normalized by mc². The time snapshots are different for 2D (ct/l = 4.6, 7.7, 10.8) and 3D (ct/l = 2.7, 4.5, 6.3), to account for the faster turbulence decay in 3D. In both cases, the initial energization occurs at current sheets, which display high values of $\langle \gamma -1{\rangle }_{\mathrm{cell}}$, and then particles propagate outside current sheets in other regions of the domain (see top panels of Figure 24). As time progresses, energized particles diffuse in the spatial domain, and the mean kinetic energy per particle becomes more uniform (middle panels of Figure 24). However, in 2D the cores of the large-scale flux tubes remain essentially unaffected, as these overdense regions with n ≫ n₀ and with higher fluctuation magnetic energy density $\delta {B}^{2}/8\pi$ (see Figure 1) are mainly populated by low-energy particles that have not been processed by reconnecting current sheets. On the other hand, the 3D domain does not present such isolated regions of low $\langle \gamma -1{\rangle }_{\mathrm{cell}}$. At quite early times, the mean kinetic energy per particle becomes fairly homogeneous across the entire 3D domain, whereas the 2D simulation preserved regions of low $\langle \gamma -1{\rangle }_{\mathrm{cell}}$ for much longer (bottom panels of Figure 24). This different behavior between 2D and 3D turbulence is also reflected in the particle energy spectrum, with 2D turbulence retaining more particles with γ ≲ γ_th0 until late times (see Figures 5 and 6).

Zoom In Zoom Out Reset image size

Particles that are stochastically scattered off the turbulent fluctuations experience a biased random walk in momentum space, which can be modeled with a Fokker-Planck approach (e.g., Blandford & Eichler 1987), provided that the fractional momentum change in a single scattering is sufficiently small. In this case, one could describe the process of stochastic acceleration from the point of view of a Fokker-Planck equation in energy space (e.g., Ramaty 1979)

$$\begin{eqnarray}&&\displaystyle \frac{\partial N}{\partial t}=-\displaystyle \frac{\partial }{\partial \gamma }\left({A}_{\gamma }N\right)+\displaystyle \frac{{\partial }^{2}}{\partial {\gamma }^{2}}\left({D}_{\gamma }N\right).\end{eqnarray} \tag{ 34 }$$

Here, as usual, N is the particle spectrum differential in energy, A_γ is the energy convection coefficient, and D_γ is the energy diffusion coefficient. Note that, in general, the convection and diffusion coefficients are time dependent (this is indeed the case for the turbulence simulations performed here). The convection coefficient A_γ represents the mean energy gain due to stochastic acceleration and is related to the diffusion coefficient in energy space as

$$\begin{eqnarray}&&{A}_{\gamma }=\displaystyle \frac{d\langle \gamma \rangle }{{dt}}=\displaystyle \frac{1}{{\gamma }^{2}}\displaystyle \frac{\partial }{\partial \gamma }\left({\gamma }^{2}{D}_{\gamma }\right).\end{eqnarray} \tag{ 35 }$$

The diffusion coefficient in energy space D_γ is also related to the diffusion coefficient in momentum space D_p, with D_γ ≃ D_p for the ultrarelativistic particles considered here. Given the fact that high-energy particles preferentially lie in the plane perpendicular to the mean field (Section 6) and that their energization is mostly contributed by perpendicular electric fields (Section 5), the momentum diffusion coefficient D_p is essentially identical to ${D}_{{p}_{\perp }}$, i.e., to the diffusion coefficient of momenta perpendicular to the mean field. The determination of this coefficient, or equivalently of the energy diffusion coefficient, establishes the properties of the stochastic acceleration phase.

We evaluate the energy diffusion coefficient directly from PIC simulations (see also Wong et al. 2019). To this aim, from each of the 2D and 3D simulations employed for this analysis, we tracked in time the positions, four-velocities, and electromagnetic field values of about 10⁷ particles that were randomly selected at the beginning of the simulation. From the time history of the particle evolution, we calculate the mean-square γ-variation

$$\begin{eqnarray}&&\langle {({\rm{\Delta }}\gamma )}^{2}\rangle =\displaystyle \frac{1}{{N}_{p}}\displaystyle \sum _{n=1}^{{N}_{p}}{\left({\gamma }_{n}(t)-{\gamma }_{n}({t}_{* })\right)}^{2}\end{eqnarray} \tag{ 36 }$$

for particles grouped in such a way that at an initial time t_* they belong to the same energy bin (N_p is the number of particles in the selected bin). The energy bin at t_* is chosen according to the particle energy calculated in the frame comoving with the drift velocity ${{\boldsymbol{v}}}_{D}=c{\boldsymbol{E}}\times {\boldsymbol{B}}/{B}^{2}$. For each particle, we perform a Lorentz boost from the observer/simulation frame to the local ${\boldsymbol{E}}\times {\boldsymbol{B}}$ frame, which results in the boosted Lorentz factor $\gamma ^{\prime} ={\gamma }_{D}\gamma \left(1-{{\boldsymbol{v}}}_{D}\cdot {\boldsymbol{v}}/{c}^{2}\right)$, where ${\gamma }_{D}=1/\sqrt{1-{({{\boldsymbol{v}}}_{D}/c)}^{2}}$ is the Lorentz factor for the drift velocity. Then, we evaluate Equation (36) by selecting particles in a small energy bin with $\gamma ^{\prime} \in [{\gamma }_{* }/\nu ,{\gamma }_{* }\nu ]$, where γ_* is the characteristic Lorentz factor of the energy bin and ν is a constant factor that should be close to unity (we choose ν = 1.1). Finally, the diffusion coefficient in energy space can be calculated as

$$\begin{eqnarray}&&{D}_{\gamma }=\displaystyle \frac{\langle {({\rm{\Delta }}\gamma )}^{2}\rangle }{2{\rm{\Delta }}t},\end{eqnarray} \tag{ 37 }$$

where ${\rm{\Delta }}t=t-{t}_{* }$ is a time interval that should be (i) long enough that the initial conditions become insignificant and particles are in the diffusive regime and (ii) short enough that the turbulence properties have not significantly changed. By using a large sample of particles in each energy bin, nonsecular variations of the particle energy are averaged out and the mean energy gain can be obtained.

The results of our analysis of the particle energy diffusion are reported in Figures 25 and 26, for 2D and 3D simulations, respectively. The top panels show the mean-square variation $\langle {({\rm{\Delta }}\gamma )}^{2}\rangle$ for particles binned according to their initial energy at ${{ct}}_{* }/l=5.25$ for the reference 2D simulation and ct_*/l = 3 for the reference 3D simulation. In both cases, at the selected time t_*, turbulence is well developed and the time-dependent magnetization calculated with the magnetic energy in turbulent fields is $\sigma ({t}_{* })\sim 1$. The plots indicate that a diffusive behavior in energy space, $\langle {({\rm{\Delta }}\gamma )}^{2}\rangle \propto {\rm{\Delta }}t$ (compare with dashed black lines), is achieved after cΔt/l ∼ 1, in both 2D and 3D reference simulations. For shorter time intervals, particles preserve memory of the initial conditions and their motion is not diffusive. The slope at late times (dashed lines) depends on particle energy, and it allows us to quantify the energy dependence of the diffusion coefficient.

Zoom In Zoom Out Reset image size

The bottom panels of Figures 25 and 26 show the particle energy dependence of the energy diffusion coefficient from simulations with different initial magnetization σ₀ (indicated with different colors in the figures). The diffusion coefficient is evaluated using Equation (37) in the time interval cΔt/l = 1.875, starting from ct_*/l = 5.25 for 2D simulations and ct_*/l = 3 for 3D. We verified that the energy dependence remains the same when taking different time intervals, or by fitting the slopes of $\langle {({\rm{\Delta }}\gamma )}^{2}\rangle$ as a function of time in the diffusive regime (as done with the dashed lines in the top panels). In order to properly compare different σ₀, we display the energy diffusion coefficient as a function of the Lorentz factor normalized by γ_σ (see Equation (8)). The energy range where stochastic acceleration occurs starts at the beginning of the power-law high-energy tail of the particle spectrum, i.e., for γ/γ_σ ≳ 1. In the stochastic acceleration range, the energy diffusion coefficient scales as D_γ ∝ γ² (compare with the dashed black lines in the bottom panels). A similar dependence on the particle energy was also found in Lynn et al. (2014), Kimura et al. (2016, 2019), and Wong et al. (2019) and is consistent with particle acceleration by nonresonant and/or broadened resonant interactions with the turbulent fluctuations (e.g., Skilling 1975; Blandford & Eichler 1987; Schlickeiser 1989; Chandran 2000; Cho & Lazarian 2006; Lemoine 2019). Then, at higher energies, near the high-energy cutoff of the power law, the energy dependence of D_γ becomes weaker as the particle Larmor radius gets closer to the energy-containing scale of the turbulence.

Zoom In Zoom Out Reset image size

The energy diffusion coefficient also depends on the actual magnetization $\sigma ({t}_{* })=\langle \delta {B}^{2}\rangle /4\pi {n}_{0}{{wmc}}^{2}$. In order to better understand this dependence, in Figure 27 we plot the energy diffusion coefficient as a function of the magnetization σ at four different times t_* (in the range ${{ct}}_{* }/l\in [4,6]$ for 2D and ${{ct}}_{* }/l\in [2,4]$ for 3D) for the four simulations having different initial magnetization. Both 2D and 3D simulations show a clear trend of increasing diffusion coefficient with increasing magnetization. The 3D simulations are well fitted by a linear relation in σ (compare with dashed black line),

$$\begin{eqnarray}&&{D}_{\gamma }\sim 0.1\sigma \left(\displaystyle \frac{c}{l}\right)\,{\gamma }^{2}.\end{eqnarray} \tag{ 38 }$$

This scaling can be understood by noting that for a stochastic process akin to the original Fermi mechanism (e.g., Blandford & Eichler 1987; Lemoine 2019), the energy diffusion coefficient is

$$\begin{eqnarray}&&{D}_{\gamma }=\displaystyle \frac{1}{3}\langle {\gamma }_{V}^{2}{\beta }_{V}^{2}\rangle \displaystyle \frac{c}{{\lambda }_{\mathrm{mfp}}}{\gamma }^{2},\end{eqnarray} \tag{ 39 }$$

where $\langle {\gamma }_{V}^{2}{\beta }_{V}^{2}{\rangle }^{1/2}$ is the typical four-velocity of the scatterers and λ_mfp is the particle scattering mean free path. Therefore, if we estimate the scattering mean free path as ${\lambda }_{\mathrm{mfp}}\sim {({B}_{0}/\delta {B}_{\mathrm{rms}})}^{2}l$ and identify γ_V β_V with the dimensionless Alfvénic four-velocity, $\langle {\gamma }_{V}^{2}{\beta }_{V}^{2}\rangle \sim \langle {B}^{2}\rangle /4\pi {{nwmc}}^{2}$, from Equation (39) we obtain D_γ ∝ σ for $\langle {B}^{2}\rangle /{B}_{0}^{2}\sim 1$, in agreement with Equation (38). Then, from these results we can also estimate the stochastic acceleration timescale

$$\begin{eqnarray}&&{t}_{\mathrm{acc}}={\left|\displaystyle \frac{1}{\gamma }\displaystyle \frac{d\langle \gamma \rangle }{{dt}}\right|}^{-1}\sim \displaystyle \frac{3}{\sigma }\displaystyle \frac{l}{c}.\end{eqnarray} \tag{ 40 }$$

In our simulations, the acceleration timescale increases in time since σ decreases in time as a combined effect of the decaying turbulent fluctuations δB_rms(t) and the increase of the enthalpy per particle mc²w(t).

Zoom In Zoom Out Reset image size

As discussed in Sections 4 and 5, a large fraction of particles are preaccelerated by magnetic reconnection before being accelerated by scattering off the turbulent fluctuations. This two-stage acceleration process is shown in Figure 28 for both 2D and 3D simulations. Here, each colored curve represents the average Lorentz factor of particles having the same injection time t_inj (within Δt_inj = 0.32c/l for 2D and Δt_inj = 0.22c/l for 3D). The linear growth from $\langle \gamma \rangle \sim 1$ up to $\langle \gamma \rangle \sim 30$ (i.e., the injection phase) is powered by field-aligned electric fields, whose magnitude is $| {E}_{\parallel }| \simeq {\beta }_{R}\delta {B}_{\mathrm{rms}}$, via

$$\begin{eqnarray}&&\displaystyle \frac{d\langle \gamma \rangle }{{dt}}=\displaystyle \frac{e}{{mc}}{\beta }_{R}\delta {B}_{\mathrm{rms}}.\end{eqnarray} \tag{ 41 }$$

The dashed black lines in Figure 28 show Equation (41) taking a reconnection rate β_R ≃ 0.05, as appropriate for relativistic reconnection with guide field comparable to the alternating fields (Werner & Uzdensky 2017). After this first acceleration phase, stochastic acceleration takes place, and, as discussed above, we can estimate

$$\begin{eqnarray}&&\displaystyle \frac{d\langle \gamma \rangle }{{dt}}=4{\kappa }_{\mathrm{stoc}}\sigma \left(\displaystyle \frac{c}{l}\right)\,\gamma ,\end{eqnarray} \tag{ 42 }$$

with κ_stoc ∼ 0.03 from the 2D simulations and κ_stoc ∼ 0.1 from the 3D ones. Taking the temporal decay of the magnetic fluctuations, as well as the temporal increase of the relativistic enthalpy, directly from our simulations, we obtain the dotted–dashed lines shown in Figure 28. For the 3D case, the decrease in time of the stochastic acceleration rate is more pronounced than the 2D case as a consequence of the faster magnetic energy decay and the corresponding decrease of the magnetization σ.

Zoom In Zoom Out Reset image size

A final remark concerns the acceleration timescales associated with magnetic reconnection and turbulence fluctuations. Fast magnetic reconnection leads to the acceleration timescale t_acc = β_R⁻¹ (ρ_L/c), where ρ_L is the particle Larmor radius. On the other hand, we have seen that stochastic acceleration by turbulent fluctuation yields t_acc = (3/σ) (l/c). Therefore, for the hypothetical case in which reconnection could drive particles up to the highest energies (ρ_L ∼ l), the acceleration timescale of fast magnetic reconnection could actually be longer than the one associated with the turbulence fluctuations for σ ≳ 1. Indeed, in this magnetically dominated regime, turbulence provides an exceptionally fast acceleration mechanism that can potentially explain the most extreme astrophysical accelerators.