Reconnection and particle acceleration in 3D current sheet evolution in moderately-magnetized astrophysical pair plasma

Before describing diagnostics for reconnection, we offer explicit definitions or clarifications of some often-used terms, so that we can avoid lengthy qualifications in the text.

Transverse: the $x$ and $y$ directions, transverse to $z$, the direction of the guide magnetic field.

Transverse magnetic field energy $U_{Bt}$: the energy in the $B_{x}$ and $B_{y}$ components of the magnetic field—the volume integral of $(B_{x}^{2}+B_{y}^{2})/8\upi$ over the entire system. Even in simulations with substantial guide magnetic field, it is mainly the transverse magnetic energy (and not energy in $B_{z}$) that is depleted during reconnection.

Guide magnetic field energy $U_{Bz}$: the energy in the $B_{z}$ components of the magnetic field. Because $B_{z}(x,y,z)$ is initially uniform and the flux $\int B_{z}dxdy$ through any transverse plane is exactly conserved in the simulation, $U_{Bz}$ can only increase from its initial value (in practice it does not increase much and the increase can often be neglected).

Plasma energy: the total kinetic energy of individual plasma particles.

Magnetic energy conversion: the conversion of (transverse) magnetic field energy $U_{Bt}$ to plasma energy. Because the guide magnetic field energy $U_{Bz}$ cannot decrease, it is essentially equivalent to refer to magnetic energy conversion or transverse magnetic energy conversion.

Plasma energization: the conversion of (transverse) magnetic energy to plasma energy.

Released magnetic energy: (transverse) magnetic energy that has been converted to plasma energy.

Magnetic energy depletion: (transverse) magnetic energy conversion, with emphasis on the reduction in magnetic energy (as it is converted to plasma energy). This is nearly a synonym for “dissipation,” except that “dissipation” refers to irreversible heating, whereas “depletion” or ”conversion” also include conversion to bulk flow kinetic energy (in this paper, however, we will not measure the difference between depletion and dissipation).

Unreconnected (magnetic) field (line): a magnetic field line that crosses the entire simulation domain in the $x$ direction (possibly diagonally), without reversing (in $x$). As we discuss in §3.3, this is an imperfect, approximate, but practical definition.

“Reconnected” (magnetic) field (line): any magnetic field line that is not unreconnected (according to the definition above). In 2D reconnection, a reconnected field line is a closed loop around one or more magnetic islands or plasmoids, and our definition of a “reconnected” field line identifies most closed loops accurately. In 3D, however, there may be many field lines that are “not unreconnected” but do not resemble closed loops or even spirals around flux ropes; we nevertheless call them “reconnected” and will often use quotation marks as a reminder that they do not necessarily resemble reconnected field lines in 2D reconnection.

Upstream: the upstream region includes all points on “unreconnected” field lines. By extension, upstream plasma is the plasma in this region, upstream magnetic energy is the magnetic energy in this region, etc.

Upstream value: When we refer to specific upstream values, such as the upstream $\sigma_{h}$, $B_{0}$, $n_{be}$, or $\theta_{b}$, we mean the asymptotic or far upstream values.

Current sheet (or sheets): the region containing currents that support the reversal in magnetic field. We often refer to the initial current sheet, which is well defined, and otherwise use the term to refer to sheet-like regions where $J_{z}$ is strong.

The layer (reconnection layer, current layer): the complement of the upstream region—i.e., the region containing “reconnected” field lines. This region contains the current sheet (or current sheets) as well as reconnection outflows and plasmoids. We think of “the layer” as the evolved current sheet. Each (double-periodic) simulation contains two layers—the upper and lower layers.

Separatrix: the surface between the “upstream” region and “the layer,” i.e., that separates “reconnected” and “unreconnected” magnetic field lines. In 2D this is a smooth, well-defined surface. In 3D that may not be the case, but none of our analysis will be sensitive to the precise location of the separatrix.

Unreconnected flux: the integral of $B_{x}$ over the $x=0$ plane intersected with the “upstream” region; i.e., the integral is over all (and only) unreconnected field. See §3.3 for more detail.

Reconnected flux: the flux of magnetic field around the major (largest) plasmoid. This concept is well defined and easily measured in 2D, where the unreconnected plus reconnected flux is conserved (cf. §3.3). It is nontrivial, however, to define and measure this in 3D, and we will not do so.

Upstream magnetic energy $U_{Bt,\rm up}$: the transverse magnetic energy in the upstream region. We exclude guide field energy from this quantity, because we are primarily interested in conversion of magnetic energy to plasma energy.

Magnetic energy in the layer $U_{Bt,\rm layer}$:, the transverse magnetic energy in the layer; $U_{Bt,\rm up}+U_{Bt,\rm layer}=U_{Bt}$. We exclude guide field energy from this quantity.

Reconnection: technically, “reconnection” should involve the conversion/rearrangement of “unreconnected field lines” to “reconnected field lines” in a way that conserves the total unreconnected plus reconnected flux. In 2D reconnection, this conservation is easily verified. However, in 3D it is nontrivial to decide precisely how to characterize reconnection, much less to measure true reconnection rates. Rather than argue for any particular precise definition of reconnection, we will focus on more tangible properties of the current layer, such as magnetic energy and plasma energy. We will use the words “3D reconnection” loosely to refer to any magnetic-energy-depleting evolution of thin current sheets that, in principle, could or would undergo 2D-like reconnection. However, any discussion of reconnecting flux or field lines will refer specifically to the 2D-like flux-conserving process.

Flux annihilation: the utter disappearance of upstream (unreconnected) magnetic flux by unspecified means. In contrast, 2D reconnection does not annihilate flux (it conserves flux as described above). For example, magnetic diffusion in a plasma with finite resistivity can directly annihilate flux; although this process is far too slow to explain observed magnetic energy releases (assuming neat laminar current sheets), it could be considerably enhanced, e.g., by effective turbulent magnetic diffusivity. Incidentally, the flux $\int B_{x}dydz$ through any $x$-plane in the simulation is exactly conserved—however, by virtue of the reversing magnetic field, this flux is zero, and it tells us nothing about how much flux is reconnected or annihilated.

(Upstream) flux depletion: the depletion of the upstream magnetic flux, without regard to the process (e.g., reconnection or annihilation).

It is useful to analyse some field quantities at the current sheet “centre” (in $y$), even as the layer kinks and deforms. In the initial state, the current sheet (or layer) central surface $y_{c}(x,z,t=0)$ is the location where $B_{x}(x,y_{c},z,t=0)=0$, i.e., where $B_{x}$ crosses through zero in the $y$ direction, reversing sign. At later times, we continue to define $y_{c}(x,z,t)$ in the same way, with the following complication. The layer may become extremely distorted by the nonlinear development of the kink instability so that it folds over on itself, resulting in multiple field reversals. Because the upstream magnetic field never changes sign, $B_{x}(x,y,z,t)$ will always reverse sign at least once along $y$; however, it may reverse sign an odd number of times at each layer. If there are multiple reversals in $y$ at some $(x,z)$ and time $t$, we take $y_{c}(x,z,t)$ to be the middle reversal [e.g., if there are three reversals, $y_{c}(x,z,t)$ will be at the second reversal].

We define the displacement of the central surface, $\Delta y_{c}(x,z,t)\equiv y_{c}(x,z,t)-y_{c}(x,z,0)$.

In §5.5 we look at $\tilde{y}_{c}(k_{z},t)$, the Fourier spectrum of $\Delta y_{c}(x,z,t)$ in $z$, averaged over $x$, at a given time $t$. To find the “Fourier spectrum in $z$ averaged over $x$,” we take the power spectrum in $z$ of $\Delta y_{c}(x,z,t)$ at every $x$ (for a given time $t$), then average the power spectra over $x$, and take the square root. We normalize the result $\tilde{y}_{c}(k_{z},t)$ so that it indicates the amplitude of the sinusoidal component with $k_{z}$. E.g., if $\Delta y_{c}(x,z,t)=A\sin(k_{z}z)$, then $\tilde{y}_{c}(k_{z},t)=A$.

Reconnection breaks upstream (unreconnected) field lines and reconnects them in a different configuration in the downstream outflows. Defining and measuring unreconnected flux is straightforward in 2D using the vector potential $A_{z}(x,y)$, which is constant along field lines. In 3D, however, it is much more difficult to define unreconnected flux precisely, and so we adopt the following practical definition. An unreconnected field line is a field line that runs from $x=0$ to $x=L_{x}$ without changing direction in $x$ ($B_{x}$ never changes sign). In 2D, this definition is unambiguous and almost always agrees with the method using $A_{z}$ (disagreement could occur in the rare case of an unreconnected field line with an $\Omega$-shaped kink; we believe such cases are negligibly rare in 2D, but we have not studied whether they are also rare in 3D). In 2D, a field line that runs from $x=0$ to $x=L_{x}$ wraps around via periodic boundary conditions exactly onto itself; therefore, this definition unambiguously determines whether the entire field line has the same sign for $B_{x}$. In 3D, this is not necessarily true; nevertheless this method captures the most obviously-unreconnected upstream field lines with $|B_{y}|\ll|B_{x}|$, that run across the simulation box with little deflection in the $y$ direction.

Specifically, we trace a field line from each grid node $(0,y,z)$ in the $x=0$ plane, following it in the $+x$ direction; if the field line ever reverses direction in $x$, then the field line is considered to be “reconnected.” If the sign of $B_{x}$ does not reverse, then the field line must eventually reach $x=L_{x}$, at which point we stop and consider the field line to be “unreconnected.” The unreconnected flux is then estimated as $\sum_{y,z}B_{x}(0,y,z)\Delta y\Delta z$ where the sum is over nodes $(y,z)$ on unreconnected field lines, and $\Delta y\Delta z$ is the cell face area. We add the flux for all unreconnected field lines between the central surfaces of the two layers (cf. §3.2) to obtain the total unreconnected flux.

During this process of tracing field lines, we note the index of every cell that is penetrated by an unreconnected field line and consider any such cell to be part of the “upstream” region; any cell not penetrated by any unreconnected field line is considered part of the “layer.” In this way we can calculate, for example, the magnetic energy in upstream and layer regions. The boundary between regions of unreconnected and reconnected field lines is (an approximation of) the separatrix, which can be clearly seen in Fig. 4; (in 2D, at least) the separatrix bounds magnetic islands (O-points or plasmoids) and goes through X-points.

This analysis is rigorously based on 2D reconnection; in 3D it still paints a useful picture, even though the notions of “unreconnected” and “reconnected” are no longer precisely well-defined, and the analysis does not distinguish between reconnection and annihilation of upstream flux (i.e., between reconnected flux and annihilated flux). In 2D, any field line that is not unreconnected is reconnected, forming a closed loop in a magnetic island (when $B_{z}$ is ignored). Moreover, we find that in 2D, flux is conserved—specifically, the flux between the major O-points of the two layers (in the double-periodic simulation) is conserved. This flux can be divided into two parts: the upstream unreconnected flux between the two layers, plus the reconnected flux between the major O-point and separatrix for each layer, and the sum of these parts remains constant, at least within our measurement precision, which is better than 1 per cent of the initial flux. Therefore, in 2D, flux is not annihilated or destroyed (in a more resistive plasma, flux would be annihilated; and even in a “collisionless” PIC simulation, numerical resistivity will eventually annihilate flux, but only over extremely long times). In 3D, field lines that are “not unreconnected” do not necessarily form a closed loop around an O-point or flux rope, but again, for brevity we will still use the term “reconnected” for such field lines. Periodically we include quotation marks around “reconnected” as a reminder that it really means “not unreconnected.” In addition (spoiler alert!), we will find that in 3D, flux is not conserved in the way it is in 2D: some flux is outright annihilated.

When we refer to a “reconnection rate,” we mean the rate at which the upstream magnetic flux decays—regardless of whether this flux depletion occurs due to reconnection (as always in 2D collisionless reconnection) or to annihilation (as might be happening in 3D simulations). The reconnection rate represents the rate at which upstream flux is changing, and therefore represents an electric field along an X-line by Faraday’s law, if a suitable X-line (or reasonable approximation) exists.

We use the symbol $\beta_{\rm rec}$ for the dimensionless reconnection rate normalized to $cB_{0}$: $\beta_{\rm rec}\equiv-(cB_{0}L_{z})^{-1}d\psi/dt$, where $\psi(t)$ is the flux upstream of one layer (or half the flux between the two layers). Usually we obtain values of $\beta_{\rm rec}$ averaged over some time—for example the time over which $\psi(t)$ falls from $0.9\psi_{0}$ to $0.8\psi_{0}$, or alternatively, 0.8–0.7$\psi_{0}$ (the choice between intervals often depends on whether all simulations being compared actually reached $0.7\psi_{0}$). Sometimes we also quote $(c/v_{A})\beta_{\rm rec}$, the dimensionless reconnection rate normalized to $B_{0}v_{A}$, where $v_{A}$ is the Alfvén velocity.

Energy is a powerful diagnostic because of its physical importance, because it is conserved (to a good approximation in these simulations), and because it is generally straightforward and unambiguous to measure. We will calculate various energy components integrated over the entire simulation volume (at any given time $t$), including total particle energy $U_{\rm plasma}(t)$ and electric field energy $U_{E}(t)$. The component of most frequent interest is the magnetic field energy $U_{B}(t)$, and more specifically the transverse magnetic field energy, $U_{Bt}(t)\equiv\int dV\,(B_{x}^{2}+B_{y}^{2})/8\upi$; it is mainly $U_{Bt}$ that gets converted to particle energy $U_{\rm plasma}$, while $U_{Bz}\equiv U_{B}-U_{Bt}$ remains relatively constant (or at least negligibly small) over the course of a simulation. For zero guide field, $U_{Bt}\approx U_{B}$ and it does not much matter which we use; when comparing simulations with different guide fields, however, it is more illuminating to compare $U_{Bt}$.

We also calculate the upstream magnetic energy $U_{Bt,\rm up}$ in transverse magnetic field components of unreconnected field lines, and then define the magnetic energy in the layer (i.e., in “reconnected” field lines) to be $U_{Bt,\rm layer}=U_{Bt}-U_{Bt,\rm up}$ (see §3.3).

Sometimes the initial current sheet configuration strongly affects the time it takes reconnection to start; when comparing time evolution of simulations with different initial configurations, it is sometimes helpful to compare different cases relative to the onset time rather than $t=0$. We define the onset time as the time $t_{\rm onset}$ when the transverse magnetic energy $U_{Bt}(t)$ has declined by 1 per cent from its initial value $U_{Bt0}$. Although this value is somewhat arbitrary, no results will depend strongly on this choice, as long as it is used consistently. Occasionally the time axis of graphs will be shifted by $t_{\rm onset}$ to allow more revealing comparison between different simulations, which may have similar time evolution apart from different onset times.

We compute electron bulk velocities $\boldsymbol{v}_{e}=\boldsymbol{J}_{e}/\rho_{e}$ (in each grid cell), where $\boldsymbol{J}_{e}$ and $\rho_{e}$ are the current and charge density due to electrons; similarly, the positron velocity is $\boldsymbol{v}_{i}=\boldsymbol{J}_{i}/\rho_{i}$. The plasma velocity is then computed as the average of electron and ion velocities, $\boldsymbol{v}\equiv(\boldsymbol{v}_{e}+\boldsymbol{v}_{i})/2$. More accurately, this average should be weighted by the mass (or $\gamma m$) of particles in each cell; however, in this paper, quasineutrality and the symmetry between electrons and positrons make the simple average a reasonable approximation.

Particle energy distributions $f(\gamma)$ are shown integrated over the entire simulation box at single time snapshots (where the Lorentz factor $\gamma$ is a proxy for particle energy $\gamma m_{e}c^{2}$). The energy spectra often display power laws $f(\gamma)\sim\gamma^{-p}$ extending to high energies (well above the average energy).

We determine the power-law index $p$ by finding the longest, straightest segment on a log-log plot (Werner et al., 2018). We compute the local slope, $p(\gamma)=-d\ln f/d\ln\gamma$ and search exhaustively for the interval $[\gamma_{1},\gamma_{2}]$ with the largest $\gamma_{2}/\gamma_{1}$ such that $p(\gamma)$ remains within a range of $\Delta p$ over the entire interval. We then choose “the” power-law index $p$ to be the median $p(\gamma)$ over $[\gamma_{1},\gamma_{2}]$. We do this separately for $\Delta p=$0.1, 0.2, and 0.4, as well as using spectra at different nearby times, considering variation in $p$—whether due to time variation or choice of $\Delta p$—as uncertainty in the measurement. In this paper, we display “error bars” on $p$ comprising the middle 68 per cent of all the values of $p$ measured.

To find the cutoff of the high-energy power law, we fit $f(\gamma)$ over $[\gamma_{1},\gamma_{2}]$ to the form $A\gamma^{-p}$, where $p$ is determined as above (thus only the normalization $A$ must be found). We then consider the cutoff $\gamma_{c}$ to be the energy at which $f(\gamma_{c})=e^{-1}A\gamma_{c}^{-p}$.

We begin by reviewing the familiar behaviour of plasmoid-dominated reconnection in 2D, for a single representative simulation with $L_{x}=1280\sigma\rho_{0}$; other initial parameters are: $B_{gz}=0$, $\eta=5$, $\beta_{d}=0.3$, $\delta=(2/3)\sigma\rho_{0}$, $a=0$. This simulation exhibits the familiar behaviour of 2D plasmoid-dominated reconnection—it is qualitatively similar to plasmoid-dominated reconnection at small and large $\sigma_{h}$. The simulation lingers in its initial state until the tearing instability triggers reconnection, resulting in a chain of plasmoids (magnetic islands or O-points) in each layer. (In the following, we describe just one of the layers in the double-periodic system; both layers behave qualitatively similarly.) Reconnection occurs at X-points in thin, elementary current sheets between plasmoids: at X-points, upstream magnetic field lines are broken, and reconnected in a different topology, with reconnected field lines wrapping around plasmoids in the reconnection outflows. In this process, some of the upstream magnetic energy is converted to particle/plasma energy, while some of it ends up as magnetic energy in the layer (i.e., in “reconnected” field in plasmoids). Secondary tearing breaks up the elementary current sheets when they become too long and too thin, detaching the small plasmoids from the X-point that fed them magnetic energy and reconnected flux. Thereafter, those plasmoids no longer grow via reconnection, but instead grow by merging with other plasmoids via the coalescence instability, which conserves the flux around O-points. A hierarchy of different-sized plasmoids thus develops, and eventually a single monster plasmoid dominates the simulation, continuing to grow as it consumes (i.e., merges with) smaller plasmoids. If $L_{y}/L_{x}$ is large enough, this major plasmoid eventually grows to a size of order $L_{x}$, and reconnection can no longer continue; the angle of the magnetic separatrix at the major X-point opens to 90 degrees, and reconnection effectively stops (actually, systems often oscillate with low amplitudes about the final state, with magnetic energy and flux sloshing between the major plasmoid and the upstream). If $L_{y}\gg L_{x}$ (even if $L_{y}\gtrsim 2L_{x}$), there will still be significant upstream “unreconnected” magnetic field that could in principle be reconnected were it not prevented from doing so by reaching a stable magnetic configuration (at least on reconnection timescales; magnetic diffusion might deplete additional magnetic energy on much, much longer timescales).

During this evolution, upstream magnetic flux is reconnected and some upstream magnetic energy is converted into plasma energy, while some ends up permanently stored in the magnetic field of the major plasmoid. To describe this process more quantitatively, we will rely heavily on several important diagnostics. These will be essential throughout the rest of the paper, as we explore reconnection over a very large parameter space, to help us compare and contrast reconnection with different parameters including $a$, $\eta$, $\beta_{d}$, $L_{x}$, $B_{gz}$, and (in 3D) $L_{z}$. We will measure, for example, the amount of magnetic flux that is reconnected over time (cf. §3.3). However, perhaps of even more importance is the behaviour of various energy components—for example, the plasma or magnetic energy versus time (cf. §3.4). In addition, we further decompose the plasma energy into the energy distribution of particles to distinguish thermal heating from NTPA (cf. §3.7). In the following, we describe these diagnostics for the single representative simulation.

The behaviour of magnetic flux and energy is qualitatively similar to that observed in other 2D regimes, such as $\sigma_{h}\gg 1$ (e.g., Guo et al., 2014; Werner et al., 2016; Werner & Uzdensky, 2017; Werner et al., 2018). Figure 1 shows unreconnected flux versus time as well as the energy (and energy changes) versus time in various components, including the energy in magnetic fields and particles. The unreconnected (upstream) magnetic flux $\psi(t)$ decreases over the course of reconnection, finally down to about 55 per cent of its initial value (i.e., reconnecting 45 per cent of upstream flux). (We note that these values are specific to this representative simulation; reconnection with other parameters may yield different values while being qualitatively similar.) During this time, the magnetic energy falls to about 70 per cent of its initial value, with the 30 per cent loss converted almost entirely to plasma energy. This occurs within a time of roughly 8–12$\>L_{x}/c$ over which the rate of reconnection continually slows (in this closed system). The dimensionless reconnection rates (cf. §3.3), averaged over the time it takes for unreconnected flux $\psi(t)$ to fall from $0.9\psi_{0}$ to $0.8\psi_{0}$, and from $0.8\psi_{0}$ to $0.7\psi_{0}$ (as measured in Werner et al., 2018), are $\beta_{\rm rec}\equiv-(cB_{0}L_{z})^{-1}\dot{\psi}=0.032$ and $0.22$, respectively (here the dimensionless rate is normalized to $B_{0}c$). Normalized to $B_{0}v_{A}$ instead of $B_{0}c$, where $B_{0}v_{A}=B_{0}c\sqrt{\sigma_{h}/(1+\sigma_{h})}=B_{0}c/\sqrt{2}$, the reconnection rates over these same intervals are $(c/v_{A})\beta_{\rm rec}=0.045$ and $0.031$. The normalized-to-$B_{0}v_{A}$ reconnection rate around 0.03–0.05 is on the order of, but less than the oft-cited nominal value of 0.1 for fast reconnection (Cassak et al., 2017)—a value that is realized for ultrarelativistic reconnection with $\sigma_{h}\gg 1$.

By way of comparison, we note that, with the same aspect ratio $L_{y}/L_{x}=2$, a similar setup with $\sigma_{h}\gg 1$ depletes roughly 40 per cent of magnetic field energy (Werner & Uzdensky, 2017) and reconnects 55 per cent of the initial flux (Werner et al., 2018). Energy conversion in high-$\sigma_{h}$ reconnection occurs within about 3$\>L_{x}/c$, with normalized reconnection rates around 0.10–0.12.

An important aspect of 2D reconnection is that the upstream field lines are indeed broken and reconnected in the following precise sense: the decrease in upstream flux equals the increase in flux around the O-points in magnetic islands (plasmoids). In 2D simulations we can measure these fluxes accurately using the $z$-component of the magnetic vector potential (see §3.3), and we find that the sum of these fluxes—which is precisely equivalent to the total flux between the major O-points in each layer—is conserved to better than 1 per cent. This is not surprising; the only way the total flux can change is (by Faraday’s law) if $E_{z}\neq 0$ at the major O-points (e.g., $E_{z}$ could be non-zero due to resistivity, which would allow annihilation of magnetic flux via magnetic diffusion, but typically the rate of annihilation is very slow, especially in a collisionless plasma).

The reconnection rate measurements quoted in this paper are defined in terms of an upstream field of $B_{0}$, although in these closed systems, the actual (far) upstream field diminishes over time (although not very much for large systems). For simulations presented in this paper, the upstream magnetic field decays only to $\approx 0.9B_{0}$, which affects the calculation less than the measurement uncertainty. There are also questions about whether “the upstream field” should be far upstream or at some upstream point closer to the current sheet (Liu et al., 2017; Cassak et al., 2017), and the difference between far upstream and near upstream might potentially depend on $\sigma_{h}$. However, we must leave a better understanding of the precise $\sigma_{h}$-dependence of the 2D reconnection rate to future work.

When $\sigma_{h}$ is not large, particles cannot gain much energy on average during reconnection. E.g., for $\sigma_{h}=1$ and 30 per cent of magnetic energy depleted, the average particle energy (averaged over the entire simulation) increases by a relative fraction of only $\Delta\gamma/\gamma_{b}=0.30(2/3)\sigma_{h}=0.20$ (where $\gamma_{b}=3\theta_{b}=0.75\sigma$ is the initial average Lorentz factor). Nevertheless, some particles reach energies well above the average energy, $\bar{\gamma}\approx 0.7$–$0.9\sigma$, as shown by the particle energy spectra in Fig. 2. Reconnection in the $\sigma_{h}=1$ regime clearly generates NTPA; a high-energy nonthermal power-law $f(\gamma)\sim\gamma^{-p}$ becomes apparent within a time of 1–2$L_{x}/c$, with approximately the same power-law index $p\approx 4$ as at later times; the power law extends up to a cutoff energy exceeding $4\sigma$ (cf. Werner et al., 2016). The power-law index of $p\approx 4$ is quite steep; to distinguish small changes (when comparing other simulations), we will usually graph $\gamma^{4}f(\gamma)$, so that if $p$ were exactly 4, the graph would be a horizontal line. The high-energy cutoff of the power law has been observed, for simulations with high $\sigma_{h}\gg 1$ and very large $L_{x}$, to grow rapidly to an energy $\sim 4\sigma$ (Werner et al., 2016), and then to continue to grow slowly (sublinearly) in time (Petropoulou & Sironi, 2018; Hakobyan et al., 2021); our results are consistent with this.

Importantly, the high-energy cutoff measured in our simulations is far below the extreme particle acceleration limit for this simulation size and duration—i.e., below the energy a particle could gain in the nominal reconnection electric field of $E_{\rm rec}$ over the simulation duration $T$ (which is usually proportional to simulation size $L_{x}$). Such an “extremely-accelerated” particle would gain maximum energy $m_{e}c^{2}\Delta\gamma_{\rm ext}=eE_{\rm rec}cT$. We can estimate $\Delta\gamma_{\rm ext}$ using $\beta_{\rm rec}$ to relate $E_{\rm rec}$ and $B_{0}$, i.e., $E_{\rm rec}=\beta_{\rm rec}B_{0}$. Since $m_{e}c^{2}/eB_{0}\equiv\rho_{0}$, $\Delta\gamma_{\rm ext}\sim\beta_{\rm rec}cT/\rho_{0}\sim 0.02(cT/L_{x})L_{x}/\rho_{0}$. In our simulations, $\beta_{\rm rec}$ decreases over time (as reconnection slows to a stop), but above we measured $\beta_{\rm rec}\approx 0.02$ in the middle of active reconnection; and in this simulation, $L_{x}=1280\sigma\rho_{0}$, so $\Delta\gamma_{\rm ext}\sim 26(cT/L_{x})\sigma$. Thus, after $T\sim 10L_{x}/c$, the maximum energy gain would be around $\Delta\gamma_{\rm ext}\sim 250\sigma$. However, rather than estimate $\beta_{\rm rec}$ and the active reconnection time $T$, it is convenient to consider that by Faraday’s law, $\int_{0}^{T}E_{\rm rec}dt=-\Delta\psi/(cL_{z})$ where $\Delta\psi$ is the change in unreconnected upstream flux; in a 2D simulation, the simulation size $L_{z}$ in the unsimulated third dimension is arbitrary, but the upstream flux is also proportional to $L_{z}$—the flux upstream of the initial current sheet is $\psi_{0}\approx B_{0}L_{z}L_{y}/4$. Therefore, $\Delta\gamma_{\rm ext}\approx(L_{y}/4L_{x})(\Delta\psi/\psi_{0})(eB_{0}L_{x}/m_{e}c^{2})=0.5(0.45)(L_{x}/\rho_{0})$, or $\Delta\gamma_{\rm ext}\approx 0.2(1280\sigma)\approx 250\sigma$.

In this simulation, we observe $\Delta\gamma\lesssim 40\sigma$, far below the extreme acceleration limit, $\Delta\gamma_{\rm ext}\approx 250\sigma$. However, some particles undergo extreme acceleration up to at least $\gamma\sim 4\sigma$ (but not beyond $\gamma\approx 20\sigma$) in a time consistent with constant acceleration in $E_{\rm rec}$. We can infer from Fig. 2 that a tiny fraction of particles reach nearly $\gamma\approx 20\sigma$ by time $0.44L_{x}/c$. At this time, $\Delta\psi/\psi_{0}=0.039$, so we estimate $\gamma_{\rm ext}=0.5(0.039)(1280\sigma)=25\sigma$, close to the observed $\gamma\approx 20\sigma$ (which we emphasize is not the power-law cutoff, but the maximum energy gained by any particle in the simulation). This suggests that particles are in fact “extremely-accelerated” up to some energy less than but on the order of $\gamma\approx 20\sigma$. However, no particle ever exceeds $\gamma\approx 40$ significantly, even after much longer times. We note that: (1) we measured $\Delta\psi$ with a time cadence of $0.44L_{x}/c$, so we cannot examine extreme acceleration on smaller timescales, and (2) ideally we would prefer to know about extreme acceleration in the middle of the simulation rather than at the beginning, but we present this result for the very beginning because it is the only time that we can know that a particle with $\gamma\approx 20$ accelerated to that energy within $0.44L_{x}/c$. Thus our results are consistent with particles experiencing extreme acceleration up to an energy around $\gamma\sim 4\sigma$ (Werner et al., 2016), and experiencing slower acceleration thereafter (as described by Petropoulou & Sironi, 2018; Hakobyan et al., 2021).

Incidentally, we also note that all particle gyroradii are much smaller than the system size. A particle with energy $\gamma$ has gyroradius $\sim\gamma\rho_{0}$ (in field $B_{0}$), and so a particle would need $\gamma=640\sigma$ for its gyroradius to equal $L_{x}/2$; therefore, the system size is much larger than the gyro-orbits of even the most energetic particles.

In the rest of this paper, we will be comparing the flux reconnected versus time, as well as the conversion of magnetic to plasma energy and the resulting particle energy spectra, for reconnection simulations with different initial parameters (including 3D simulations with different lengths $L_{z}$). In the following subsections, we consider 2D simulations all with the same upstream plasma conditions as in this subsection, but with different initial current sheet configurations (§4.2), and different system sizes (§4.3); then, we consider the addition of different guide field strengths (§4.4). We will find that for large systems, the initial current sheet configuration and system size $L_{x}$ have relatively weak effects on reconnection, including on energy conversion and NTPA; in contrast, (sufficiently strong) guide magnetic field will significantly slow energy conversion and inhibit NTPA.

In this subsection we show that 2D reconnection behaviour is largely (but not completely) determined by the background (upstream) plasma and not the configuration of the initial current sheet. We will compare results from 2D simulations with varying initial magnetic perturbation strength $a$, and varying initial current sheet parameters (i.e., density $n_{d0}$ or, equivalently, $\eta=n_{d0}/n_{b0}$, temperature $\theta_{d}$, drift speed $\beta_{d}c$, and half-thickness $\delta$). All simulations in this subsection have identical background plasma (i.e., $\sigma=10^{4}$, $\sigma_{h}=1$, hence upstream $\beta_{\rm plasma}=0.5$, and $B_{gz}=0$, as well as $L_{y}/L_{x}=2$).

The initial magnetic field perturbation $a$ [cf. Eq. (4)], can be varied independently of all other parameters, as in the following substudy (2D-a); however, increasing $a$ shifts the system away from the Harris equilibrium. In contrast, we always maintain the initial equilibrium when altering $\eta$, $\theta_{d}$, $\beta_{d}$, and $\delta$. The equilibrium satisfies pressure balance, or $\eta\theta_{d}/\gamma_{d}=\sigma/2$, and Ampere’s law, which requires $\delta=\sigma\rho_{0}/(\eta\beta_{d})$. Thus, the initial current sheet, though fully specified through four parameters, has only two independent degrees of freedom in this study; usually we specify $\eta$ and $\beta_{d}$. For reference we have listed the relevant parameters corresponding to different values of $\eta$ and $\beta_{d}$ in table 1.

It is useful to compare the sheet half-thickness $\delta$ with the characteristic gyroradii of the upstream plasma and the drifting (initial current sheet) plasma. The characteristic gyroradius of the upstream plasma is $\rho_{b}=3\theta_{b}\rho_{0}=(3/4)\sigma\rho_{0}\sim\sigma\rho_{0}$; the average gyroradius of the drifting plasma is $\rho_{d}=3\theta_{d}\rho_{0}=3\gamma_{d}/(2\eta)$. The ratio of sheet thickness to upstream plasma scales is $\delta/\sigma\rho_{0}=1/(\eta\beta_{d})$, and the ratio to drifting plasma scales is $\delta/\rho_{d}=2/(3\gamma_{d}\beta_{d})$. In substudy (2D-b), below, we vary $\eta$ while keeping $\beta_{d}$ constant, which changes $\delta$ with respect to upstream scales but not with respect to drifting plasma scales. On the other hand, in substudy (2D-c), we vary $\beta_{d}\propto\eta^{-1}$, changing $\delta/\rho_{d}$ but leaving $\delta/\sigma\rho_{0}$ constant.

(2D-a) Varying the initial perturbation. For $a=0$, the magnetic field lines are initially entirely parallel to $x$. Increasing $a$ perturbs the field so that there is a magnetic island and an X-point in the layer. As described in §2, when $a=1$, the initial island half-height $s$ (the maximum height of the separatrix above the initial midplane) equals the current sheet half-thickness, $s=\delta$; for $a\lesssim 1$, the separatrix that bounds the island is contained within the initial layer. We investigate the effect of varying $a$ over a wide range, for two different initial current sheets: the first [$\beta_{d}=0.3$, $\eta=5$, $\delta=(2/3)\sigma\rho_{0}$, $\theta_{d}=0.10\sigma$] was chosen to be like most simulations in this paper, marginally resolving the initial sheet with $\delta/\Delta x=2$; the second [$\beta_{d}=0.3$, $\eta=1$, $\delta=(10/3)\sigma\rho_{0}$, $\theta_{d}=0.52\sigma$] was chosen to resolve the current sheet with $\delta/\Delta x=10$. In both cases we used modest system sizes with $L_{x}=320\sigma\rho_{0}$. To help distinguish between the effect of $a$ and stochastic variation, we ran two simulations for each $a$, identical except for random initial particle velocities (the source of randomness is discussed later in detail in §5.2).

For the denser, thinner sheet ($\beta_{d}=0.3$, $\eta=5$, used in most simulations here), we ran simulations with $a\in\{$0, 0.22, 0.55, 1.1, 2.2, 5.5, 11, 22, 55, 110$\}$, corresponding to $s/\delta\in\{$0, 0.44, 0.72, 1.1, 1.6, 3.1, 5.4, 10, 23, 43$\}$, respectively. For $a<0.3$, we note that $s<\Delta x$, so the simulation can hardly tell the difference between $a=0$ and $a=0.22$; and for $a=110$, the separatrix extends out to $s=43\delta=29\sigma\rho_{0}\approx 0.1L_{x}$. In Fig. 3(a) we observe that the evolution of magnetic energy versus time is nearly the same for the two simulations with $a<0.3$; it is also the same, up to stochastic variation, within the group of simulations with $0.55\lesssim a\lesssim 22$. The simulations with $a<0.3$ exhibit slightly slower magnetic energy conversion for $1\lesssim tc/L_{x}\lesssim 4L_{x}/c$ than the simulations with $0.55\lesssim a\lesssim 22$, but by $t\gtrsim 5L_{x}/c$ they have all converted the same amount of energy, so while noticeable, this is a very minor difference. For $a\gtrsim 55$ (very large perturbations), we observe an increasingly quick onset of reconnection and more rapid magnetic depletion. The particle energy spectra in Fig. 3(b) are all very similar [considering that differences are enhanced by showing $\gamma^{4}f(\gamma)$ instead of $f(\gamma)$], although very large perturbations result in an almost imperceptibly harder spectrum with a slightly smaller fraction of particles reaching the very highest energies. We conclude that for all but very large perturbations, the initial perturbation does not have a very substantial effect on energy conversion or NTPA.

For the less dense, initially-thick current sheet, [$\beta_{d}=0.3$, $\eta=1$, $\delta=(10/3)\sigma\rho_{0}$, $\theta_{d}=0.52\sigma$], we explored $a\in\{$0, 0.04, 0.11, 0.22, 0.44, 1.1, 2.2, 4.4, 11, 22$\}$, corresponding to separatrix heights $s/\delta\in\{$0, 0.2, 0.31, 0.44, 0.63, 1, 1.6, 2.5, 5.1, 9$\}$. The results are shown in Fig. 3(c,d). Here, $s>\Delta x$ as long as $a>0.011$ (i.e., for all the nonzero $a$ used); for $a=22$, the separatrix extends to $s=9\delta=30\sigma\rho_{0}\approx 0.1L_{x}$. The less dense, thicker sheet leads to a slower onset of reconnection; increasing the initial perturbation strength tends to hasten reconnection onset. However, once reconnection is triggered, the energy evolves very similarly for zero and large perturbations, as seen in Fig. 3(c), which shows magnetic energy versus time, shifted relative to $t_{\rm onset}$, the time at which the magnetic energy has fallen by 1 per cent (cf. §3.5). With this time shift, the magnetic energy curves are almost identical in all cases except for the largest perturbation, $a=22$. Not surprisingly, the particle spectra in Fig. 3(d) are similarly identical (within stochastic variation). We again conclude that the strength of perturbation (at least for $a\lesssim 20$) has negligible effect on energy conversion and NTPA—aside from having a strong effect on the time to reconnection onset.

For both overdensities $\eta=1$ and $\eta=5$, the initial perturbation has significant consequences for the evolution of the plasmoid hierarchy, even though this does not seem to alter the overall energy conversion rate and NTPA. Without perturbation (Fig. 4, left), a number of roughly equal-sized plasmoids form at the onset of reconnection, with sizes determined by the fastest-growing wavelength of the tearing instability. These plasmoids undergo the coalescence instability and merge in pairs, resulting in a smaller number of larger plasmoids, each of which contains a similar fraction of the initially-drifting particles that made up the initial current sheet. In contrast, even a very small initial perturbation will favour one plasmoid (the “major” plasmoid), which is always larger than the others (Fig. 4, right); it may contain many or most of the initially-drifting particles. In this case, subsequent plasmoid mergers tend to be between different-sized plasmoids; plasmoids form at the major X-point, and move toward the major plasmoid (with mergers between smaller plasmoids sometimes occurring before they reach the major plasmoid). We note that the final state has just a single monster plasmoid, regardless of initial perturbation.

Summarizing (2D-a), we see that (for $L_{x}\gtrsim 320\sigma\rho_{0}$) the initial perturbation (if not terribly large, i.e., $a\lesssim 20$, or $s/\delta\lesssim 9$) does not significantly alter global energy evolution (after reconnection onset) or NTPA. However, the initial perturbation can reduce the time it takes reconnection to start, and it can also have a significant effect on the nature of plasmoid formation and subsequent evolution. Similar differences in plasmoid evolution caused by artificial triggering with a locally reduced pressure have been previously studied in 2D semirelativistic electron-ion reconnection (Ball et al., 2019); there NTPA was found to be similarly insensitive to triggering, at least for weak guide field.

(2D-b) Varying current sheet density $\eta$ with fixed $\beta_{d}=0.3$; i.e., varying $\delta/\sigma\rho_{0}$ for fixed $\delta/\rho_{d}=2.1$. We now look at the effects of different initial current sheets, keeping $\beta_{d}=0.3$ constant and varying $\eta\in\{$0.5, 1, 3.1, 5, 10$\}$ (considering $\eta<1$ is unusual in reconnection research, but we include $\eta=0.5$ to show that it is not that different from $\eta=1$, although it does start to indicate a trend for the limit $\eta\rightarrow 0$). Table 1 shows how other quantities such as $\delta$ vary as $\eta$ is varied in this study. Increasing $\eta$ directly increases the current sheet density $n_{d0}$ (relative to the upstream plasma density $n_{b0}$, which remains the same). As $\eta$ increases, the temperature $\theta_{d}$ decreases to maintain pressure balance across the current sheet against the upstream magnetic field, $\theta_{d}=\gamma_{d}\sigma/(2\eta)\approx 0.52\sigma/\eta$; this reduces the fundamental length scales (and hence also time scales) of the current sheet. (We note, for context, that if background plasma with $n_{b0}$ and $\theta_{b}$ were adiabatically compressed with adiabatic index 4/3, appropriate for relativistically-hot plasma, to a density $n_{\rm ad}$ and temperature $\theta_{\rm ad}$ sufficient to balance the upstream magnetic pressure $B_{0}^{2}/8\upi$ as well as the upstream plasma pressure $n_{b0}\theta_{b}m_{e}c^{2}$, this would yield $\eta_{\rm ad}\equiv n_{\rm ad}/n_{b0}=2.28$ and $\theta_{\rm ad}/\theta_{b}=1.3$.)

As $\eta$ varies (with fixed $\beta_{d}$), we have $\theta_{d}\propto\eta^{-1}$ and so $\delta$ as well as the fundamental length scales of the current sheet plasma (e.g., both the gyroradius $\rho_{d}$ and collisionless skin depth $d_{ed}$) all scale $\propto\eta^{-1}$ (as long as $\gamma_{d}\approx 1$). Therefore, $\delta/\rho_{d}=2.1$ and $\delta/d_{ed}=2.7$ remain constant (see table 1), although the current sheet thickness changes with respect to upstream plasma scales (e.g., $\rho_{b}=0.75\sigma\rho_{0}$) as well as $L_{x}$ and $\Delta x$.

For this investigation we use a larger system size, $L_{x}=1280\sigma\rho_{0}$, so that $L_{x},L_{y}\gg\delta$ for all $\delta$ explored; other fixed parameters are: $\beta_{d}=0.3$ and $B_{gz}=0$. We used the exact same initial magnetic field $\boldsymbol{B}(x,y)$ in all cases, with the same small perturbation: $a\delta/L_{x}=0.0052$ — the dependence $a\sim\delta^{-1}$ ensures that the initial field is independent of $\delta$—i.e., $s/L_{x}$ is fixed [cf. Eq. (4)]. Based on (2D-a) above, we expect our conclusions for this subsection would be the same for other small perturbations, including zero perturbation.

Figure 5 shows (left) the magnetic energy evolution and (right) NTPA spectra for simulations with a range of $\eta=0.5$–$10$, corresponding to $\delta/\sigma\rho_{0}=6.7$–$0.33$. We first note, as above, that reconnection onset is delayed by less dense, thicker current sheets, but that once magnetic energy depletion begins, energy evolution and NTPA are very similar and almost independent of $\eta$. However, closer examination shows a tentative trend of decreasing NTPA efficiency with decreasing $\eta$; this trend is very small for $\eta>1$, but becomes noticeable for the $\eta=0.5$ case (with $\theta_{d}/\theta_{b}=4.2$ and $\delta\approx 7\sigma\rho_{0}$), which reconnects a little slower and generates slightly less efficient NTPA. Specifically, $\eta=0.5$ yields a power law $p=4.6\pm 0.1$, whereas denser layers with $\eta\geq 1$ yield $p=4.4\pm 0.1$. Also, denser initial sheets seem to lead to slightly more overall magnetic energy depletion. This weak $\eta$-dependence hints at a more severe behaviour in the limit of underdense current sheets: in a simulation (not shown here) with $\eta=0.1$ and $\delta=33\sigma\rho_{0}$, reconnection did not even start, despite a simulation duration longer than $10L_{x}/c$—hence no magnetic energy was depleted.