---
abstract: 'A model of global magnetic reconnection rate in relativistic collisionless plasmas is developed and validated by the fully kinetic simulation. Through considering the force balance at the upstream and downstream of the diffusion region, we show that the global rate is bounded by a value $\sim 0.3$ even when the local rate goes up to $\sim O(1)$ and the local inflow speed approaches the speed of light in strongly magnetized plasmas. The derived model is general and can be applied to magnetic reconnection under widely different circumstances.'
author:
- 'Yi-Hsin Liu'
- Michael Hesse
- Fan Guo
- William Daughton
- Hui Li
title: A model of global magnetic reconnection rate in relativistic collisionless plasmas 
---

[*Introduction–*]{} Magnetic fields often serve as the major energy reservoirs in high energy astrophysical systems, such as pulsar wind nebulae [@coroniti90a; @arons12a; @lyubarsky01a; @devore15a], gamma-ray bursters [@thompson94a; @zhangB11a; @mckinney12a] and jets from active galactic nuclei [@beckwith08a; @giannios10a; @jaroschek04a], where relativistic cosmic rays and gamma rays of energies up to TeV are generated explosively [@abdo11a; @bottcher13a]. Among the proposed physics processes (e.g.,[@sironi15a; @uhm14a; @zweibel09a]) that could unleash the magnetic energy, magnetic reconnection is considered to be a promising mechanism. For comparison, collisionless shocks, regarded to be efficient for particle acceleration in weakly magnetized plasmas, are inefficient in dissipating energy and accelerating non-thermal particles in magnetically dominated flows [@sironi15a]. Hence the study of magnetic reconnection in these exotic systems continues to be an interesting topic in high energy astrophysics.

One of the most important issues in relativistic reconnection studies is how fast magnetic energy can be dissipated in the reconnection layer, which determines the time scale of the explosive energy release events. Another related problem is the mechanism of non-thermal particle acceleration [@YYuan16a; @werner16a; @FGuo16a; @FGuo15a; @FGuo14a; @melzani14b; @sironi14a; @cerutti12a; @bessho12a; @zenitani01a]. Proposed mechanisms include the direct acceleration by the reconnection electric field at the diffusion region [@zenitani01a; @uzdensky11a], the Fermi mechanism at the outflow regions that involves particles bouncing back and forth between reconnection outflows emanated from different x-lines [@FGuo14a; @dahlin14a; @drake06a], and many other ideas (e.g., [@zank14a; @drury12a; @pino05a]). In collisionless plasmas, the energy gain of a particle must come from the work done by the electric field $\sim q\int{{\bf E}\cdot {\bf v} dt}$. Thus, determining the reconnection electric field in the relativistic limit is crucial to determine the acceleration rate and efficiency.

In such magnetically-dominated plasmas, the magnetic energy density is much larger than the rest mass energy density and the Alfvén speed approaches the speed of light. Early theoretical work suggested that the magnetic reconnection rate in the relativistic limit may increase compared to the non-relativistic case due to the enhanced inflow arising from the Lorentz contraction of plasma passing through the diffusion region [@blackman94a; @lyutikov03a]. However, it was later pointed out that the thermal pressure within a pressure-balanced current sheet will constrain the outflow to mildly relativistic conditions, where the Lorentz contraction is negligible [@lyubarsky05a] and a relativistic inflow is therefore impossible.

Recently, fully kinetic simulations by Liu et al. [@yhliu15a] showed that the local inflow speed approaches the speed of light, and the reconnection rate normalized to the immediately upstream condition of the diffusion region can be enhanced to $\sim O(1)$ in strongly magnetized plasmas. However, the global reconnection rate normalized to the far upstream asymptotic value remains $\lesssim 0.3$ [@FGuo15a; @yhliu15a; @melzani14a; @sironi14a; @bessho12a; @sironi16a] and this discrepancy is not understood. While the relativistic resistive-Petschek model [@petschek64a] suggests a similar value for the global rate [@lyubarsky05a], to realize a Petschek solution requires an [*ad hoc*]{} localized resistivity [@biskamp86a; @sato79a], otherwise, the current sheet collapses to the long Sweet-Parker layer [@sweet58a; @parker57a]. A mechanism for the localized diffusion region is therefore essential to model the reconnection rate. In this Letter, we derive the relation between the global rate and the degree of localization through considering the force balance at the upstream and downstream of the diffusion region. We then propose a mechanism that naturally leads to the localization in such collisionless plasmas.

[*Simulation setup–*]{} The kinetic simulation is performed using a Particle-in-Cell code- VPIC [@bowers09a], which solves the fully relativistic dynamics of particles and electromagnetic fields. The relativistic Harris sheet [@yhliu15a; @kirk03a; @zenitani07a; @wliu11a; @bessho12a; @melzani14a] is employed as the initial condition. The initial magnetic field ${\bf B}=B_{x0} \mbox{tanh}(z/\lambda) \hat{\bf x}$ corresponds to a layer of half-thickness $\lambda$. Each species has a distribution $f_h \propto \mbox{sech}^2(z/\lambda)\mbox{exp}[-\gamma_d(\gamma_Lmc^2+ mV_d u_y)/T']$ in the simulation frame, which is a component with a peak density $n'_0$ and temperature $T'$ boosted by a drift velocity $\pm V_d$ in the y-direction for positrons and electrons, respectively. In addition, a non-drifting background component $f_b \propto \mbox{exp}(-\gamma_L m c^2/T_b)$ with a uniform density $n_b$ is included. Here ${\bf u}=\gamma_L {\bf v}$ is the the space-like components of 4-velocity, $\gamma_L=1/[1-(v/c)^2]^{1/2}$ is the Lorentz factor of a particle, and $\gamma_d \equiv 1/[1-(V_d/c)^2]^{1/2}$. The drift velocity is determined by Ampére’s law $cB_{x0}/(4\pi\lambda)=2 e\gamma_d n'_0 V_d $. The temperature is determined by the pressure balance $B_{x0}^2/(8\pi)=2 n'_0 T'$. The resulting density in the simulation frame is $n_0=\gamma_d n'_0$. In this Letter, the primed quantities are measured in the fluid rest (proper) frame, while the unprimed quantities are measured in the simulation frame unless otherwise specified. Densities are normalized by the initial background density $n_b$, time is normalized by the plasma frequency $\omega_{pe}\equiv(4\pi n_b e^2/m_e)^{1/2}$, velocities are normalized by the light speed $c$, and spatial scales are normalized by the inertial length $d_e\equiv c/\omega_{pe}$.

The domain size is $L_x\times L_z=384d_e \times 384d_e$ and is resolved by $3072\times6144$ cells. We load 100 macro-particles per cell for each species. The boundary conditions are periodic in the x-direction, while in the z-direction the field boundary condition is conducting and the particles are reflected at the boundaries. The half-thickness of the initial sheet is $\lambda=d_e$, $n_b=n'_0$, $T_b/m_ec^2=0.5$ and $\omega_{pe}/\Omega_{ce}=0.05$ where $\Omega_{ce}\equiv eB_{x0}/(m_e c)$ is a cyclotron frequency. The upstream magnetization parameter is $\sigma_{x0}=B_{x0}^2/(4\pi w)$ with enthalpy $w=2n'_b m_ec^2+[\Gamma/(\Gamma-1)]P'$. Here $\Gamma$ is the ratio of specific heats and $P'\equiv 2n'_b T'_b$ the total thermal pressure. For $\Gamma=5/3$ [@weinberg72a; @synge57a], $\sigma_{x0}=89$ in this run. A localized perturbation with amplitude $B_z=0.03B_{x0}$ is used to induce a dominant x-line at the center of simulation domain.

![The evolution of measured global reconnection rate $R_G$, local rate $R_L$, local inflow speed $V_{in,L}/c$ and $B_{xL}/B_{x0}$ in a plasma of $\sigma_{x0}=89$. The blue circle marks the deviation of $R_L$ from $R_G$. The grey dashed line at value $0.0$ is for reference. The orange vertical line marks the time for the analyses in Fig. \[feature\] and \[force\_balance\]. []{data-label="rate"}](rate){width="8.5cm"}

[ *Simulation results–*]{} In this Letter, we define the global reconnection rate as $R_G \equiv cE_y/(B_{x0}V_{A0})$ and the local reconnection rate as $R_L\equiv cE_y/(B_{xL}V_{AL})$. Subscripts “0” and “L” indicate quantities far from, and immediately upstream of, the diffusion region where the frozen-in condition ${\bf E}+{\bf V_e}\times {\bf B}=0$ breaks ($|z|\lesssim3.5d_e$ [@yhliu15a]). $E_y$ is the reconnection electric field at the x-line and the Alfvén speed in the relativistic limit [@sakai80a; @anile89a; @lichnerowicz67a; @yhliu14a] is $V_{A0}=c[\sigma_{x0}/(1+\sigma_{x0})]^{0.5}$ and $V_{AL}=c[\sigma_{xL}/(1+\sigma_{xL})]^{0.5}$ with $\sigma_{xL}\simeq (B_{xL}/B_{x0})^2\sigma_{x0}$. The evolution of reconnection rates are plotted in Fig. \[rate\], along with the local electron inflow speed, $V_{in,L}$, and the ratio of magnetic fields $B_{xL}/B_{x0}$. Before a quasi-steady state is reached, both the local and global rates increase as the simulation progresses. The deviation of the local rate from the global rate occurs at time $t\simeq 250/\omega_{pe}$ and $B_{xL}/B_{x0}\simeq 0.8$. $R_G$ reaches a plateau of value $\simeq 0.15$ at $t \gtrsim 300/\omega_{pe}$ while $R_L$ continues to grow and $B_{xL}/B_{x0}$ continues to drop. The local rate $R_L$ eventually reaches a plateau of value $\simeq 0.6$ and $B_{xL}/B_{x0}$ reaches a plateau of value $\simeq 0.22$ at time $t\gtrsim 600/\omega_{pe}$. The local inflow speed basically traces the local rate because of the frozen-in condition $E_y\simeq V_{in,L}B_{xL}/c$ and $V_{AL}\simeq c$ in this case, which leads $R_L=V_{in,L}/V_{AL}\simeq V_{in,L}/c$ . The values of these two quantities can approach $\sim O(1)$ with a larger $\sigma_{x0}$, as reported before [@yhliu15a; @lyutikov16a; @zenitani09a].

![The morphology of relativistic magnetic reconnection at $t=600/\omega_{pe}$. In (a), the $V_{ez}$ and a cut at $x=0$; In (b), the $|B_x|$ and a cut of $B_x$ at $x=0$. The white contour is the in-plane magnetic flux. To better illustrate the variation of the upstream field in (b), we have put an upper limit $B_{x0}$ in the color scale, which artificially reduces the $|B_x|$ around the magnetic islands at outflow exhausts.[]{data-label="feature"}](feature){width="9cm"}

To get a better idea of the spatial variation of the inflow velocity and magnetic fields at the quasi-steady state, the $V_{ez}$ and $B_x$ at time $t=600/\omega_{pe}$ are shown in Fig. \[feature\] with the in-plane magnetic flux overlaid. Immediately upstream of the intense thin current sheet, the $|V_{ez}|$ peaks at $|z|\simeq d_e$ with value $\simeq 0.65$, where $B_x$ drops to a value $\simeq 3$. Because of the thin current sheet, $d_e$-scale secondary tearing modes [@yhliu15a] are generated repeatedly, which can be seen in Fig. \[feature\]. Note that $R_G$ reaches the plateau in Fig. \[rate\] long before the generation of secondary tearing modes. The enhancement of $V_{ez}/c$ closer to the diffusion region is anti-correlated with the reduction of $B_x$ because $E_y\simeq V_{ez} B_x/c$ should be spatially uniform in a quasi-steady state under the 2D constraint, per Faraday’s law.

![In (a), the force balance in the z-direction along $x=0$ in Fig. \[feature\]; In (b), the pressure balance along $x=0$.[]{data-label="force_balance"}](force_balance){width="8.5cm"}

To get a clue of how the $B_{xL}$ drops from $B_{x0}$, we examine the force balance across the x-line at $x=0$. By combining the momentum equations for electrons and positrons [@yhliu15a; @hesse07a], the equation of force balance can be derived as $$\sum_j^{e,p}mn_j{\bf V}_j\cdot\nabla{\bf U}_j+\nabla\frac{B^2}{8\pi}+\nabla\cdot {\bf P}-\frac{{\bf B}\cdot\nabla {\bf B}}{4\pi}=-\sum_j^{e,p}mn_j\frac{\partial}{\partial t}{\bf U}_j
\label{force}$$ Here the pressure tensor ${\bf P} \equiv \sum_j^{e,p} \int d^3u {\bf v u}f_j-n_j{\bf V}_j{\bf U}_j$, and subscripts “e” and “p” stand for electrons and positrons respectively. ${\bf U} \equiv (1/n)\int d^3u {\bf u} f$ is the first moment of the space-like components of 4-velocity, and ${\bf V} \equiv (1/n)\int d^3u{\bf v} f$ as usual. On the left hand side of Eq. (\[force\]), the terms represent the inertial force, magnetic pressure gradient force, plasma thermal gradient force and magnetic tension, respectively. In the upstream region the magnetic pressure is balanced by the tension force as shown in Fig. \[force\_balance\](a). The thermal pressure is negligible because of the small plasma $\beta\equiv P/(B^2/8\pi)\simeq 0.005$. The time-derivative of inertia is negligible in the quasi-steady state. Therefore, the force balance results in a significant reduction of $B_x$ from the value far upstream at $|z| \gtrsim 150d_e$ to the value immediately upstream of the diffusion region at $|z|\simeq 3.5 d_e$ as shown in the profile of $B^2/8\pi$ in Fig. \[force\_balance\](b).

[ *Simple model–*]{} When the current sheet pinches locally, it implies a curved upstream magnetic field as illustrated in Fig. \[model\](a). The local magnetic field immediately upstream of the diffusion region, $B_{xL}$, becomes smaller than $B_{x0}$, so that the magnetic pressure gradient force balances the magnetic tension. A larger degree of localization implies a larger curvature, and a smaller $B_{xL}$, as indicated by the “line-density” of the in-plane flux in both Fig. \[model\](a), and the upstream region in Fig. \[feature\]. Hence, even though the local reconnection rate can be enhanced significantly due to the normalization, the global reconnection rate may not increase much.

To estimate this effect in the $\beta \ll 1$ limit, we analyze the force balance, $\nabla B^2/8\pi \simeq {\bf B}\cdot\nabla {\bf B}/4\pi$, at point 1 marked in Fig. \[model\](b): $$\frac{B_{x0}^2-B_{xL}^2}{8\pi\Delta z}\simeq \left(\frac{B_{x0}+B_{xL}}{2}\right)\frac{2B_z}{4\pi\Delta x}.
\label{up_force}$$ Note that $\nabla\cdot {\bf B}=0$ is also satisfied at point 1. The $B_x$ at point 1 is linearly interpolated from $B_{x0}$ and $B_{xL}$. The upstream inertial force can be formally ordered out, and it is also negligible in Fig. \[force\_balance\].

A curved upstream magnetic field naturally implies an flaring angle, and that is measured by $\Delta z/\Delta x\simeq B_z/[(B_{x0}+B_{xL})/2]$. For the proof of principle, we match it to the opening angle of the reconnection exhaust just outside of the diffusion region: $\Delta z/\Delta x\simeq B_{zL}/B_{xL}$. We obtain the relation, $$\frac{B_{zL}}{B_{xL}}\simeq \sqrt{\frac{1-B_{xL}/B_{x0}}{1+B_{xL}/B_{x0}}}.
\label{BzL_BxL}$$ This expression suggests that a larger opening angle requires a further reduction of $B_{xL}/B_{x0}$. In this sense, $B_{xL}/B_{x0}$ gauges the localization of sheet pinch. When $B_{xL}/B_{x0} \rightarrow 0$, the opening angle approaches $45^\circ$ in this model.

![The cartoons of magnetic field lines upstream of the diffusion region ($z >0$) in (a), and the geometry of reconnection in (b). The dimension of the diffusion region in (c). The predictions with $\sigma_{x0}=89$ in (d), the dashed lines use $V_{out,L}=V_{AL}$. The orange vertical line corresponds to that in Fig. \[rate\].[]{data-label="model"}](model_all){width="8cm"}

Combined with $E_y\simeq B_{zL} V_{out,L}/c$, the reconnection rates are $$R_G\simeq \left(\frac{B_{xL}}{B_{x0}}\right)\left(\frac{B_{zL}}{B_{xL}}\right)\left(\frac{V_{out,L}}{V_{A0}}\right);\ R_L\simeq \left(\frac{B_{zL}}{B_{xL}}\right)\left(\frac{V_{out,L}}{V_{AL}}\right)
\label{RG}$$ and the local inflow speed is $$V_{in,L}\simeq R_L V_{AL}.
\label{vin}$$

Using Eq. (\[BzL\_BxL\]) and the outflow speed $V_{out,L}\sim V_{AL}$, the predicted $R_G$, $R_L$ and $V_{in,L}/c$ as functions of $B_{xL}/B_{x0}$ are plotted in Fig. \[model\](d) as dashed-lines. If $B_{xL}/B_{x0}=1$, the opening angle is zero and reconnection is not expected. In the limit of $B_{xL}/B_{x0}\rightarrow 0$, the reconnecting component vanishes and reconnection ceases (i.e., $R_G=0$).

However, the geometrical constrain can reduce the outflow speed from $V_{AL}$ when the opening angle approaches $45^\circ$. This correction can be modeled through analyzing the force-balance in the x-direction at point 2 of Fig. \[model\](c): $n'mU_{out}^2/L+B_{zL}^2/8\pi L\simeq (B_{zL}/2)2(B_{xL}/2)/4\pi\delta$, where the inertial force becomes important. The outflow can be relativistic, $U_{out}\sim \gamma_{out}V_{out,L}\sim V_{out,L}^2/(1-V_{out,L}^2/c^2)$. Assuming the incompressibility of plasmas, then the aspect ratio of the diffusion region $\delta/L\sim B_{zL}/B_{xL}$, and the outflow speed becomes $$V_{out,L}\simeq c\sqrt{\frac{(1-B_{zL}^2/B_{xL}^2)\sigma_{xL}}{1+(1-B_{zL}^2/B_{xL}^2)\sigma_{xL}}}.
\label{vout}$$ This expression suggests that when $\delta/L \ll 1$ (i.e., $B_{zL}/B_{xL} \ll 1$) then $V_{out,L}\sim V_{AL}$. When $\delta/L\rightarrow 1$ (i.e., $45^\circ$), the outflow tension is balanced by the magnetic pressure and the outflow vanishes. Plugging Eqs.(\[BzL\_BxL\]) and (\[vout\]) back to Eqs. (\[RG\])-(\[vin\]), we get the solid curves in Fig. \[model\](d). This correction further constrains the reconnection rate when the opening angle is larger and $B_{xL}/B_{x0}$ is smaller.

This model suggests that during the pinching of the current sheet, a weak localization with $B_{xL}/B_{x0} \lesssim 0.9$ is enough to lead $R_G$ to $\sim 0.2$, then it varies slowly over a wide range of $B_{xL}/B_{x0}$. The local rate $R_L$ and local inflow speed $V_{in,L}/c$ can reach $\sim O(1)$ under stronger localization. The evolution of reconnection rates in Fig. \[rate\] can be qualitatively described by this model through decreasing $B_{xL}/B_{x0}$. The rates in the quasi-steady state at time $t=600/\omega_{pe}$ of Fig. \[rate\] also compares well with the prediction at $B_{xL}/B_{x0}\simeq 0.22$ with the predicted $R_G \simeq 0.14$, $R_L \simeq 0.69$ and $V_{in,L} \simeq 0.62c$. Given the simplicity of this model, this agreement is quite remarkable.

While the localization mechanism may vary in different systems, we point out a natural tendency that can lead to the $B_{xL}/B_{x0}$ reduction in such plasmas: A diffusion region sandwiched by a large $B_{xL}\simeq B_{x0}$ at $d_e$-scale (i.e., where the frozen-in condition is broken) requires the current sheet plasma to have a huge thermal pressure to balance the magnetic pressure, and a high drift speed to support the current. For instance, the initial $d_e$-scale current sheet has $T'=100m_ec^2$, $n_0\simeq 10$ and $\gamma_dV_d\simeq 10$. However, the maximum possible reconnection electric field may not be efficient enough in heating and accelerating the cold non-drifting inflowing plasma before they exit the diffusion region [@hesse11a], hence the $B_{xL}$ drops significantly until the $d_e$-scale current sheet becomes sustainable in the quasi-steady state. If this drop continues with a larger $\sigma_{x0}$, reconnection in the more extreme limit is prone to choke itself off in the quasi-steady state.

[*Discussion–*]{} Knowing the magnitude of electric field is essential for estimating the acceleration of super-thermal particles in highly magnetized astrophysical systems. This study suggests that the magnitude of the reconnection electric field is bounded by $\sim 30\%$ of the reconnecting component of magnetic field, even in the large-$\sigma_{x0}$ limit. While a weak localization of the diffusion region is required, the global reconnection rate $R_G \sim 0.1-0.3$ is not sensitive to a further increase of localization over a wide range of $B_{xL}/B_{x0}$, but the local rate and local inflow speed are. This explains the large difference between the local and global reconnection rates observed in the simulation. In this model, a larger $\sigma_{x0}$ has little effect on the profile of the global rate $R_G$, but it could make the local inflow speed closer to the speed of light [@yhliu15a]. In addition, the effect of a guide field can be included by making the relevant Alfvén speed $V_A=c[\sigma_x/(1+\sigma_x+\sigma_g)]^{0.5}$ with $\sigma_g\equiv (B_g/B_{x0})^2\sigma_{x0}$ accounting for the effect a guide field $B_g$. This expression is basically the projection of the total Alfvén speed in the outflow direction [@yhliu14a; @melzani14a; @yhliu15a]. A guide field also has little effect on $R_G$, but it significantly reduces the local inflow speed and the magnitude of the reconnection electric field through reducing the speed of Alfvénic outflows, as observed in Liu et al. [@yhliu15a]. The prediction in the non-relativistic and low-$\beta$ limit can be obtained by taking $\sigma_x \ll  1$, and the resulting $R_G$ has a slightly smaller amplitude.\
Y.-H. Liu thanks for helpful discussions with J. Dorelli. This research was supported by an appointment to the NASA Postdoctoral Program at the NASA-GSFC, administered by Universities Space Research Association through a contract with NASA. Simulations were performed with LANL institutional computing, NASA Advanced Supercomputing and NERSC Advanced Supercomputing.

[56]{} natexlab\#1[\#1]{}bibnamefont \#1[\#1]{}bibfnamefont \#1[\#1]{}citenamefont \#1[\#1]{}url \#1[`#1`]{}urlprefix\[2\][\#2]{} \[2\]\[\][[\#2](#2)]{}

, ****, ().

, ****, ().

, ****, ().

, , , , , , , ****, ().

, ****, ().

, ****, ().

, ****, ().

, , , ****, ().

, ****, ().

, , , ****, ().

, , , , ****, ().

, , , , ****, ().

, , , ****, ().

, ****, ().

, ****, ().

, , , , , ().

, , , , , ****, ().

, , , , , , , , , ****, ().

, , , , ****, ().

, , , , ****, ().

, , , , , ****, ().

, ****, ().

, , , , ****, ().

, ****, ().

, ****, ().

, , , ****, ().

, , , ****, ().

, , , , ****, ().

, , , , , ****, ().

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, , , , , ****, ().

, , , , , ****, ().

, , , ().

, in ** (), vol.  of **, pp. .

, ****, ().

, ****, ().

, in ** (), p. .

, ****, ().

, , , , , , , ****, ().

, ****, ().

, ****, ().

, , , , , , **** ().

, ** (, ), chap. , p. .

, ** (, ), chap. , p. .

, ****, ().

, ** (, ), chap. , p. .

, ** (, ), chap. , p. .

, , , , , ****, ().

, , , , ().

, , , ****, ().

, ****, ().

, , , , , ****, ().
