A Model for Nonthermal Particle Acceleration in Relativistic Magnetic Reconnection

TL;DR

This paper introduces a semi-analytical model for nonthermal particle acceleration in relativistic magnetic reconnection, validated by kinetic simulations and applicable to astrophysical systems.

Contribution

It presents a systematic, semi-analytical framework for understanding particle energy distributions in magnetic reconnection, incorporating injection, acceleration, and escape processes.

Findings

Model accurately predicts spectral index and bounds of power-law distributions.

Simulation results agree with the model's predictions.

Application to astrophysical systems yields realistic power-law indices and break energies.

Abstract

The past decade has seen an outstanding development of nonthermal particle acceleration in magnetic reconnection in magnetically-dominated systems, with clear signatures of power-law energy distributions as a common outcome of first-principles kinetic simulations. Here we propose a semi-analytical model for systematically investigating nonthermal particle acceleration in reconnection. We show particle energy distributions are well determined by particle injection, acceleration, and escape processes. Using a series of kinetic simulations, we accurately evaluate the energy- and time-dependent model coefficients. The resulting spectral characteristics, including the spectral index and lower and upper bounds of the power-law distribution, agree well with the simulation results. Finally, we apply the model to predict the power-law indices and break energies in astrophysical reconnection systems.

Equations

$${\partial }_{t}f+{\partial }_{\epsilon }({\alpha }_{\text{acc}}\epsilon f)={\partial }_{\epsilon }^2(D_{\epsilon \epsilon }f)-{\alpha }_{\text{esc}}f+\frac{f_{\text{inj}}}{{\tau }_{\text{inj}}},$$

$${\alpha }_{\text{acc}\parallel }$$

$$\approx ec{E}_R{⟨{\mu }_E⟩}_{\text{X}}{\epsilon }^{-1}{F}_{\text{X}},$$

$$D_{0\parallel }\approx e^2{c}^2{E}_R^2{\epsilon }^{-2}{F}_{\text{X}}{\tau }_{\text{dec}\parallel }.$$

$${\alpha }_{\text{acc}\perp }=⟨-e{\mathbf{E}}_{\perp }\cdot \mathbf{v}/\epsilon ⟩,$$

$${\alpha }_{\text{acc}\perp }\sim c\frac{{\Gamma }_{Ax}^2\left(1+2{\beta }_{Ax}(1+b_g^2{)}^{-1/2}+{\beta }_{Ax}^2\right)-1}{L_{\text{O}}(1+b_g^2{)}^{1/2}},$$

$$D_{0\perp }\sim \frac{1}{3}{\Gamma }_{Ax}^2{\beta }_{Ax}^2{L}_{\text{O}}^{-1}(1+b_g^2{)}^{-1/2}c.$$

$${\alpha }_{\text{esc}}$$

$$\sim V_{Ax}/L_{\text{O}}.$$

$${\alpha }_{\text{acc}\parallel }\sim 2.5\times 10^{-5}\frac{{\bar{\epsilon }}_t^3(\bar{\epsilon }/{\bar{\epsilon }}_t{)}^2}{{\left(1+(\bar{\epsilon }/{\bar{\epsilon }}_t{)}^2\right)}^{3/2}},$$

$${\alpha }_{\text{esc}}\sim 6\times 10^{-4}\frac{{\bar{\epsilon }}^2{\bar{\epsilon }}_t^2}{{\bar{\epsilon }}^2+{\bar{\epsilon }}_t^2}\left[\left(1-\frac{{\bar{\epsilon }}_{t0}^2}{{\bar{\epsilon }}_t^2}\right)\frac{10{\bar{\epsilon }}_t^2}{10{\bar{\epsilon }}_t^2+{\bar{\epsilon }}^2}+\frac{{\bar{\epsilon }}_{t0}^2}{{\bar{\epsilon }}_t^2}\right],$$

$$p\approx \frac{1}{{\sigma }_x+0.2(1+\tanh (b_g))}+0.04\tanh (b_g){\sigma }_x+1.7b_g+2.1,$$

$$\frac{{\rho }_b}{L}=0.007(1-e^{-0.6{\sigma }_x})e^{-0.06\tanh (b_g){\sigma }_x}(1+b_g^2{)}^{-1}$$

$$d\epsilon =({\alpha }_{\text{acc}\parallel }+{\alpha }_{\text{acc}\perp })\epsilon dt+\sqrt{2(D_{0\parallel }+D_{0\perp }){\epsilon }^2}d{W}_t,$$

$${\epsilon }_{\text{O}}=e{B}_0\sqrt{1+b_g^2}{L}_{\text{O}}/(10m_{e}c{\Gamma }_{Ax}),$$

Full text

High energy astrophysics (739); Plasma astrophysics (1261)

††journal: ApJL

1 Introduction

Nonthermal particle acceleration processes during magnetic reconnection are plausible mechanisms responsible for high-energy emissions observed in magnetically dominated systems, such as relativistic jets from gamma-ray bursts (GRBs) (Zhang & Yan, 2011; McKinney & Uzdensky, 2012) and active galactic nuclei (AGN) (Giannios et al., 2009; Zhang et al., 2015, 2018), pulsar wind nebulae (PWNe) (Uzdensky & Spitkovsky, 2014), and solar flares (Lin et al., 2003). Recently, interest has surged in studying particle acceleration and resulting particle energy spectra, mainly via particle-in-cell simulations (Guo et al., 2014; Sironi & Spitkovsky, 2014). These studies have successfully obtained power-law energy spectra $f(\varepsilon)\propto\varepsilon^{-p}$ with spectral index $p$ decreasing with the plasma magnetization $\sigma$ and approaching 1 when $\sigma\gtrsim 10$ (Sironi & Spitkovsky, 2014; Guo et al., 2014, 2016; Werner et al., 2016). However, the origin of the nonthermal energy spectrum and how it varies in different situations is still undergoing active debate.

Despite the advances in kinetic simulations (Hoshino et al., 2001; Zenitani & Hoshino, 2001; Sironi & Spitkovsky, 2014; Liu et al., 2011; Dahlin et al., 2014; Guo et al., 2014, 2016; Li et al., 2015, 2018), theoretical explorations on particle acceleration have limited success in predicting the resulting energetic particle spectrum. Although a simple derivation shows $p\sim 1$ when particle escape is ignored (Guo et al., 2014, 2015), numerous simulations have consistently obtained $p>2$, indicating some key physics is still missing (Werner et al., 2018; Ball et al., 2018; Uzdensky, 2022). In addition, how the resulting energy spectra vary with the guide field, and how the results from kinetic simulations can be extrapolated to astrophysical scales are not clear.

In this Letter, we present a model for particle acceleration in magnetic reconnection. In contrast to previous studies, we consider particle injection, acceleration, and escape together as key ingredients for understanding the energy spectra and evaluate them in a reconnection system. In a series of fully kinetic simulations, we quantify the acceleration (both first order and second order) and escape processes with their time and energy dependence. We use the model coefficients to predict the nonthermal power-law characteristics (spectral indices, lower and higher bounds of the power-laws) based on a Fokker-Planck approach and find that the prediction matches the simulations very well. While our model is capable of explaining nonthermal acceleration in kinetic simulations in general, it also reveals essential physical factors that determine energetic particle spectrum, making a crucial step toward a comprehensive understanding of those processes in large-scale astrophysical systems.

2 A Model

Fig. 1(a) illustrates our model. The reconnection inflow continuously brings thermal particles into the reconnection layer, where a fraction of them can be accelerated out of the thermal pool through particle injection processes (Ball et al., 2019; Kilian et al., 2020; Sironi, 2022; Guo et al., 2023; French et al., 2023). A primary acceleration phase then leads to the formation of power-law distribution. While most reconnection studies focus on the first-order acceleration (Guo et al., 2014, 2015), our model includes the second-order Fermi acceleration, which is found to be comparable to the first-order mechanism (see below). This second-order acceleration can be due to inhomogeneous distribution of island acceleration, or turbulent like behavior, similar to relativistic turbulence (Wong et al., 2020; Comisso & Sironi, 2018, 2019). In addition, we find that particle escape plays an important role in determining the spectral shape. Since particles do not participate in the acceleration once they are trapped by the largest magnetic islands/flux ropes, they are considered separately as an effective “escape” process. Together, we consider these processes in a Fokker-Planck approach (Blandford & Eichler, 1987) describing the evolution of the particle energy distribution $f(\varepsilon)$, according to

[TABLE]

where $\varepsilon=(\gamma-1)m_{s}c^{2}$ is the kinetic energy, $\alpha_{\text{acc}}$ is the acceleration rate, $D_{\varepsilon\varepsilon}\coloneqq D_{0}\varepsilon^{2}$ is the energy diffusion coefficient, $\alpha_{\text{esc}}\equiv\tau_{\text{esc}}^{-1}$ is the escape rate, $f_{\text{inj}}$ is the injected thermal particle distribution, and $\tau_{\text{inj}}$ is particle injection time scale. $\alpha_{\text{acc}}\coloneqq(\partial_{t}\varepsilon+\partial_{\varepsilon}D_{\varepsilon\varepsilon})\varepsilon^{-1}$ describes a combination of the first-order Fermi processes and the accompanying first-order term associated with second-order Fermi mechanisms. Instead of separating the injection and nonthermal acceleration, we will treat them as a continuous process and determine $\alpha_{\text{acc}}$ and $D_{0}$ due to parallel or perpendicular electric fields ($\bm{E}_{\parallel}$ and $\bm{E}_{\perp}$).

We assume that reconnection starts from a single elongated current sheet (with a length of $L_{x}$), which breaks into a series of magnetic islands (or flux ropes in 3D) (Loureiro et al., 2007; Bhattacharjee et al., 2009). These islands tend to merge to form larger islands, and secondary islands are continuously generated. The entire process lasts a few Alfvén-crossing times $\tau_{A}\coloneqq L_{x}/V_{Ax}$ until most magnetic flux is reconnected, where $V_{Ax}$ is the reconnection outflow speed. The total number of particles in the reconnection layer due to the reconnection inflow $dN/dt\propto RV_{Ax}$ (where the reconnection rate $R\sim 0.1$ (Liu et al., 2017; Goodbred & Liu, 2022)). The largest islands grow in the size of $L_{\text{O}}=d_{i}+RV_{Ax}t$, assuming they are initially in the $d_{i}$ (ion inertial length) scale. For the rest of the discussion, $n_{\text{acc}}(\varepsilon)$ and $n_{\text{esc}}(\varepsilon)$ are the numbers of the accelerating and escaped particles in different energy bins, respectively.

For a single electron, the acceleration rate due to the parallel electric field is $-e\bm{E}_{\parallel}\cdot\bm{v}/\varepsilon$. Since $E_{\parallel}\approx E_{R}=RV_{Ax}B_{0}/c$ near the X-lines and $\approx 0$ rest of the reconnection layer, the acceleration rate due to $\bm{E}_{\parallel}$ is

[TABLE]

where $\left<\cdots\right>$ indicates the ensemble average over $n_{\text{acc}}(\varepsilon)$ particles, $\left<\cdots\right>_{\text{X}}$ is the average over $n_{\text{X}}(\varepsilon)$ particles near the X-lines, and $\mu_{E}\coloneqq-\bm{v}\cdot\bm{E}_{\parallel}/(vE_{\parallel})$. The filling factor of these electrons $F_{\text{X}}\sim d_{e}/L_{\text{O}}\sim t^{-1}$ for large $t$. $\left<\mu_{E}\right>_{\text{X}}\approx 0$ for low-energy thermal electrons with nearly isotropic distributions and increases when they are accelerated by $\bm{E}_{\parallel}$ or the Fermi mechanism (Drake et al., 2006). We will determine the exact energy dependence using the simulation results. The energy diffusion rate due to $\bm{E}_{\parallel}$ is $D_{0\parallel}\coloneqq\left<(\delta\alpha_{\text{acc}\parallel})^{2}\right>\tau_{\text{dec}\parallel}\approx\left(e^{2}c^{2}E_{R}^{2}\left<\mu_{E}^{2}\right>_{\text{X}}\varepsilon^{-2}F_{\text{X}}-\alpha_{\text{acc}\parallel}^{2}\right)\tau_{\text{dec}\parallel}$, where $\delta\alpha_{\text{acc}\parallel}=-e\bm{E}_{\parallel}\cdot\bm{v}/\varepsilon-\alpha_{\text{acc}\parallel}$ is the fluctuation of $\alpha_{\text{acc}\parallel}$ and $\tau_{\text{dec}\parallel}\sim d_{i}/c$ is the decorrelation time for $\bm{E}_{\parallel}$ (the timescale on which particles see decorrelated $\bm{E}_{\parallel}$, le Roux et al., 2015). It indicates how fast the broadening of the energy distribution is due to $\bm{E}_{\parallel}$. Since $\alpha_{\text{acc}\parallel}^{2}\propto F_{\text{X}}^{2}\propto t^{-2}$ for large $t$, the first term dominates. For particles with an anisotropy along the magnetic field, $\left<\mu_{E}^{2}\right>_{\text{X}}$ changes slowly between 0.5 and 1. Therefore,

[TABLE]

The acceleration rate due to $\bm{E}_{\perp}$ is

[TABLE]

as particles gaining energy through the Fermi mechanism when colliding with the reconnection outflow (Guo et al., 2014). The energy gain of each collision is about $\Delta\varepsilon=\left(\Gamma_{Ax}^{2}\left(1+2V_{Ax}v_{x}/c^{2}+V_{Ax}^{2}/c^{2}\right)-1\right)\varepsilon$, and the collision time scale is about $L_{\text{O}}/v_{x}$, where $\Gamma_{Ax}=(1-V_{Ax}^{2}/c^{2})^{-1/2}$, and $v_{x}$ is the particle velocity along the reconnection outflow direction. When there is a guide field $B_{g}=b_{g}B_{0}$, $v_{x}\sim c/(1+b_{g})^{1/2}$ on average, and $V_{Ax}=c\sqrt{\sigma_{x}/(\sigma+1)}$ (Liu et al., 2015). $\sigma_{x}=B_{0}^{2}/4\pi w$ is plasma magnetization based on the reconnection magnetic field $B_{0}$ (Liu et al., 2015), where $w$ is the enthalpy density. $\sigma=B^{2}/4\pi w$ is total plasma magnetization, where $B=(B_{0}^{2}+B_{g}^{2})^{1/2}$. Thus, $V_{Ax}$ decreases with $b_{g}$ and increases with $\sigma_{x}$. We obtain

[TABLE]

where $\beta_{Ax}=V_{Ax}/c$. When the power-law index $p<2$ in highly magnetized plasmas, most of the kinetic energy is in the high-energy tail (Sironi & Spitkovsky, 2014; Guo et al., 2014), and the enthalpy density $w$ keeps increasing. The resulting momentum flux density $w\Gamma_{Ax}^{2}V_{Ax}^{2}$ will keep increasing if $V_{Ax}$ stays constant, which is not sustainable. Thus, we expect the reconnection outflow to decrease gradually in highly magnetized plasmas [See Supplemental Material A], resulting in an approximately constant momentum flux density. The energy diffusion rate due to the Fermi mechanism is $D_{0\perp}\sim\Gamma_{Ax}^{2}\beta_{Ax}^{2}c\lambda_{\text{mfp}}^{-1}/3$ (Lemoine, 2019; Comisso & Sironi, 2019), where $\lambda_{\text{mfp}}\sim L_{\text{O}}(1+b_{g}^{2})^{1/2}$ is scattering mean-free-path. Thus,

[TABLE]

Since energetic particles escape the acceleration region when advected by the reconnection outflow, the escape rate is

[TABLE]

3 Numerical Simulations

We test the proposed model through a series of 2D particle-in-cell kinetic simulations of magnetic reconnection in a proton-electron plasma in the $xz$–plane. The simulations all start from a forcefree current sheet (Guo et al., 2014). The boundaries along $x$ are periodic for both fields and particles. The boundaries along $z$ are reflective for particles and perfectly conducting for fields. A long-wavelength perturbation with $\delta B_{z}=0.03B_{0}$ is included to initiate reconnection (Birn et al., 2001). We performed the simulations in proton-electron plasmas ($m_{i}/m_{e}=1836$) similar to earlier studies (Werner et al., 2018; Kilian et al., 2020; Ball et al., 2019). The cold ion magnetization parameters $\sigma_{ic}\coloneqq B_{0}^{2}/(4\pi n_{i}m_{i}c^{2})\in[0.4,1.6,6.4,25.6,102.4]$. The guide fields $b_{g}\in[0.0,0.5,1.0,2.0]$. The plasmas are initially uniform with density $n_{0}$ and follow the Maxwell–Jüttner distributions with a normalized electron temperature $\theta_{e}=kT_{e}/m_{e}c^{2}=10\sigma_{ic}$ and $T_{i}=T_{e}$, resulting in $\sigma_{x}=0.39,1.51,5.2,13.0,19.7$ for the runs with different $\sigma_{ic}$. The domain sizes are $L_{x}\times L_{z}=L_{0}\times L_{0}/2$, where $L_{0}=4096\sqrt{\theta_{e}}d_{e0}$ and $d_{e0}\coloneqq c\sqrt{m_{e}/4\pi n_{0}e^{2}}$ is the nonrelativistic electron inertial length. The domain is resolved with a grid $n_{x}\times n_{z}=8192\times 4096$, and the resulting cell sizes $\Delta x=\Delta z=0.5\sqrt{\theta_{e}}d_{e0}=0.5d_{e}$, where $d_{e}$ is the electron inertial length including the relativistic correction. With $b_{g}=0.0$, we performed three additional simulations with $L_{x}\in[L_{0}/4,L_{0}/2,2L_{0}]$ and the same cell sizes for each $\sigma_{ic}$ to examine the system size dependence. We use 100 particles/cell/species and track about 1 million tracer particles at a high cadence to calculate the model coefficients. These tracer particles are a small, uniformly selected subset of all electrons, and they evolve self-consistently along with the rest of the particles.

4 Model-Simulation Comparison

In all simulations, electron energy spectra develop power-law tails with various power-law indices. Fig. 2(a) shows a sustainable power-law spectrum with a stable $p=2.5$ and cutoff energy $\bar{\varepsilon}_{b}$ linearly increasing with time. This result suggests that the cutoff energy increases with the box sizes, confirmed in Fig. 2(b). As the magnetization increases, the reconnection outflow and the acceleration rate due to the motional electric field will become faster, resulting in a harder spectrum (see Fig. 5). Similarly, the guide field will slow down the outflow and the Fermi acceleration, leading to a softer spectrum (see Fig. 5). We will demonstrate that the model can capture these dependences.

To evaluate the model coefficients in simulations, we must separate particles into the accelerating and the escaped populations. Fig. 1(b) & (c) shows two electron trajectories to illustrate the escape mechanism. The two electrons gain substantial energies before being confined inside the large magnetic island. Their kinetic energies fluctuate without much change after that. Thus, they are treated as escaped particles after their acceleration stops. After separating the two populations, one can calculate acceleration and escape rates using Eqs. (2), (5), and  (8). The process is repeated at every step of the tracer particles, allowing us to measure how the rates change over time. We evaluate $D_{0}$ by tracking the energy spread of the particles in different energy bins (Comisso & Sironi, 2019; Wong et al., 2020). Fig. 3 shows the energy dependence of the rates in the runs with $\sigma_{ic}=102.4$. $\alpha_{\text{acc}\parallel}$ peaks around $\bar{\varepsilon}=10$ and decreases when $\bar{\varepsilon}>10$, following the scaling in Eq. (3). Fig. 3(b) shows that $\alpha_{\text{acc}\parallel}$ peaks at higher energies when $b_{g}>0$, due to the increasing $\left<\mu_{E}\right>_{\text{X}}$ with $b_{g}$ (Li et al., 2018). The transition energy $\bar{\varepsilon}_{t}=\varepsilon_{t}/(kT_{e})\approx 12+10\tanh(2b_{g})$ according to the simulations, and the dashed lines in Fig. 3(b) follow the energy dependence

[TABLE]

which is $\sim\bar{\varepsilon}^{2}$ for low-energy particles and $\sim\bar{\varepsilon}^{-1}$ for high-energy particles (Eq. 3). This energy dependence is determined by fitting the energy-dependent $\alpha_{\text{acc}\parallel}$ to achieve a smooth transition between the two parts. As the rates exhibit significant temporal fluctuations, as illustrated in Fig. 4, the coefficient should be treated as an approximation rather than an exact value. $\alpha_{\text{acc}\perp}$ is nearly constant, as expected from Eq. (6). In Fig. 3(a), we adopt a $\beta_{Ax}=0.8<\sqrt{\sigma_{x}/(\sigma+1)}$ due to the decelerating reconnection outflow. Fig. 3(a) shows that while the primary acceleration mechanism at high energies is the Fermi process, the direct and Fermi mechanisms are comparable at low energies, consistent with prior findings by Guo et al. (2019). Regarding the escape, $\alpha_{\text{esc}}$ increases with particle energy until it saturates at a constant $\sim V_{Ax}/L_{\text{O}}$, consistent with Eq. (9). When $b_{g}$ is finite, the larger anisotropy indicates that more particles can escape from the acceleration region by streaming along the magnetic field lines. Consequently, $\alpha_{\text{esc}}$ peaks around $\varepsilon_{t}$ (Fig. 3(c)), at which the anisotropy is the largest. We use the following function to model the energy dependency of $\alpha_{\text{esc}}$ from low to high energies.

[TABLE]

where $\bar{\varepsilon}_{t0}$ is $\bar{\varepsilon}_{t}(b_{g}=0)$. $\alpha_{\text{esc}}\propto\varepsilon^{2}$ for low-energy particles and approaches $\sim V_{Ax}/L_{\text{O}}$ for high-energy particles (Eq. 9). Similar to Eq. 10, the energy dependence is established by fitting energy-dependent $\alpha_{\text{esc}}$ to achieve a smooth transition between the low-energy and high-energy parts. The complex form accounts for the peak that arises near the transition energy $\varepsilon_{t}$. Fig. 3 (c) shows some discrepancies at low energies, especially for $b_{g}=2$ (purple curve) case. Such discrepancies are caused by the strong temporal fluctuations in the rates (Fig. 4), likely due to the initial condition or boundary conditions in the simulations. Fig. 3(a) shows that $D_{0}$ is dominated by $D_{0\parallel}\propto\varepsilon^{-2}$ (Eq. 4) at low energies and $D_{0\perp}\sim$ constant (Eq. 7) at high energies.

Fig. 4 shows how the rates change with time. Fig. 4(a) shows that $\alpha_{\text{acc}\parallel}\propto t^{-1}$, following the scaling in Eq. (3). It is larger when there is a finite guide field, agreeing with Fig. 3(b). Although $\alpha_{\text{acc}\perp}$ strongly fluctuates (Fig. 4(b)), it decreases with time, as expected from Eq. (6). $\alpha_{\text{acc}\perp}$ can be $\propto t^{-2}$ when $\sigma$ is large enough due to the decelerating reconnection outflow [See Supplemental Material A]. Fig. 4(c) shows that $\alpha_{\text{esc}}\sim t^{-1}$, consistent with Eq. (9). For more model-simulation comparison of the time dependence, see Supplement Material B.

We then solve the Fokker-Planck equation using the energy- and time-dependent rates [See Supplemental Material C for details].

Fig.5 (a) and (b) show that the modeled spectra (thick lines) agree well with those in the simulations (thin lines). The model can capture the spectra’ $\sigma$- and $b_{g}$- dependences and the power-law break energies $\varepsilon_{b}$. Despite some notable agreements between the simulation and model results, as depicted by the overlaying orange curves in Fig. 5 (a) and the green and blue curves in Fig. 5 (b), there exist discrepancies between the simulation results and model results. Such discrepancies can be attributed to the settings of the simulations, such as the initial and boundary conditions, as well as the constrained system sizes. Consequently, the model results are anticipated to reveal the overall patterns but not necessarily replicate the simulation results exactly.

Since the model is simple enough, we can predict the spectra for much larger systems and study their long-term evolution. The insets of panels (a) and (b) show that $p$ gradually increases until saturation at a constant value in large enough systems. Since the guide field plays a more important role in the high-$b_{g}$ reconnection, it takes a larger box (up to $10^{3}L_{0}$ when $b_{g}=2$) for the spectra to saturate. The model also shows that the break energy linearly increases with the system size up to $10^{6}$ times the PIC simulation sizes. Thus, the model can predict the spectra in astrophysical reconnection sites. Fig. 5 (c) and (d) show the spectral indices and the gyroradius of the break energies, normalized by the system sizes $L$. When $b_{g}=0$, $p$ approaches 2.1 as the $\sigma_{x}$ increases, indicating that the power-law extension can keep growing without causing the kinetic energy to diverge. When $b_{g}$ is finite, the spectra are the hardest when the $\sigma_{x}$ is between 1 and 10. $p$ can be approximated as

[TABLE]

where the first two terms capture the variation with $\sigma_{x}$ when the guide field is finite. When $b_{g}=0$, $p\approx(\sigma_{x}+0.2)^{-1}+2.1$. For a fixed magnetization $\sigma_{x}$, the spectrum becomes softer as the guide field becomes stronger. Fig. 5 (d) shows the gyroradius of the electrons with the break energies ($\rho_{b}\coloneqq(\varepsilon_{b}/m_{e}c^{2}+1)m_{e}c/eB_{0}$) is a fraction of the system size $L$. When $b_{g}=0$, $\rho_{b}$ increases with $\sigma_{x}$ and approaches 0.007 times the system sizes at $2\tau_{A}$. When $b_{g}$ is finite, $\rho_{b}$ peaks when $\sigma_{x}$ is a few, suggesting the acceleration is most efficient when the $\sigma_{x}$ is between 1 and 10. $\rho_{b}$ can be approximated as

[TABLE]

at $2\tau_{A}$. When $b_{g}=0$, $\rho_{b}/L=0.007(1-e^{-0.6\sigma_{x}})$. When $b_{g}$ is finite, $\rho_{b}$ peaks when the $\sigma_{x}$ is between 1 and 10. For a fixed $\sigma_{x}$, $\rho_{b}$ monotonically decreases with the guide field. Note that Eqs. 12 and 13 rely only on the hot magnetization parameter $\sigma_{x}$ and the guide field strength $b_{g}$. Nonetheless, the equations are applicable to proton-electron plasmas under the conditions where $0.1<\sigma_{x}<100$ and electrons are relativistic. Although the equations can be potentially used for pair plasmas as they do not depend on the ion-to-electron mass ratio, they have not been thoroughly verified against fully kinetic simulations using pair plasmas.

5 Discussion and Conclusions

The present study builds upon the prior research by Guo et al. (2014, 2015) and makes improvements in the following aspects: firstly, instead of separating the injection and the nonthermal acceleration, it treats them as a continuous process and incorporates the energy-dependence of acceleration rates; secondly, it incorporates 2nd-order acceleration, including those due to the parallel electric field and the Fermi mechanism; thirdly, it accurately evaluates particle escape rate in the simulations and compares it with the model results; fourthly, it extends the model to proton-electron plasmas; and lastly, besides plasma magnetization, it considers the effect of guide fields on particle acceleration. While these processes are more or less studied in previous studies individually, this study is the first attempt to put these processes together in the same model.

Some important physics is still missing in the current model and could be the key to addressing the particle acceleration in the astrophysical reconnection sites. First, previous research has demonstrated that radiative cooling, such as synchrotron or inverse Compton cooling, significantly impacts the acceleration mechanisms and the dynamics of the current layer, even when cooling is only marginally important (Zhang et al., 2018; Werner et al., 2019; Hakobyan et al., 2019; Sironi & Beloborodov, 2020; Sridhar et al., 2021). Second, 3D reconnection can self-generate plasma turbulence (Daughton et al., 2011; Guo et al., 2014), which could affect particle transport and acceleration (Dahlin et al., 2015; Li et al., 2019; Zhang et al., 2021b), and change particle acceleration and escape mechanisms (Zhang et al., 2021a). However, other studies have shown that 3D physics is not as crucial in the relativistic regime as in the non-relativistic regime (Guo et al., 2021). Our preliminary results indicate that the spectra are not saturated to the end of the simulation due to the limited system sizes of the 3D simulations. Thus, incorporating 3D effects in the model requires more studies using massive 3D kinetic simulations, which are currently not affordable for revealing the scalings in system sizes. Third, the simulation setup, including boundary conditions and initialization (Ball et al., 2018), may impact the energy dependence and the time evolution of the rates, as shown in Figs. 3 and 4. These could lead to the deviation of the modeling predictions from the simulation results, given that the semi-analytical model still needs empirical fitting of the kinetic simulation results. These simulations are essential in determining the energy dependencies at low energies and quantifying the constants in each rate for different simulation setups. Addressing all these factors using fully kinetic simulations is critical to making the model more reliable for extending the kinetic simulation results to astrophysical reconnection systems.

In this Letter, we presented a model for determining the main characteristics of power-law particle energy spectra in magnetic reconnection in magnetically dominated plasmas. The power-law spectrum produced by this model is controlled by particle injection, acceleration, energy diffusion, and escape processes. Using a series of first-principles 2D kinetic simulations, we evaluate these transport coefficients and obtained nice agreements with simulation results. By solving the Fokker-Planck equation describing these processes using the modeled transport coefficients, we demonstrate that the modeled spectra agree well with the simulation results for a broad range of magnetization and guide fields. We then use the simulation-verified model to predict the power-law spectral indices and break energies for astrophysically relevant reconnection systems and provide their empirical expressions for different magnetization and guide fields. Our results have strong implications to understanding nonthermal particle acceleration and emissions in high-energy astrophysical systems.

We thank the anonymous referee for a constructive review. We gratefully acknowledge our discussions with William Daughton and Patrick Kilian. X.L. acknowledges the support from NASA through Grant 80NSSC21K1313, National Science Foundation Grant No. AST-2107745, and Los Alamos National Laboratory through subcontract No. 622828. The simulations used resources provided by the National Energy Research Scientific Computing Center (NERSC, a U.S. Department of Energy Office of Science User Facility at Lawrence Berkeley National Laboratory) and the Texas Advanced Computing Center (TACC, at the University of Texas at Austin).

Appendix A Decelerating Reconnection Outflow

Fig. 6 shows that the reconnection outflow is decelerating as reconnection evolves. The deceleration is the strongest when the guide field is weak and the particle acceleration is the most efficient.

Appendix B Time Dependence of the Rates

Fig. 7 shows how the rates change with time for runs with different $b_{g}$. Although $\alpha_{\text{acc}\perp}$ strongly fluctuates (Fig. 7(a)), it decreases with time and $b_{g}$, as expected from Eq. (6) in the main text. When $b_{g}=0$, due to the decelerating reconnection outflow [See Supplemental Material C], $\alpha_{\text{acc}\perp}$ decreases faster than $t^{-1}$. Fig. 7(b) shows that $\alpha_{\text{esc}}$ is higher when there is a finite guide field. $\alpha_{\text{esc}}\sim t^{-1}$ (consistent with Eq. (9) in the main text) even as the reconnection outflow slows down when $b_{g}=0$ and $\sigma$ is large (Fig. 7(f)) because $V_{Ax}$ can be close to $c$ while $\Gamma_{Ax}$ decreases significantly. The energy diffusion rate $D_{0}$ strongly fluctuates but gradually decreases with $t$. When $b_{g}=0$, it decreases faster than $t^{-1}$ due to the decelerating outflow.

Appendix C Solving the Fokker-Planck Equation

We solve the Fokker-Planck equation (1) using the energy- and time-dependent rates. The acceleration part of Eq. (1) is equivalent to the stochastic differential equation (SDE) of the Itô type,

[TABLE]

where $W_{t}$ is the standard Wiener process, and $dW_{t}$ is the normalized distributed random number with mean 0 and variance $\Delta t$. The SDE can be solved using pseudo particles. We apply a taper function $\varepsilon_{\text{O}}/(\varepsilon_{\text{O}}+\varepsilon)$ to $\alpha_{\text{acc}\perp}$ and $D_{0\perp}$ when solving the equation to model the reduction of the rates with particle energies when their gyroradius is close to the largest island size. Specifically, we use

[TABLE]

which includes the correction when the outflow is relativistic. Motivated by the simulation results, we model the reconnection rate $R\approx 0.1V_{Ax}/V_{A0}$, where $V_{A0}=\sqrt{\sigma_{x}/(\sigma_{x}+1)}$. As a result, the number of injected particles and the parallel electric field have a strong dependence on the guide field. After updating the energies of the pseudo particles at each time step according to Eq. (C1), we determine whether the particles escape according to the escape rate and inject new pseudo particles. Before the next cycle, we evaluate the enthalpy density of the particles and update the outflow speed $V_{Ax}$ assuming a constant momentum flux density.