Low-energy Injection and Nonthermal Particle Acceleration in Relativistic Magnetic Turbulence

Magnetic turbulence in plasmas reveals itself through fluctuating magnetic fields, bulk velocity, and density over a broad range of spatial and temporal scales. It is commonly found and studied in astrophysical environments such as pulsar wind nebulae (Porth et al., 2014; Lyutikov et al., 2019; Cerutti & Giacinti, 2020; Lu et al., 2021), stellar coronae and flares (Matthaeus et al., 1999; Cranmer et al., 2007; Liu et al., 2006; Fu et al., 2020; Pongkitiwanichakul et al., 2021), black hole accretion disks (Balbus & Hawley, 1998; Brandenburg & Subramanian, 2005; Sun & Bai, 2021), radio lobes (Vogt & Enßlin, 2005; O’Sullivan et al., 2009), and jets from active galactic nuclei (Marscher et al., 2008; Zhang et al., 2023). All of these systems exhibit high-energy emissions that suggest nonthermal particle acceleration. In turbulent plasmas, the kinetic energy from large-scale motion cascades to smaller and smaller scales, which is eventually dissipated through turbulence-particle interactions. Understanding how particles in turbulent plasmas get accelerated to high energy is an unsolved problem in high-energy astrophysics.

Turbulence is often invoked as a particle acceleration mechanism that leads to nonthermal particle spectra. Recently, several studies have used kinetic particle-in-cell (PIC) simulations to gain insight into nonthermal particle acceleration mechanisms in its relativistic regime (Zhdankin et al., 2017, 2018; Comisso & Sironi, 2018, 2019; Wong et al., 2020; Hankla et al., 2021; Vega et al., 2022). The most commonly discussed acceleration mechanism in magnetic turbulence is stochastic Fermi acceleration (Fermi, 1949; Petrosian, 2012; Lemoine & Malkov, 2020), where particles can gain energy by scattering back and forth in the turbulent fluctuations. Magnetic reconnection (Biskamp, 2000; Zweibel & Yamada, 2009; Yamada et al., 2010; Ji et al., 2022; Yamada, 2022), which occurs naturally as magnetic turbulence generates thin current sheets, may also support strong particle acceleration (Sironi & Spitkovsky, 2014; Guo et al., 2014, 2015; Werner et al., 2016; Guo et al., 2020). More interestingly, magnetic reconnection can have an intriguing relation with turbulence and their interplay during particle acceleration is not completely clear (Loureiro & Boldyrev, 2017; Dong et al., 2018, 2022; Comisso & Sironi, 2019; Li et al., 2019; Zhang et al., 2021, 2024a; Guo et al., 2021). Nevertheless, these recent numerical simulations and theoretical models suggest that magnetic turbulence, especially in its relativistic limit ($\sigma\equiv B^{2}/4\pi h\gg 1$; i.e. the magnetic enthalpy $B^{2}/4\pi$ greatly exceeds the plasma enthalpy $h$), plays a major role in nonthermal particle acceleration.

In general, Fermi acceleration requires particle injection mechanism(s) to accelerate particles to energies that enable them to participate in a continual acceleration process. This process naturally defines an injection energy, beyond which injected particles enter the power-law range of the particle spectrum (French et al., 2023). The injection problem has recently been studied in the context of relativistic magnetic reconnection (Guo et al., 2019; Ball et al., 2019; Kilian et al., 2020; Sironi, 2022; French et al., 2023; Guo et al., 2023). While it has been suggested that during magnetic reconnection the parallel electric field $\textbf{E}_{\parallel}\equiv(\textbf{E}\cdot\textbf{B})\textbf{B}/|\textbf{B}|^{2}$ dominates the injection (Ball et al., 2019), studies have shown that perpendicular electric fields ($\textbf{E}_{\perp}\equiv\textbf{E}-\textbf{E}_{\parallel}$) can play an important, if not dominant role (Kilian et al., 2020; French et al., 2023). Meanwhile, X-points with $|E|>|B|$ are shown to be negligible for particle injection and high-energy acceleration (Guo et al., 2019, 2023). Particle injection has also been investigated in relativistic magnetic turbulence (Comisso & Sironi, 2019), where parallel electric fields in reconnection diffusion regions were concluded to dominate the injection process. Meanwhile, the subsequent particle energization in the power law was shown to be dominated by perpendicular electric fields ($\textbf{E}_{\perp}$) from stochastic scattering off turbulent fluctuations. However, Comisso & Sironi (2019) focused only on a small population of high energy particles with final energies many times greater than the injection energy. Since the importance of $\textbf{E}_{\perp}$ has been demonstrated in magnetic reconnection, it is worthwhile to investigate whether $\textbf{E}_{\perp}$ is important in magnetic turbulence as well.

In a recent study, French et al. (2023) analyzed particle injection and further acceleration in relativistic magnetic reconnection with emphasis on the influence of guide field and domain size. They measured the injection energy of each nonthermal particle spectrum using a spectral fitting procedure. They decompose the work done by parallel and perpendicular electric field components and quantify the contributions by different mechanisms, thereby illuminating which mechanism dominates the initial energization and the subsequent nonthermal acceleration. In this paper, we employ a similar methodology to study collisionless relativistic turbulence by carrying out two-dimensional (2D) PIC simulations and calculating the shares of work done by parallel ($W_{\parallel}$) and perpendicular ($W_{\perp}$) electric fields. We find that, similar to magnetic reconnection, the contribution of $W_{\perp}$ to particle injection grows with increasing domain size until the largest simulation domain, and may all exceed $50\%$ contribution for macroscale systems. However, in contrast to magnetic reconnection, the relative contributions of $W_{\parallel}$ vs $W_{\perp}$ to subsequent energization of particles of energies $\varepsilon>\varepsilon_{\rm inj}$ is relatively insensitive to domain size.

The rest of the paper is organized as follows: Section 2 describes our simulation setup. In Section 3 we present the simulation results and analyses for understanding the particle injection and nonthermal particle acceleration. Section 4 discusses implications for observations and summarize the conclusions.

We use the Vectorized Particle-In-Cell (VPIC) simulation code to investigate nonthermal particle acceleration in relativistic magnetic turbulence. VPIC solves the relativistic Maxwell-Vlasov equations to self-consistently evolve kinetic plasmas and their interaction with electromagnetic fields (Bowers et al., 2008a, b, 2009). We simulate magnetically-dominated decaying turbulence in a two-dimensional (2D) square domain ($x$-$y$) of size $L^{2}$. The initial setup is similar to earlier work (Comisso & Sironi, 2019; Pongkitiwanichakul et al., 2021; Zhang et al., 2023), where an electron-positron pair plasma is initialized with a turbulent magnetic field $\bm{B}=B_{0}\bm{\hat{z}}+\delta\bm{B}$. $B_{0}$ is the magnitude of the uniform component and $\delta\bm{B}$ is the fluctuating component, which is given by

Here, $\delta B(\bm{k})$ is the Fourier amplitude of the mode with wavevector $\bm{k}$, $\bm{\hat{\xi}}(\bm{k})=i\bm{k}\times\bm{B}_{0}/|\bm{k}\times\bm{B}_{0}|$ are the Alfvénic polarization unit vectors, and $\phi_{\bm{k}}$ expresses random phases. $\bm{k}$ represents the wavevector such that $\bm{k}=(k_{x},k_{y})$, where $k_{x}=2m\pi/L$ and $k_{y}=2n\pi/L$ with $m\in\{-N,\dots,-1,1,\dots,N\}$ and $n\in\{-N,\dots,-1,1,\dots,N\}$. $N$ is the number of modes along each dimension, which is set to be $8$ in this paper. We also define wavenumber $k=|\bm{k}|=\sqrt{k_{x}^{2}+k_{y}^{2}}$ as the amplitude of the wavevector. The boundary conditions are periodic for both particles and fields. The initial electric field E is set to 0.

We initialize the plasma and magnetic fields with magnetization parameter $\sigma_{0}\equiv B_{0}^{2}/(4\pi n_{0}m_{e}c^{2})=\omega_{\rm ce}^{2}/2\omega_{\rm pe}^{2}=20$, where $\omega_{\rm pe}\equiv\sqrt{4\pi n_{e}e^{2}/m_{e}}$ is the plasma electron frequency and $\omega_{\rm ce}\equiv eB_{0}/m_{e}c$ is the electron cyclotron frequency defined using the uniform background magnetic field $B_{0}$. Here, $m_{e}$ is the electron mass, $c$ is the speed of light, $e$ is the electron charge, and $n_{0}=n_{p}+n_{e}$ is the number density of the pair plasma in the simulation domain. The turbulence amplitude ${\delta B}_{\rm rms0}/B_{0}=1$, where ${\delta B}_{\rm rms0}$ is the space-averaged root-mean-square value of the initial magnetic field fluctuations. The domain size $L$ is normalized by the electron skin depth $d_{e}\equiv c/\omega_{\rm pe}$ and each $d_{e}$ is resolved to 4 grid cells (i.e., $d_{e}=4\Delta x$). To allow most of the turbulent magnetic energy to be converted to the particles, the simulations are run for two light crossing times $2L/c$. To independently examine the influence of domain size on our results, we run an array of otherwise identical simulations with $L/d_{e}\in\{512,1024,1440,2048,2880,4096\}$.

In all our simulations, we use 100 particles of each species per cell that are initialized with a Maxwellian distribution with dimensionless temperature $\theta_{0}\equiv k_{B}T_{0}/m_{e}c^{2}=0.3$. Here, $k_{B}$ is the Boltzmann constant and $T_{0}$ is the initial plasma temperature. We also have done some test simulations with a larger number of particles per cell and/or higher spatial resolution and found that the results described below still hold.

For each simulation, we trace $\sim 200,000$ particles of each species and save the electric and magnetic fields $\bm{E}$ and $\bm{B}$ as well as velocities v at their positions at every time step, to understand their injection and nonthermal particle acceleration (Li et al., 2023).

Figure 1 shows the evolution of the magnitude of electric current density $|J/J_{0}|$ in the simulation domain for the simulation with $L/d_{e}=1440$ at times $\omega_{\rm pe}t=$ (a) $20$, (b) $200$, (c) $960$, and (d) $2880$, normalized to $J_{0}\equiv n_{0}ec/2$. The initial perturbation seen in panel (a) generates fluctuations across different scales, after a brief initial phase. As turbulence develops, many plasmoids¹¹1Note that many of the large plasmoids are due to the initial evolution of the initial perturbation, whereas during the evolution of the simulation small-scale plasmoids are generated during the reconnection process, which indicates that the energy is transferred to smaller scales. and current sheets are produced in 2D turbulence, where magnetic reconnection is likely to happen (panel b).

Figure 2 zooms in on a reconnection site occurring in the simulation at $\omega_{\rm pe}t=960$ and displays colormaps of (a-b) the absolute current density $\lvert J/J_{0}\rvert$, (c-d) the parallel electric field $E_{\parallel}$, and (e-f) the perpendicular electric field $E_{\perp}$. Here, $E_{\parallel}$ and $E_{\perp}$ are plotted in units of $B_{0}/\sqrt{2\sigma_{0}}$. From inspecting these figures we see that $E_{\perp}\gg E_{\parallel}$ on a global scale, and it becomes clear that $E_{\parallel}$ is well-localized to reconnection X-points at plasmoid interfaces. However, $E_{\perp}$ can still have a substantial strength at reconnection regions owing to the reconnection outflow immediately downstream of these X-points (French et al., 2023).

In Figure 3, we show how the fractions of energy stored in particles, magnetic fields, and electric fields evolve as the simulation proceeds. The total energy is well conserved. As the turbulence decays and reconnection events begin liberating magnetic field energy into nearby particles, the fraction of energy stored by particles grows from $\sim 2.5\%$ at $t=0$ to $\sim 35\%$ by the final time. This corresponds to the decrease of magnetic field energy. Since the initial electric field is set to be zero and induced rapidly due to the changing magnetic field, its energy experience a strong, transit growth in the initial stage $\omega_{\rm pe}t<500$.

Figure 4 shows the power spectra of magnetic field fluctuations $\delta\bm{B}$ for various domain sizes at 2 light crossing time. The power $P(k)$ is normalized by the total power for that simulation at that time. In all the cases, we observe that a Kolmogorov-like $k^{-5/3}$ scaling quickly established and last until the end of the simulation. For larger domains, the fluctuations extends to larger spatial scales (lower $k$), and the small scale fluctuations have lower amplitude. Meanwhile, the amplitude of the fluctuation $\delta B_{rms}$ decays from 1.0 to about 0.5 in the end of the simulation, quite consistently in all simulations.

We analyze the nonthermal spectra for all of our simulations, and quantify several spectral features: power-law index $p$ that represents the slope in the nonthermal region of the spectrum, the injection energy $\varepsilon_{\rm inj}$ and cut-off energy $\varepsilon_{c}$ that mark the lower and upper energy bounds of the nonthermal region respectively, and the power-law extent $R\equiv\varepsilon_{c}/\varepsilon_{\rm inj}$. From these nonthermal particle spectra, we perform a fitting procedure at the end of the simulation to obtain the characteristic parameters ($\varepsilon_{\rm inj}$, $\varepsilon_{c}$, $p$) of our particle spectra (Werner et al., 2017; French et al., 2023), from which we also calculate the power-law extent $R$. The procedure begins by smoothing a particle spectrum $f$ via isotonic regression so that the local power-law index $p_{\varepsilon}\equiv-\,d\log{f(\varepsilon)}/d\log{\varepsilon}$ can be defined. Here, $\varepsilon$ refers to the particle energy. Then all “valid” power-law segments are obtained by brute force, where validity is determined by a predefined power-law tolerance and minimum power-law extent, yielding a list of power-law indices, injection energies, and cutoff energies (see French et al. (2023) for details). Finally, after removing duplicates (e.g., identical power-law segments resulting from different $p$-tolerances) and outliers (i.e., data points beyond $\pm$ 2 standard deviations from the mean) from each collection of values, each characteristic parameter ($p$, $\varepsilon_{\rm inj}$, $\varepsilon_{c}$) is defined by the mean of its collection and its error by one standard deviation of its collection.

Figure 5(a) shows the time evolution of particle energy spectra for the simulation with domain size $L/d_{e}=1440$. As the simulation starts, the turbulent magnetic fluctuations (Figure 1) lead to strong particle acceleration and the development of a clear nonthermal power-law spectrum within $1$-$2$ light crossing times. The spectral index $p\sim 2.8$ and does not appreciably change in the late stage of the simulation. Figure 5(b) shows the nonthermal spectra obtained at final times for simulations with $L/d_{e}\in\{512,1440,4096\}$ (normalized to the total number of particles in each simulation). By performing the aforementioned fitting procedure on these spectra, we find that the injection energy $\varepsilon_{\rm inj}$ is insensitive to the domain size $L$, whereas the cutoff energy $\varepsilon_{c}$ steadily increases with $L$. The power-law index $p$ steepens slightly with increasing domain size (see discussions below).

The spectral properties ($\varepsilon_{\rm inj}$, $\varepsilon_{c}$, $p$) are plotted against domain size $L$ for all of our simulations in Figure 6. We find that the simulation with $L/d_{e}=512$ was too small to yield precise measurements of these quantities (yielding a relatively large uncertainty), and therefore is not included. By inspecting Figure 5(b), we find the injection energy $\varepsilon_{\rm inj}$ to be insensitive to domain size, the power-law index $p$ to be slightly larger for larger domain sizes, and the cutoff energy $\varepsilon_{c}$ to be larger for larger domain sizes, in accordance with the trends in Figure 6.

Figure 6(a) shows that $p$ only weakly depends on $L$ and reaches $p\simeq 2.9$ for the largest $L/d_{e}=4096$, similar to Zhdankin et al. (2018). This weak dependence could be due to the decay of turbulence, leading to weaker acceleration in the late stage. The injection energy $\varepsilon_{\rm inj}$ shown in Figure 6(b) follows a similar trend, converging around $\varepsilon_{\rm inj}\simeq 10.5\,m_{e}c^{2}$ ($\simeq(\sigma_{0}/2)m_{e}c^{2}$) with an error $\pm\,0.5\,m_{e}c^{2}$. In contrast, $\varepsilon_{c}$ increases linearly with $L$ (Figure 6(c)), suggesting that particles can be accelerated to higher energies in simulations with larger domain sizes. Hence the power-law extent $R$ grows linearly with increasing domain size (Figure 6(d)), owing to the invariance of $\varepsilon_{\rm inj}$ and linear rise of $\varepsilon_{c}$ with increasing $L$.

To better understand particle acceleration mechanisms, we analyze the energy gains of individual tracer particles and break them down into the work done by parallel ($W_{\parallel}$) and perpendicular ($W_{\perp}$) electric fields. This is done by first using the tracked particle data to calculate the electric field parallel to the local magnetic field $\bm{E_{\parallel}}=(\bm{E}\cdot\bm{B}/{B}^{2})\bm{B}$ and perpendicular to it $\bm{E_{\perp}}=\bm{E}-\bm{E_{\parallel}}$. Then we can then calculate the work done by each component, i.e. $W_{\parallel}(t)\equiv q\int_{0}^{t}\textbf{v}(t^{\prime})\cdot\textbf{E}_{\parallel}(t^{\prime})\,dt^{\prime}$ and $W_{\perp}(t)\equiv q\int_{0}^{t}\textbf{v}(t^{\prime})\cdot\textbf{E}_{\perp}(t^{\prime})\,dt^{\prime}$.

Four examples of such tracer particles are shown in Figure 7, with horizontal dashed lines indicating the injection threshold $W_{\rm inj}\equiv\varepsilon_{\rm inj}-\varepsilon_{0}$ for each particle, which represents the energy gain necessary for the particle to cross the injection energy $\varepsilon_{\rm inj}$. Since the initial energy $\varepsilon_{0}$ of each particle is sub-relativistic (i.e., $\lesssim 1$), the injection thresholds $W_{\rm inj}$ hover just below the injection energy; in particular, $W_{\rm inj}\simeq$ (a) $9.5\,m_{e}c^{2}$, (b) $10.1\,m_{e}c^{2}$, (c) $10.3\,m_{e}c^{2}$ and (d) $10.1\,m_{e}c^{2}$ whereas $\varepsilon_{\rm inj}\simeq 10.5$ for the case $L/d_{e}=1440$.

In Figure 7(a), we see that for a high energy particle, the energy gain during injection is dominated by $W_{\parallel}$. Later, $W_{\parallel}$ flattens out, and $W_{\perp}$ dominates the energy gain. The pattern is similar to examples shown in Comisso & Sironi (2019) and has been seen in reconnection simulations (Guo et al., 2015; Kilian et al., 2020; French et al., 2023). Hence, the subsequent acceleration for this particle to high energies is a result of the perpendicular electric fields via a Fermi-like mechanism. Figure 7(b) shows a different high energy particle for which $W_{\parallel}$ flattens out at a much lower energy and $W_{\perp}$ dominates both the injection and post-injection phases. We also find relatively rare cases with $W_{\parallel}$ dominating the post-injection phase, shown in Figure 7(c) and (d).

Since every particle experiences a different evolution, our analysis is performed statistically over an ensemble of tracer particles (about 10-20% of all the tracers) whose final energy exceeds $\varepsilon_{\rm inj}$. Further, we monitor particles that cross certain energy thresholds $\varepsilon_{\rm threshold}$ separately. We break the energization process of each monitored particle into two phases: the energy gain up to the injection energy $\varepsilon_{\rm inj}$ termed pre-injection, and subsequent energy gain termed post-injection. The “pre-injection parallel share” is defined as the fraction of monitored particles which have $W_{\parallel}(t_{\rm inj})>W_{\perp}(t_{\rm inj})$ (where $t_{\rm inj}$ is the time step whereupon $\varepsilon=\varepsilon_{\rm inj}$ is reached). Similarly, the “post-injection parallel share” is defined as the fraction of monitored particles whose post-injection parallel energization exceeds perpendicular energization (i.e., $W_{\parallel}(t_{\rm final})-W_{\parallel}(t_{\rm inj})>W_{\perp}(t_{\rm final})-W_{\perp}(t_{\rm inj})$, where $t_{\rm final}$ is the final time step of the simulation). Figure 8 shows the parallel share for particles with final energy $\varepsilon_{\rm final}\geq\varepsilon_{\rm threshold}\in\{\varepsilon_{\rm inj},4\,\varepsilon_{\rm inj},16\,\varepsilon_{\rm inj}\}$.

We ran all of our simulations twice using the random number generator seeds to be 1 and 2. The values shown in Figure 8 are the average of these two simulations and the error bars end points are the actual values of the two simulations.

For $\varepsilon_{\rm threshold}=\varepsilon_{\rm inj}$ (blue line in Figure 8(a)), the pre-injection parallel share decreases with increasing domain size and drops to $\sim 50\%$ for the largest domain, implying that $W_{\parallel}$ and $W_{\perp}$ play a comparable role in the initial particle energization. However, this curve has not yet saturated with increasing domain size, suggesting that $W_{\perp}$ could dominate the injection stage for larger systems. As $\varepsilon_{\rm threshold}$ increases, the pre-injection parallel share also increases. For very high energy particles ($\varepsilon_{\rm threshold}=16\,\varepsilon_{\rm inj}$), the energy gain for most ($>90\%$) particles is dominated by $W_{\parallel}$ for small $L$. For larger $L$, the parallel share declines to $\simeq 75\%$. This decreasing trend again indicates that the pre-injection parallel share fraction for high-energy particles could be even smaller for larger systems.

The post-injection shares are converged with system size $L$ for each $\varepsilon_{\rm threshold}$. For $\varepsilon_{\rm threshold}=\varepsilon_{\rm inj}$, the parallel share is $\sim 50\%$, indicating that $W_{\parallel}$ and $W_{\perp}$ contribute comparably to particle energization in the post-injection phase. As $\varepsilon_{\rm threshold}$ increases, the post-injection parallel share decreases: When $\varepsilon_{\rm threshold}=4\,\varepsilon_{\rm inj}$, $W_{\parallel}$ contributes 20%, and for $\varepsilon_{\rm threshold}=16\,\varepsilon_{\rm inj}$, the $W_{\parallel}$ contribution is negligible. This indicates that for very high energy particles, $W_{\perp}$ dominates the post-injection energy gain for almost all particles.

In this paper, we have presented results from 2D PIC simulations with $\sigma_{0}=20$ and $L/d_{e}$ varying from $512$ to $4096$ to investigate the mechanisms of nonthermal particle acceleration in turbulent plasma.

We find that for $\varepsilon_{\rm threshold}=16\,\varepsilon_{\rm inj}$, the smaller domain sizes pre-injection parallel shares are higher than $90\%$, indicating that $W_{\parallel}$ dominates the pre-injection phase for most particles. This is in alignment with the results of Comisso & Sironi (2019), where they claim that initial particle acceleration is caused by $W_{\parallel}$. In the post-injection case for the same $\varepsilon_{\rm threshold}$, we find that the parallel share is close to $0\%$, which indicates that almost all high energy particles get most of their energy from $W_{\perp}$. This finding also aligns with Comisso & Sironi (2019), which shows $W_{\perp}$ dominates late-stage energization. However, it must be noted that the particles analyzed by Comisso & Sironi (2019) are all very high energy with $\varepsilon_{\rm threshold}=18\sigma_{0}$. Even for high energy particles, we find that the pre-injection parallel share starts to decrease and drops to $75\%$, indicating $W_{\parallel}$ only dominates the initial energization of three-quarters of the tracer particles. Given the decreasing trend continues at the largest box size (green line in Figure 8(a)), it is likely that the contribution by $W_{\parallel}$ in the pre-injection phase might be even smaller for astrophysical scale systems. Furthermore, when we look at the full picture by analyzing all injected tracer particles ($\varepsilon_{\rm threshold}=\varepsilon_{\rm inj}$), we recognize that $W_{\perp}$ plays a greater role in particle energization during the pre-injection phase, and $W_{\parallel}$ also a plays a more significant role in post-injection particle energization, especially particles with energy close to the lower bound of the power-law distribution.

We find strong agreement with Zhdankin et al. (2018) in how the power-law index $p$ depends on domain size $L$ (c.f., Figure 6). In particular, we find the power-law index to steadily steepen with increasing domain size, with $p\simeq 2.9$ when $L/d_{e}=4096$. However it is still unclear at which domain size $L/d_{e}$ and at what value $p$ will converge. Simulations with continuous driving may help resolve this issue.

Our simulations use a constant magnetization $\sigma_{0}=20$ and turbulence amplitude ${\delta B}_{\rm rms0}/B_{0}=1$ in an electron-positron plasma. If the mechanisms that underlie injection in relativistic turbulence are the same as those for relativistic magnetic reconnection (French et al., 2023; Vega et al., 2024), then the share of work done by $E_{\parallel}$ ($E_{\perp}$) could increase with magnetization (c.f., Fig. 29 of Zhdankin et al. (2020)), but decrease with the turbulence amplitude. While electrons and positrons undergo identical injection processes, protons may undergo significantly different processes and requires a future study. Recent studies show that proton injection and acceleration in turbulence and magnetic reconnection are dominated by perpendicular electric field (Comisso & Sironi, 2022; Zhang et al., 2024b). Further studies are needed to resolve these important issues.