3.1. Terminology

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{\rm \pi}$ 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 \,{\textrm {d}x}\,{\textrm {d}y}$ 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 2-D 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 three dimensions, 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 2-D 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 two dimensions this is a smooth, well-defined surface. In three dimensions 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 two dimensions, where the unreconnected plus reconnected flux is conserved (cf. § 3.3). It is non-trivial, however, to define and measure this in three dimensions, 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 2-D reconnection, this conservation is easily verified. However, in three dimensions it is non-trivial 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 ‘3-D reconnection’ loosely to refer to any magnetic-energy-depleting evolution of thin current sheets that, in principle, could or would undergo 2-D-like reconnection. However, any discussion of reconnecting flux or field lines will refer specifically to the 2-D-like flux-conserving process.

Flux annihilation: the utter disappearance of upstream (unreconnected) magnetic flux by unspecified means. In contrast, 2-D 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\,{\textrm {d}y}\,\textrm {d}z$ 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).

3.2. Current sheet (layer) central surface

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$.

3.3. Unreconnected flux and reconnection rate

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 two dimensions using the vector potential $A_z(x,y)$, which is constant along field lines. In three dimensions, 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 two dimensions, 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 two dimensions, but we have not studied whether they are also rare in three dimensions). In two dimensions, 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 three dimensions, 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 figure 4; (in two dimensions, at least) the separatrix bounds magnetic islands (O-points or plasmoids) and goes through X-points.

This analysis is rigorously based on 2-D reconnection; in three dimensions 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 two dimensions, 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 two dimensions, 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 two dimensions, 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 three dimensions, 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 three dimensions, flux is not conserved in the way it is in two dimensions: 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 2-D collisionless reconnection) or to annihilation (as might be happening in 3-D 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 _\textrm {rec}$ for the dimensionless reconnection rate normalized to $cB_0$: $\beta _\textrm {rec} \equiv -(cB_0L_z)^{-1}\,\textrm {d}\psi /\textrm {d}t$, where $\psi (t)$ is the flux upstream of one layer (or half the flux between the two layers). Usually we obtain values of $\beta _\textrm {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 _\textrm {rec}$, the dimensionless reconnection rate normalized to $B_0 v_A$, where $v_A$ is the Alfvén velocity.

3.4. Energy

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_\textrm {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 \, \textrm {d}V\, (B_x^2+B_y^2)/8{\rm \pi}$; it is mainly $U_{Bt}$ that gets converted to particle energy $U_\textrm {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).

3.5. Onset time $t_\textrm {onset}$

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_\textrm {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_\textrm {onset}$ to allow more revealing comparison between different simulations, which may have similar time evolution apart from different onset times.

3.6. Bulk Velocity

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.

3.7. Particle energy spectra and power-law fitting

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. Reference Werner, Uzdensky, Begelman, Cerutti and Nalewajko2018). We compute the local slope, $p(\gamma )=-\textrm {d}\ln f/\textrm {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}$.

4.1. Overview of 2-D reconnection with $\sigma _h=1$

We begin by reviewing the familiar behaviour of plasmoid-dominated reconnection in two dimensions, 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 2-D 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$^{\circ }$, 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 time scales; magnetic diffusion might deplete additional magnetic energy on much, much longer time scales).

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 three dimensions) $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 vs 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 2-D regimes, such as $\sigma _h\gg 1$ (e.g. Guo et al. Reference Guo, Li, Daughton and Liu2014; Werner et al. Reference Werner, Uzdensky, Cerutti, Nalewajko and Begelman2016; Werner & Uzdensky Reference Werner and Uzdensky2017; Werner et al. Reference Werner, Uzdensky, Begelman, Cerutti and Nalewajko2018). Figure 1 shows unreconnected flux vs time as well as the energy (and energy changes) vs 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 approximately 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 approximately 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. Reference Werner, Uzdensky, Begelman, Cerutti and Nalewajko2018), are $\beta _\textrm {rec}\equiv -(cB_0 L_z)^{-1} \dot {\psi }=0.032$ and $0.022$, 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 _\textrm {rec}=0.045$ and $0.031$. The normalized-to-$B_0 v_A$ reconnection rate around 0.03–0.05 is of the order of, but less than the oft-cited nominal value of 0.1 for fast reconnection (Cassak, Liu & Shay Reference Cassak, Liu and Shay2017) – a value that is realized for ultrarelativistic reconnection with $\sigma _h \gg 1$.

Figure 1.

Panel (a) shows the unreconnected flux (normalized to the initial unreconnected flux $\psi _0$) diminishing over time as reconnection occurs. Horizontal dashed lines mark the levels $\psi =0.9\psi _0$, $0.8\psi _0$ and $0.7\psi _0$; the dimensionless reconnection rates averaged over each of these $0.1\psi _0$ drops are $\beta _\textrm {rec}=$0.03 and 0.02 (cf. § 3.3), and dashed lines with those slopes are drawn. Panel (b) shows the energy in various components vs time (normalized to $\int \,\textrm {d}V \, B_0^2/8{\rm \pi}$), and panel (c) shows the changes in those components from their initial values. Energy components shown are total energy ($U_\textrm {total}$, dash-dotted black line), magnetic ($U_B$, solid blue), plasma/particle ($U_\textrm {plasma}$, dashed green) and electric ($U_E$, dotted red – in panel (b), $U_E$ is magnified by a factor of 20).

By way of comparison, we note that, with the same aspect ratio $L_y/L_x=2$, a similar set-up with $\sigma _h\gg 1$ depletes roughly 40 per cent of magnetic field energy (Werner & Uzdensky Reference Werner and Uzdensky2017) and reconnects 55 per cent of the initial flux (Werner et al. Reference Werner, Uzdensky, Begelman, Cerutti and Nalewajko2018). Energy conversion in high-$\sigma _h$ reconnection occurs within approximately 3$\:L_x/c$, with normalized reconnection rates around 0.10–0.12.

An important aspect of 2-D 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 2-D 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 (Cassak et al. Reference Cassak, Liu and Shay2017; Liu et al. Reference Liu, Hesse, Guo, Daughton, Li, Cassak and Shay2017), 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 2-D reconnection rate to future work.

When $\sigma _h$ is not large, particles cannot gain much energy on average during reconnection. For example, 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 figure 2. Reconnection in the $\sigma _h=1$ regime clearly generates NTPA; a high-energy non-thermal 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. Reference Werner, Uzdensky, Cerutti, Nalewajko and Begelman2016). 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. Reference Werner, Uzdensky, Cerutti, Nalewajko and Begelman2016), and then to continue to grow slowly (sublinearly) in time (Petropoulou & Sironi Reference Petropoulou and Sironi2018; Hakobyan et al. Reference Hakobyan, Petropoulou, Spitkovsky and Sironi2021); our results are consistent with this.

Figure 2.

The particle energy spectra $f(\gamma )$ at logarithmically spaced times, compensated by $\gamma ^2$ to show the relative amount of energy stored in different energy ranges (a) and by $\gamma ^4$ to enhance differences (b). Here, particle energy is $\gamma m_e c^2$. The average particle energy $\bar {\gamma }$ is indicated at $t=0$ and $tc/L_x=9.3$. A power law $\gamma ^{-4.4}$ is also shown (in panel (b), $f\sim \gamma ^{-4}$ would be horizontal).

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_\textrm {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_ec^2 \Delta \gamma _\textrm {ext} = eE_\textrm {rec}cT$. We can estimate $\Delta \gamma _\textrm {ext}$ using $\beta _\textrm {rec}$ to relate $E_\textrm {rec}$ and $B_0$, i.e. $E_\textrm {rec}=\beta _\textrm {rec}B_0$. Since $m_e c^2/eB_0\equiv \rho _0$, $\Delta \gamma _\textrm {ext} \sim \beta _\textrm {rec} cT/\rho _0 \sim 0.02 (cT/L_x) L_x/\rho _0$. In our simulations, $\beta _\textrm {rec}$ decreases over time (as reconnection slows to a stop), but above we measured $\beta _\textrm {rec}\approx 0.02$ in the middle of active reconnection; and in this simulation, $L_x=1280\sigma \rho _0$, so $\Delta \gamma _\textrm {ext} \sim 26 (cT/L_x)\sigma$. Thus, after $T\sim 10L_x/c$, the maximum energy gain would be around $\Delta \gamma _\textrm {ext}\sim 250 \sigma$. However, rather than estimate $\beta _\textrm {rec}$ and the active reconnection time $T$, it is convenient to consider that by Faraday's law, $\int _0^T E_\textrm {rec}\,\textrm {d}t = -\Delta \psi /(cL_z)$ where $\Delta \psi$ is the change in unreconnected upstream flux; in a 2-D 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 _\textrm {ext} \approx (L_y/4 L_x) (\Delta \psi /\psi _0) (eB_0 L_x/m_e c^2) = 0.5(0.45)(L_x/\rho _0)$, or $\Delta \gamma _\textrm {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 _\textrm {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_\textrm {rec}$. We can infer from figure 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 _\textrm {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 of the order of $\gamma \approx 20\sigma$. However, no particle ever exceeds $\gamma \approx 40$ significantly, even after much longer times. We note that: (i) we measured $\Delta \psi$ with a time cadence of $0.44 L_x/c$, so we cannot examine extreme acceleration on smaller time scales, and (ii) 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. Reference Werner, Uzdensky, Cerutti, Nalewajko and Begelman2016), and experiencing slower acceleration thereafter (as described by Petropoulou & Sironi Reference Petropoulou and Sironi2018; Hakobyan et al. Reference Hakobyan, Petropoulou, Spitkovsky and Sironi2021).

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 vs 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 3-D simulations with different lengths $L_z$). In the following subsections, we consider 2-D 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.

4.2. Two-dimensional reconnection: dependence on initial current sheet configuration

In this subsection we show that 2-D 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 2-D 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 _\textrm {plasma}=0.5$, and $B_{gz}=0$, as well as $L_y/L_x=2$).

The initial magnetic field perturbation $a$ (cf. (2.2)), 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.

Case (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 figure 3(a) we observe that the evolution of magnetic energy vs 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 figure 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.

Figure 3.

(a) Transverse magnetic energy vs time, and (b) electron energy spectra $f(\gamma )$, compensated by $\gamma ^4$, at times $t$ when $U_{Bt}(t)=0.75U_{B0}$ – marked by the horizontal grey line in (a) – for different initial magnetic field perturbations in 2-D simulations with an initially dense, thin current sheet ($\eta =5$, $\delta /\sigma \rho _0=2/3$). Panels (c,d) show the same for a less dense, thicker layer ($\eta =1$, $\delta /\sigma \rho _0=10/3$), with spectra at $U_{Bt}(t)=0.77U_{B0}$.

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 figure 3(c,d). Here, $s>\Delta x$ as long as $a > 0.011$ (i.e. for all the non-zero $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 figure 3(c), which shows magnetic energy vs time, shifted relative to $t_\textrm {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 figure 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 (figure 4a,c), 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 (figure 4b,d); 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.

Figure 4.

Magnetic field lines superimposed over the density $n_{de}$ of the initially drifting electrons that form the initial current sheet (i.e. excluding background electrons), for cases (a,c) without an initial magnetic perturbation ($a=0$), and (b,d) with a small initial perturbation $a=5.5$ or $s/\delta =3$ (cf. figure 3a), at times $t=0.44L_x/c$ (a,b) and $t=1.2L_x/c$ (c,d). In both cases, the drifting particles are well contained in plasmoid cores. With zero perturbation (a,c), the first plasmoids formed are roughly equal in size, and they generally merge in pairs so that at any time, there is some number of equal-sized largest plasmoids, each containing a similar fraction of the drift particles. In contrast, a small perturbation (b,d) favours one plasmoid, which remains always larger than all other plasmoids, and rapidly gathers most and then all of the drifting particles. Magnetic field lines in the layer (i.e. within the separatrix) are solid green, with the separatrix being solid cyan; upstream of the separatrix, field lines are dashed blue. These simulations have size $L_x=320\sigma \rho _0$, with an initially thin current sheet, $\delta =(2/3) \sigma \rho _0$.

Summarizing case (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 2-D semirelativistic electron–ion reconnection (Ball et al. Reference Ball, Sironi and Özel2019); there NTPA was found to be similarly insensitive to triggering, at least for weak guide field.

Case (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_\textrm {ad}$ and temperature $\theta _\textrm {ad}$ sufficient to balance the upstream magnetic pressure $B_0^2/8{\rm \pi}$ as well as the upstream plasma pressure $n_{b0}\theta _b m_ec^2$, this would yield $\eta _\textrm {ad}\equiv n_\textrm {ad}/n_{b0}=2.28$ and $\theta _\textrm {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. (2.2)). 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 (a) the magnetic energy evolution and (b) 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.

Figure 5.

For simulations with varying $\eta$ (hence varying initial current sheet half-thicknesses $\delta$; cf. table 1): (a) transverse magnetic energy vs time (shifted by $t_\textrm {onset}$; cf. § 3.5), and (b) electron energy spectra times $t$ when $U_{Bt}(t)=0.75 U_{B0}$ (indicated by the horizontal grey line in the left panel). Dotted grey-on-yellow lines indicate power-law slopes of 4.1 and 4.8. All simulations have size $L_x=1280\sigma \rho _0$.

We thus conclude that the initial current sheet density has a negligible effect on total energy conversion and NTPA efficiency as long as $\eta \gg 0.1$; however, for small $\eta$, NTPA efficiency and total energy conversion decrease as $\eta$ decreases, and the thick, low-density current sheet may eventually become effectively stable against tearing and reconnection.

Case (2D-c): varying $\beta _d\propto \eta ^{-1}$; i.e. varying $\delta /\rho _d$ for fixed $\delta /\sigma \rho _0=(10/3)\sigma \rho _0$. The motivation for this particular parameter space exploration is a little complicated. Essentially, we would like to study the influence of $\beta _d$ on reconnection. After substudy (2D-b), which varied $\eta$ but not $\beta _d$, it would be natural to hold $\eta$ constant while varying $\beta _d$. However, we will instead vary $\eta \propto \beta _d^{-1}$ to keep $\eta \beta _d$ constant. Importantly, this fixes $\delta =\sigma \rho _0/(\eta \beta _d)$ with respect to upstream plasma scales. In particular, $\delta /\Delta x=10$ and $L_x/\delta =384$ will be the same in every case, so that, as $\beta _d$ and $\eta$ vary, the initial current sheet will remain well resolved but very thin compared with $L_y=2L_x$. Satisfying these two criteria will be important in three dimensions (§§ 5.3–5.4) when we will focus on the evolution of the current layer at early times as it undergoes large deformations due to RDKI, requiring $\delta /\Delta x\gg 1$ and $\delta /L_y \ll 1$; in addition, it will be more straightforward to compare layer deformation and its interaction with the upstream plasma for different simulations when the current sheets begin with the same thickness. Although neither RDKI nor extreme layer deformation will occur in two dimensions, this substudy provides a baseline against which to compare the 3-D simulations.

Therefore, we simultaneously vary $\eta$ along with $\beta _d$ to keep the product $\eta \beta _d$ constant. Maintaining constant $\eta \beta _d$ fixes $\delta /\sigma \rho _0$, while essentially varying the drifting plasma length scales ( such as gyroradius $\rho _d=3\theta _d\rho _0$) with respect to $\delta$ – i.e. varying $\delta /\rho _d = 2/(\gamma _d \beta _d)$.

We compare five configurations with the same background plasma (with zero guide field), system size $L_x=1280\sigma \rho _0=384\delta$ and zero initial magnetic perturbation ($a=0$). We vary $\beta _d\in \{$0.075, 0.15, 0.3, 0.6, 0.9$\}$; to keep $\eta \beta _d$ constant, the overdensities $\eta$ are $\{$4, 2, 1, 0.5, 0.33$\}$, yielding the same $\delta =(10/3)\sigma \rho _0=10\Delta x$ for all five cases. Other parameters are (see also table 1): $\theta _d/\theta _b \in \{0.5$, 1.0, 2.1, 5.0, 14$\}$ and $\delta /\rho _d \in \{8.9$, 4.4, 2.1, 0.89, 0.32$\}$, respectively.

At $\beta _d=0.075$, $\eta =4$, the gyroradius and collisionless skin depth of the drifting plasma are at their smallest values of $\rho _d=0.38\sigma \rho _0$ and $d_{ed}=0.31\sigma \rho _0$, both near the grid resolution $\Delta x$, and the half-thickness $\delta$ is 10 times larger than those scales. At the other end of the parameter scan, $\beta _d=0.9$, $\eta =0.33$, the microphysical scales of the initial current sheet become somewhat larger than $\delta$: $\delta \approx 0.3 \rho _d \approx 0.6 d_{ed}$.

To characterize the effects of varying $\beta _d$ (with $\eta$), we focus on the time history of magnetic and plasma energy, and on the resulting NTPA; and again, we find little variation among this set of simulations (see figure 6). All cases show very similar magnetic depletion rates and NTPA (within stochastic variation), except for the case $\beta _d=0.9$ ($\eta =0.33$), which is a little different; this case ends up with a bit more energy in plasmoids, although the rate at which upstream magnetic energy and unreconnected flux decrease is similar to the others. The $\beta _d=0.9$ case has significantly hotter (although lower density) drift particles (cf. table 1), which results in correspondingly hotter and less dense plasmoids cores, since the initially drifting particles, being close to the initial midplane, are among the first to be swept into plasmoids and remain near the plasmoid cores (the original drift kinetic energy may also be converted to thermal energy, contributing to even hotter cores).

Figure 6.

(a) Transverse magnetic energy vs (shifted) time, and (b) electron energy spectra when 23 per cent of the initial magnetic energy has been converted (see the horizontal grey line in panel a), for simulations with varying initial current drift $\beta _d$, but constant initial sheet thickness $\delta =(10/3)\sigma \rho _0$. (Table 1 shows how initial temperatures and gyroradii vary with $\beta _d$.) Dotted grey-on-yellow lines indicate power-law slopes of 4.2 and 4.6. All simulations have size $L_x=1280\sigma \rho _0$.

Summary of 2D-a,b,c. We have seen that in two dimensions, the details of the initial current sheet can affect the time delay before the onset of reconnection, as well as the evolution of the plasmoid hierarchy. Nevertheless, the initial current sheet configuration has relatively little effect on the long-term magnetic energy conversion and NTPA in 2-D reconnection; after a short initial stage, reconnection energetics are dominated by the upstream plasma, regardless of the very early current sheet evolution. This conclusion may seem to be at odds with that from Ball et al. (Reference Ball, Sironi and Özel2019), which found subtle but significant differences in NTPA due to initial perturbation and varying initial sheet thickness; however, Ball et al. (Reference Ball, Sironi and Özel2019) studied transrelativistic electron–proton reconnection, with low upstream plasma beta and non-zero guide magnetic field ($B_{gz}/B_0=$0.1 and 0.3). They used a small pressure imbalance as an initial perturbation, and found that its influence diminished with diminishing guide field, potentially consistent with our zero-guide-field results here. The impact of varying sheet thickness (keeping $\eta$ constant) was also seen at the larger guide field, $B_{gz}=0.3B_0$, so there are several potentially important differences between these two studies: transrelativistic electron–proton vs ultrarelativistic pairs, low upstream plasma beta vs beta of 0.5 ($\sigma _h=1$), weak-to-moderate vs zero guide field, and pressure vs magnetic field perturbation.

4.3. Two-dimensional reconnection: system-size dependence

In this subsection we investigate how reconnection depends on system size $L_x$, to estimate the applicability of our kinetic simulations to astrophysical systems, which may be much larger than even the largest supercomputers can possibly simulate while including full kinetics. We present results for system sizes ranging over a factor of 30, from $L_x/\sigma \rho _0=80$ to 2560. These simulations achieve a substantial scale separation between the system size $L_x$ and all the microphysical scales, which are of order $\sigma \rho _0$ (cf. § 2). Such large scale separations are feasible because, for $\sigma _h=1$, all the kinetic scales are comparable, hence they can all be resolved by a cell size $\Delta x$ only marginally smaller than $\sigma \rho _0$. All simulations presented in this subsection are two-dimensional with zero guide field, starting with an initial current sheet half-thickness $\delta =(2/3)\sigma \rho _0$, $\eta =5$, $\beta _d=0.3$ and zero perturbation.

Before discussing how $L_x$ affects typical reconnection behaviour, we note that stochastic evolution of the plasmoid hierarchy affects energy conversion and NTPA, especially over short time scales. For example, in the $L_x=1280\sigma \rho _0$ simulation, the magnetic energy depletion slows (compared with earlier times) between $t\approx 3L_x/c$ and $t \approx 5L_x/c$, so that $U_{Bt}(t)$ is distinctly higher at $t=5L_x/c$ than in all the other simulations (figure 7a, inset). At $t\approx 3L_x/c$, the lower layer in the simulation domain for $L_x=1280\sigma \rho _0$ has two large, roughly equal-sized plasmoids, separated by $\approx L_x/2$ (while the upper layer already has only one final plasmoid at this time). In principle, one of the plasmoids could either move left to merge with the other, or (since the simulation is periodic in $x$) it could move right. Thus they reach a ‘Buridan's Ass’ metastable state that seems to slow reconnection evolution for a short time. However, by $4.4L_x/c$, the two plasmoids have decided which way to move, and by $4.8L_x/c$ they have started to merge, causing a burst of magnetic energy depletion; by the time the merger is nearly complete at $6.1L_x/c$, $U_{Bt}(t)$ has ‘caught up’ to all the other simulations. Over short times, the reconnection dynamics thus depends on the detailed behaviour of large plasmoids, which has an element of randomness; however, the random plasmoid behaviour does not alter the long-term behaviour of reconnection.

Figure 7.

For different system sizes $L_x$: (a) transverse magnetic energy vs time, with an inset showing $L_x=1280\sigma \rho _0$ alone, and (b) unreconnected flux vs time. The horizontal dashed lines indicate the periods over which the average reconnection rate is given in figure 8, and slopes corresponding to $\beta _\textrm {rec}=0.04$ and 0.02 (cf. § 3.3) are shown as grey dotted lines.

Figure 7 shows the evolution of transverse magnetic energy $U_{Bt}(t)$ and unreconnected flux $\psi (t)$ for a range of $L_x$. Although $U_{Bt}(t)$ often decreases fitfully, $\psi (t)$ tends to decrease more smoothly. Initially, $\psi (t)$ falls very rapidly as reconnection is dominated by the initial current sheet, and the rate of fall slows a little as the upstream plasma begins to dominate reconnection (the larger $L_x/\delta$ is, the sooner this happens relative to the total duration of reconnection). Eventually, the system reaches a stable magnetic configuration and $\psi (t)$ approaches a constant value; therefore, the reconnection rate must decrease over time, eventually toward zero (in a closed system). The reconnection rates over two different time intervals are shown in figure 8(a), namely for the two intervals $0.9 > \psi (t)/\psi _0 > 0.8$ and $0.8 > \psi (t)/\psi _0 > 0.7$, i.e. between the horizontal dashed lines in figure 7(b). As the system size $L_x/\sigma \rho _0$ increases from 160 to 2560, $\beta _\textrm {rec}$ (cf. § 3.3) falls from $\approx 0.04$ to $\approx 0.03$ in the earlier interval, and from $\beta _\textrm {rec}\approx 0.03$ to $\beta _\textrm {rec}\approx 0.02$ in the later interval. However, $\beta _\textrm {rec}$ stabilizes within the uncertainty of measurement by $L_x\gtrsim 640\sigma \rho _0$.

Figure 8.

(a) The average reconnection rate (normalized to $B_0 c$ on the left scale, and to $B_0 v_A$ on the right scale) vs system size $L_x$, over the period when the unreconnected flux falls from 0.9–0.8$\psi _0$ (blue squares), and also from 0.8–0.7$\psi _0$ (red circles) – i.e. between the horizontal dashed lines in figure 7(b). (b) The fraction of initial (transverse) magnetic field energy converted to plasma energy (blue diamonds) and of upstream flux depleted (red plusses) during reconnection vs system size $L_x$ (for aspect ratio $L_y/L_x=2$, with two layers). The largest system size, $L_x=1280\sigma \rho _0$ had not quite finished reconnection when the simulation halted around $10L_x/c$; it would have depleted slightly more energy and flux if given more time.

Reconnection continues until approximately $0.27U_{B0}$ and $0.42\psi _0$ have been depleted; the fractions of lost energy and flux are practically independent of $L_x$ (figure 8b), although we note that they are specific to the fixed aspect ratio $L_y/L_x=2$. The largest simulation, $L_x=2560\sigma \rho _0$ shows slightly less loss in $U_{Bt}$ and $\psi$, likely because it had not quite finished reconnecting when it halted around $10L_x/c$ (cf. figure 7).

Importantly, the upstream magnetic energy and flux are both depleted by the same amount – 42 per cent. However, of that $0.42U_{B0}$ initially associated with magnetic field lines that eventually undergo reconnection, only $0.27U_{B0}$ is converted to plasma energy; the remaining $0.15 U_{B0}$ ends up stored in the ‘reconnected’ magnetic field around and within plasmoids.

Turning our attention to NTPA, figure 9 shows the electron energy spectra $f(\gamma )$ at $t=6L_x/c$ for simulations with different $L_x$ (panel a shows $\gamma ^2 f(\gamma )$, and b shows $\gamma ^4 f(\gamma )$ to enhance small differences). Aside from the smallest, all simulations have a clear high-energy non-thermal power law $f(\gamma )\sim \gamma ^{-p}$; for the largest simulation, it extends over almost a decade in energy. At $t=6L_x/c$, the power-law index $p$ steepens slightly as $L_x/\sigma \rho _0$ increases from 160 to 2560, from $p\approx 4.0$ to $p\approx 4.5$ (figure 9c—cf. § 3.7). This steepening may be caused by a different acceleration mechanism that operates over long times (hence only in very large simulations), as we will discuss presently.

Figure 9.

For different system sizes $L_x$: (a) electron energy spectra $f(\gamma )$, compensated by $\gamma ^2$, at time $t=6L_x/c$, with power-law indices of 4 and 4.4 indicated by dotted grey-on-yellow lines, (b) $\gamma ^{4} f(\gamma )$ at $t=6L_x/c$, to enhance differences, (c) the median power-law index during the time period $5 < tc/L_x < 7$, with ‘error’ bars containing 68 per cent of the values fitted to spectra within this period. For $L_x/\sigma \rho _0=$256 and 320, results from two simulations are shown (slightly offset in $L_x$ for clarity in panel c).

The time evolution of power-law index $p$ and high-energy cutoff $\gamma _c$ (cf. § 3.7), for all the $L_x$, are shown in figure 10(a,b). The initial development of the power law, as the first particles reach high energies, results in $p(t)$ rapidly falling (i.e. the power law hardens/flattens) to approach an index between 3.9 and 4.5. From its minimum value, the slope appears to steepen slowly over time, mostly notably for simulations with $L_x\geq 640\sigma \rho _0$, all of which show $p(t)$ rising by approximately 0.2 from its minimum.

Figure 10.

For different system sizes $L_x$, the fitted values for the power-law index $p$ (panel a) and high-energy cutoff $\gamma _c$ (b) vs time. (c) For $L_x=2560\sigma \rho _0$, $p(t)$ (blue, left axis) and $\gamma _c(t)$ (red, right axis) are essentially the same for two different resolutions, $\Delta x/\sigma \rho _0 = 1/3$ and $1/4$, even though for this large system size, $\Delta x/\sigma \rho _0=1/3$ resulted in unphysical total energy growth of 1 per cent, while $\Delta x / \sigma \rho_0 = 1/4$ saw a more satisfactory 0.05 per cent increase. Error bars indicate uncertainty, encompassing the middle 68 per cent of values measured over the associated time interval.

Correspondingly, the cutoff $\gamma _c(t)$ grows rapidly as the power law develops, attaining after only 1–2$L_x/c$ a value ranging from $\gamma _c=6\sigma$ for $L_x=80\sigma \rho _0$ to $\gamma _c=15\sigma$ for $L_x=2560\sigma \rho _0$. For $L_x\geq 1280\sigma \rho _0$, $\gamma _c(t)$ continues to grow over time, and even the case $L_x=640\sigma \rho _0$ experiences some limited growth after approximately $t=7L_x/c$. Thus at late time $t=9L_x/c$, $\gamma _c$ ranges from $6\sigma$ to $33\sigma$ as $L_x$ increases by a factor of 32, consistent with $\gamma _c \sim \sqrt {L_x}$; importantly, although $\gamma _c$ increases with $L_x$, the increase is clearly sublinear. This is consistent with very rapid (extreme) acceleration up to $\gamma _c\sim 4\sigma$ (Werner et al. Reference Werner, Uzdensky, Cerutti, Nalewajko and Begelman2016; Uzdensky Reference Uzdensky2020), at which energy particles become trapped in plasmoids and cease extreme acceleration, but undergo much slower acceleration inside plasmoids (Petropoulou & Sironi Reference Petropoulou and Sironi2018; Uzdensky Reference Uzdensky2020; Hakobyan et al. Reference Hakobyan, Petropoulou, Spitkovsky and Sironi2021). Interestingly, the slow growth in $\gamma _c(t)$ seems to correspond to the slow steepening in $p(t)$, suggesting that the slower acceleration mechanism may yield a steeper power law. The increase in $p$ with $L_x$ seen in figure 9(c) may simply reflect the increasing influence of the slower acceleration mechanism after longer times (accessible in large systems).

Very large reconnection simulations sometimes suffer at late times from slow-growing instabilities. The largest case, $L_x=2560\sigma \rho _0$, exhibited potentially troubling energy non- conservation, with the total energy growing by almost 1 per cent between $t=7L_x/c$ and $t=10L_x/c$. We find that this typically results in unphysical heating of the lowest-energy particles but does not affect the high-energy spectrum. To rule out this possibility, we compare the simulation with fiducial resolution $\Delta x=\sigma \rho _0/3$ to one with higher resolution, $\Delta x=\sigma \rho _0/4$ in figure 10(c); both show essentially similar evolution of $p(t)$ and $\gamma _c(t)$, although the higher-resolution case suffered only a 0.05 per cent increase in total energy.

In summary, the system-size dependence of 2-D reconnection with $\sigma _h=1$ is weak but not entirely insignificant for $L_x \gtrsim 160\sigma \rho _0$. The reconnection rate (normalized to $cB_0$) is $\beta _\textrm {rec}\approx 0.03$ or $(c/v_A)\beta _\textrm {rec}\approx 0.04$. During reconnection, some particles are very rapidly accelerated into a power-law energy spectrum $\sim \gamma ^{-p}$ with $p\approx 3.9$–$4.4$, extending up to a cutoff around $\gamma _c \gtrsim 6\sigma$. This rapid acceleration presumably ends when particles become trapped in plasmoids. Over time, a much slower acceleration mechanism continues to accelerate particles trapped in plasmoids (Petropoulou & Sironi Reference Petropoulou and Sironi2018; Hakobyan et al. Reference Hakobyan, Petropoulou, Spitkovsky and Sironi2021), resulting in a slight steepening of $p$ (up to $\approx 4.6$ in our largest simulation) and slow but significant growth in $\gamma _c$ (up to $\approx 33\sigma$ in our largest simulation).

4.4. Two-dimensional reconnection: guide magnetic field

We now investigate the effect of adding an initial uniform guide magnetic field $B_{gz}\hat {\boldsymbol {z}}$, considering $B_{gz}/B_0 \in \{$0, 0.25, 0.5, 0.75, 1, 1.5, 2, 3, 4$\}$. All simulations presented in this subsection are two-dimensional with $L_x=1280\sigma \rho _0$, starting with an initial current sheet half-thickness $\delta =(2/3)\sigma \rho _0$ ($\eta =5$, $\beta _d=0.3$), and zero perturbation ($a=0$).

A weak guide field $B_{gz} \lesssim 0.25 B_0$ has very little effect on the magnetic energy evolution. A strong guide field, however, slows reconnection and inhibits overall magnetic energy conversion, as shown in figure 11. Figure 11(a) shows the transverse magnetic energy $U_{Bt}$ vs time, for simulations with $B_{gz}/B_0$ ranging from 0 to 4; as $B_{gz}$ increases above $\simeq 0.5B_0$, $U_{Bt}(t)$ falls more slowly, and ultimately decreases by a smaller amount. Correspondingly (cf. figure 11b), the unreconnected flux decreases more slowly, and less flux is reconnected overall.

Figure 11.

For several guide-field strengths $B_{gz}$: (a) transverse magnetic energy vs time, and (b) unreconnected flux vs time, with dot-dashed lines showing the slopes corresponding to $\beta _\textrm {rec}=$0.04 and 0.008 (cf. § 3.3). To estimate stochastic variability, two simulations are shown for $B_{gz}=0$, two for $B_{gz}=0.25B_0$ and three for $B_{gz}=0.5B_0$.

This slowing of reconnection has been previously attributed to the decrease in the effective component ($x$-projection) of the Alfvén velocity (Liu et al. Reference Liu, Daughton, Karimabadi, Li and Peter Gary2014, Reference Liu, Guo, Daughton, Li and Hesse2015; Werner & Uzdensky Reference Werner and Uzdensky2017). The upstream Alfvén velocity, including $B_0$ and $B_{gz}$, is $v_A=c\sqrt {(B_0^2+B_{gz}^2)/(B_0^2+B_{gz}^2 + 4{\rm \pi} h)}$, where $h$ is the plasma enthalpy density ($h=4n_b\theta _b m_e c^2$ in the ultrarelativistic limit). The projection of $v_A/c$ along $B_0$ is therefore

(4.1)\begin{equation} \frac{v_{A,x}}{c} = \sqrt{ \frac{B_0^2}{B_0^2 + B_{gz}^2}} \frac{v_A}{c} = \sqrt{ \frac{B_0^2}{B_0^2 + B_{gz}^2 + 4{\rm \pi} h}} = \sqrt{ \frac{B_0^2}{B_0^2 + 4{\rm \pi} h_\textrm{eff}}}, \end{equation}

where, following Werner & Uzdensky (Reference Werner and Uzdensky2017), we define $h_\textrm {eff}\equiv h+B_{gz}^2/4{\rm \pi}$ to be an effective enthalpy density including a contribution from the guide field. Then, defining $\sigma _{h,{\rm eff}}\equiv B_0^2/(4{\rm \pi} h_\textrm {eff})$, we have $v_{A,x} = (\sigma _{h,{\rm eff}}^{-1}+1)^{-1/2}$, which reduces to $v_{A}/c = (\sigma _h^{-1}+1)^{-1/2}$ for $B_{gz}=0$.

Figure 12 quantifies the reconnection rates (panel a) and amount of magnetic energy and upstream flux depleted (panel b). The left panel shows the reconnection rates, averaged over the time interval during which $\psi (t)$ falls from $0.9\psi _0$ to $0.8\psi _0$. The reconnection rate, normalized to $B_0 c$ (blue circles) falls from roughly 0.03–0.04 for $B_{gz}\lesssim 0.75B_0$ to around 0.01 for $B_{gz} = 4 B_0$. The red squares show reconnection rates normalized to $B_0 v_{A,x}$. For $B_{gz}/B_0\gtrsim 1$, $(c/v_{A,x})\beta _\textrm {rec}$ is relatively constant around 0.03; however, $(c/v_{A,x})\beta _\textrm {rec}$ is closer to 0.04 or even 0.05 for $B_{gz}/B_0 < 1$. As a compromise, figure 12 shows a dashed line at $(c/v_{A,x})\beta _\textrm {rec}=0.038$; the dotted line shows the corresponding $\beta _\textrm {rec}$.

Figure 12.

(a) Reconnection rate vs $B_{gz}$: blue circles show the reconnection rate $\beta _\textrm {rec}$ (normalized to $B_0 c$), and red squares show $(c/v_{A,x})\beta _\textrm {rec}$ – cf. (4.1). The blue dotted line shows what $\beta _\textrm {rec}$ would be if $(c/v_{A,x})\beta _\textrm {rec}$ were constant at 0.038 (the red dashed line). The reconnection rates are averages over the period during which $\psi (t)$ falls from $0.9\psi _0$ to $0.8\psi _0$ (i.e. between the horizontal grey, dashed lines in figure 11b). (b) The fractional amount of transverse magnetic energy $U_{Bt}$ (blue x symbols) and upstream magnetic flux (red plusses) lost during reconnection, vs $B_{gz}$, along with empirical fits for both.

According to (4.1), the guide field should not significantly suppress reconnection for $B_{gz}/B_0 \ll \sigma _h^{-1/2} = 1$. We note that a much stronger guide field would be needed to suppress non-relativistic reconnection (where $\sigma _h \ll 1$). On the other hand, our previous work (Werner & Uzdensky Reference Werner and Uzdensky2017) considered $\sigma _h \simeq 25$, and correspondingly observed that even $B_{gz} = B_0/2$ significantly suppressed reconnection, slowing the reconnection rate and reducing the overall magnetic depletion (although $B_{gz}=B_0/4$ had only a weak effect).

Figure 12(b) shows that the total fraction of $U_{Bt}$ converted to $U_\textrm {plasma}$ decreases with stronger guide field – from around 27 per cent for $B_{gz}=0$ to only approximately 6 per cent for $B_{gz}=4 B_0$; similarly, the amount of flux reconnected decreases from 42 per cent to around 21 per cent. We find that the fractional loss in $\psi$ is fit fairly well by $0.86(L_x/L_y)(1+2B_{gz}^2/B_0^2)^{-1/5}$, and the loss in $U_{Bt}$ by $0.56(L_x/L_y)(1+4B_{gz}^2/B_0^2)^{-2/5}$; these functional forms are meant for convenient comparison with future results and we offer no justification for them.

Even a weak guide appears to enhance the reconnection rate slightly. However, because stochastic variability can yield different reconnection rates, more work is needed to determine whether this trend is statistically significant.

With stronger guide field, reconnection ends (or reaches a stable configuration) with a smaller amount of magnetic energy and flux depletion. This may be caused by the increased pressure of the guide magnetic field in plasmoids. The guide magnetic field pressure resists compression with adiabatic index $\varGamma =2$, compared with only $\varGamma =4/3$ for relativistic plasma alone (e.g. compressing a flux tube with uniform magnetic field from volume $V_i$ to $V_f$ while preserving the magnetic flux inside requires work $\propto (V_i/V_f)^{\varGamma -1}-1$ where $\varGamma =2$; whereas compressing the plasma in the tube requires work $\propto (V_i/V_f)^{4/3-1}-1$). As a result, for the same amount of (transverse magnetic field) flux reconnected, the major plasmoid will be larger when the guide field is stronger. Therefore, the growth of the major plasmoid reaches the system-size scale earlier (in terms of amount of flux reconnected, not time) and shuts down reconnection.

Perhaps not surprisingly, the slowing of reconnection inhibits NTPA, as shown in figure 13, which displays the final particle energy spectra (panel a) and the fitted power-law indices (panel b). Overall, less energy is transferred from the magnetic field to the plasma, and at any given well-above-average energy, stronger guide field results in a smaller fraction of particles at that energy, hence NTPA is less efficient. As $B_{gz}$ increases, the power law steepens from $p\simeq 4$ to $p\gtrsim 10$ at $B_{gz}=4B_0$. The line $p\approx 4.4 + 2B_{gz}/B_0$ very roughly captures this trend, although it is important to remember that the uncertainty in measuring steep power laws ($p\gtrsim 4$) can be quite high (our fitting method returns consistent values characterizing the spectra, but fundamentally one should question what it means to measure a steep power law, e.g. $\sim \gamma ^{-8}$ as $\gamma$ varies over just a half-decade). This general trend was also observed for $\sigma _h=25$ reconnection in Werner & Uzdensky (Reference Werner and Uzdensky2017), where the steepening was fairly well described by $p=1.9+0.7\sigma _{h,{\rm eff}}^{-1/2}$; however, this formula provides a very poor fit for the $\sigma _h=1$ simulations, perhaps because we are in the regime of very steep power laws that cannot be reliably measured at higher energies without orders of magnitude more simulated particles.

Figure 13.

(a) Particle energy spectra $f(\gamma )$ at $t=12L_x/c$, for different imposed guide magnetic fields, show that increasing guide field inhibits NTPA, resulting in steeper power laws. Slopes with power-law indices of 4.4, 4.8 and 7 are shown for comparison. (b) The steepening is quantified by power-law indices fit to the high-energy part of $f(\gamma )$, with error bars encompassing the middle 68 per cent of fitted indices over times from 4 to 10$L_x/c$. The line $4.4+2B_{gz}/B_0$ is shown as a guide, but it is important to remember that precise measurement of steep power laws is very difficult, and in fact the cases with $B_{gz}/B_0 \geq 1.5$ are not very different from the initial Maxwellian at $t=0$ (shown in dashed grey).

5.1. Overview of 3-D reconnection with $\sigma _h=1$, and comparison with two dimensions

In this section we begin with an overview of basic 3-D current sheet behaviour for one specific configuration, namely: $\sigma _h=1$, zero guide field ($B_{gz}=0$), zero initial perturbation ($a=0$), an initially thin current sheet $\delta =(2/3)\sigma \rho _0$ ($\eta =5$) and $L_x=512\sigma \rho _0 = L_y/2$ (the largest $L_x$ of any 3-D simulations in this paper). We will show results from simulations with $L_z/L_x=0$ (a 2-D simulation), 1/32, 1/8 and 1; we will show in particular that $L_z=L_x/32$ behaves essentially like the 2-D simulation, while the two largest simulations are similar to each other but quite different from the smaller two, presumably because of 3-D effects that appear for sufficiently large $L_z$.

Before looking at the evolution of 2-D and 3-D current sheets, it is helpful to see how global energies and unreconnected flux vary in time. We begin by looking at global field energies vs time, normalized to initial magnetic energy $U_{B0}$, in figure 14(a–c). In these simulations without guide field, almost all the field energy is contained in $U_{Bt}$, the energy in the transverse magnetic field $\boldsymbol {B}_t = B_x \hat {\boldsymbol {x}} + B_y \hat {\boldsymbol {y}}$. In all cases $U_{Bz}(t) < 0.003U_{B0}$, and the energy in the electric field $\boldsymbol {E}$ (not shown) also remains below $0.003 U_{B0}$. The $B_z$ and $E$ field components can nevertheless be substantial within the (the small volume of) the layer. Further investigation of $B_z$ is left for future work, but we note that there should be no $B_z$ at all in a 2-D magnetohydrodynamic description (for pair plasma with no initial guide field, Zenitani, Hesse & Klimas Reference Zenitani, Hesse and Klimas2009), while 2-D PIC simulations do develop non-zero $B_z$ (Zenitani & Hesse Reference Zenitani and Hesse2008).

Figure 14.

Energy in various magnetic field components vs time for 4 different $L_z/L_x$, normalized to $U_{B0} \approx \int \,\textrm {d}V \, B_0^2/8{\rm \pi}$: (a) $U_{Bt} = \int \,\textrm {d}V \, (B_x^2+B_y^2)/8{\rm \pi}$, (b) $U_{Bz} = \int \,\textrm {d}V \, B_z^2/8{\rm \pi}$ and (c) $U_{Bt,{\rm layer}}$, the magnetic energy in the layer (cf. § 3.1). The electric field energy (not shown) is always less than $0.003 U_{B0}$. Finally (d), the unreconnected magnetic flux $\psi (t)$; dotted grey lines indicate dimensionless reconnection rates, $(cB_0L_z)^{-1} (\textrm {d}\psi /\textrm {d}t)$, of 0.045 and 0.015, which are roughly the rates of nearly 2-D and 3-D simulations, respectively, between the grey dashed lines, and also 0.002, roughly the 3-D rate at later times.

Figure 14(a) shows $U_{Bt}$ decreasing over time in all cases as it is converted to plasma energy. The simulations with larger $L_z/L_x=1/8,1$ show similar magnetic depletion, and those with smaller $L_z/L_x=0,1/32$ closely resemble each other; in this and other diagnostics we will see that the former ($L_z/L_x=1/8,1$) exhibit 3-D behaviour, while the latter are essentially two-dimensional. We find that magnetic energy depletes roughly half as fast in three as in two dimensions, at early times up to $\sim 4L_x/c$ (figure 14a; see also table 2, which compares energy changes and reconnection rates). While the 2-D simulations have mostly finished reconnection by $t\sim 4L_x/c$ (approaching a final state a few $L_x/c$ later), the 3-D simulations embark upon a second stage of prolonged slow energy conversion (cf. Liu et al. Reference Liu, Li, Yin, Albright, Bowers and Liang2011), converting magnetic energy to plasma energy an order of magnitude still more slowly for many tens of $L_x/c$. Were these simulations to run many times longer (which they did not, because of prohibitive computational cost), the 3-D simulations would eventually convert more magnetic energy than the essentially 2-D simulations; we will see in § 5.5, using smaller simulations ($L_x=341\sigma \rho _0$), that 3-D simulations convert more magnetic energy over $\sim 50 L_x/c$ than 2-D simulations.

Table 2.

The fractional loss in transverse magnetic field energy $U_{Bt}$ in the entire simulation (nearly equal to the gain in plasma energy $U_\textrm {plasma}$), divided into upstream and layer regions ($U_{Bt,{\rm up}}+U_{Bt,{\rm layer}}=U_{Bt}$), and the loss in unreconnected magnetic flux $\psi$, over the first $10L_x/c$, for simulations with $L_x=512\sigma \rho _0$ and $B_{gz}=0$. Changes are shown normalized to initial magnetic energy $U_{B0}$ and initial unreconnected flux $\psi _0$. The fully 3-D simulations would continue to deplete magnetic energy if run longer (cf. § 5.5). Also shown is $\beta _\textrm {rec}$ (cf. § 3.3), averaged over the time during which $\psi (t)$ drops from 0.95$\psi _0$ to 0.90$\psi _0$.

More magnetic energy conversion is possible in three dimensions because the magnetic energy pumped into plasmoids can be depleted in three dimensions, whereas in 2-D plasmoids remain stable. To demonstrate how little magnetic energy remains stored in plasmoids in three dimensions, we separately measure the transverse magnetic energy in the upstream and in the layer regions (cf. § 3.1): $U_{Bt} = U_{Bt,{\rm up}} + U_{Bt,{\rm layer}}$. With no guide field and no initial perturbation, the initial magnetic energy $U_{B0}=U_{Bt0}$ is entirely in unreconnected transverse field. During magnetic reconnection in two dimensions, part of this gets converted into plasma energy $U_\textrm {plasma}$, and part into $U_{Bt,{\rm layer}}$ in plasmoids. Over the first $10L_x/c$, the (nearly) 2-D simulations take the initial upstream magnetic energy $U_{B0}$ and convert about $\Delta U_\textrm {plasma}\approx 0.27U_{B0}$ to plasma energy while building up a comparable amount, $\Delta U_{Bt,{\rm layer}}\approx 0.20_{B0}$, in plasmoids, where it remains trapped. Since the energized particles are also trapped in plasmoids, we estimate that the magnetic energy in plasmoids is approximately 2/3 of the plasma energy in plasmoids. Over the same time, fully 3-D simulations convert $\Delta U_\textrm {plasma}\approx 0.2 U_{B0}$ to plasma energy, but magnetic energy in the layer quickly decays, leaving a strikingly smaller upper bound on the amount of energy, $\Delta U_{Bt,{\rm layer}}\approx 0.01 U_{B0}$, that could be in plasmoids (see figure 14c and table 2). Although 3-D simulations deplete upstream energy $U_{Bt,{\rm up}}$ more slowly, they thus convert it more completely to plasma energy.

Along with magnetic energy, we can also measure the unreconnected magnetic flux, which decreases over time, faster in two than in three dimensions. In three dimensions, as in two dimensions, some of the upstream flux is promptly reconnected, ending up in plasmoids, as in two dimensions. Unlike in two dimensions, however, in three dimensions some upstream flux may be promptly annihilated, and also (only in three dimensions) reconnected flux is subsequently annihilated as plasmoids or flux-rope structures become unstable and decay. It is non-trivial to measure how much upstream flux is annihilated promptly without first undergoing reconnection, but we can conclude from the relatively small amount of magnetic energy in the layer that ultimately, in three dimensions, most of the lost upstream flux is annihilated (vs being stored permanently in plasmoids). Figure 14(d) shows the unreconnected flux $\psi (t)$ diminishing as reconnection occurs for the four different values of $L_z/L_x$. Table 2 shows that up to $t=10L_x/c$, the nearly 2-D simulations lose 40 per cent of the initial unreconnected flux, and fully 3-D simulations lose approximately only 15 per cent (though, over time, they would reconnect and/or annihilate more upstream flux, whereas the 2-D simulations have essentially finished). The reconnection rate $\beta _\textrm {rec}$ at early times, measured between $\psi (t)=0.95\psi _0$ and $\psi (t)=0.90\psi _0$ (the horizontal dashed lines in figure 14d) and normalized to $B_0 c$, is approximately 0.045 for the nearly 2-D simulations, and 0.015 for three dimensions (see table 2); later, around $t=8L_x/c$, $\beta _\textrm {rec}$ is of the order of $\sim 0.002$ for the 3-D simulations (and fluctuating around zero for two dimensions).

Having shown significant differences (as well as qualitative similarities) between 2-D and 3-D energy and upstream flux depletion, we will describe the evolution of the current sheet. First, we will describe the kinking of the sheet, which cannot occur in two dimensions at all; then we will spend some time describing 3-D behaviour that resembles 2-D reconnection.

Perhaps the most expected (by now) new effect in three dimensions is the kinking of the current sheet in the $z$ direction due to the RDKI (Zenitani & Hoshino Reference Zenitani and Hoshino2007, Reference Zenitani and Hoshino2008). Simulating the $z$ dimension enables the RDKI, which can occur because of the relative movement of electrons and positrons in the $z$-direction and causes rippling of the sheet. Figure 15 shows (for the case with $L_z=L_x$) a current sheet rippling due to RDKI with an amplitude $\Delta y \approx \sigma \rho _0$ at $t=0.3L_x/c$, growing to $3\sigma \rho _0$ at $t=0.6L_x/c$ (for comparison, the initial current sheet thickness was $\delta =(2/3)\sigma \rho _0$). It had originally been thought that RDKI – at least in cases with weak guide field – might out-compete the tearing instability, thereby suppressing reconnection (Zenitani & Hoshino Reference Zenitani and Hoshino2008). However, later, larger 3-D simulations (in the high-$\sigma _h$ regime) found that, despite active RDKI, reconnection ultimately does not seem to be inhibited, even with zero guide field (Guo et al. Reference Guo, Li, Daughton and Liu2014, Reference Guo, Liu, Daughton and Li2015; Sironi & Spitkovsky Reference Sironi and Spitkovsky2014; Werner & Uzdensky Reference Werner and Uzdensky2017).

Figure 15.

The displacement $\Delta y_c(x,z,t)$ of the current sheet central surface from the original midplane (cf. § 3.2), for $L_z=L_x$ at two different times, $t=0.3L_x/c$ (a) and $0.6 L_x/c$ (b), shows increasing rippling due to the RDKI.

Indeed, 3-D simulations do exhibit behaviour that, at least locally, shares many characteristics with 2-D reconnection. Although energy is clearly transferred from magnetic fields to plasma, determining whether and where reconnection occurs is non-trivial in three dimensions. We tentatively suggest, based on our 3-D simulations, that reconnection occurs in patches within the layer – patches that may be the analogues of elementary current sheets in 2-D reconnection (i.e. the smallest inter-plasmoid current sheets, which are not further broken up by secondary tearing). It is possible that reconnection is slower (overall) in three dimensions because the area of actively reconnecting regions (patches) is smaller, while the local reconnection rate remains as high as in two dimensions (such patchiness with locally high reconnection rates has been observed for $\sigma _h\sim 1$ by Yin et al. Reference Yin, Daughton, Karimabadi, Albright, Bowers and Margulies2008). Our results are consistent with this hypothesis, but so far we have been unable to measure local reconnection rates with enough certainty to conclude this definitively; we leave a dedicated investigation of this problem to future studies. However, we offer evidence that reconnection does in fact occur (though in patches), and that the patches evolve over time, and continue to exhibit signatures of reconnection even at late times during the ‘slow 3-D reconnection stage’ (e.g. $t>4L_x/c$). Specifically, we will see that actively reconnecting regions are areas where the current sheet is thin, that they exhibit outflows (roughly) in the $\pm x$ directions, and that they are the areas with strongest parallel (to $\boldsymbol {B}$) electric fields and $\boldsymbol {E}\boldsymbol {\cdot }\boldsymbol {J}>0$.

Figure 16 shows the electron density $n_e$, normalized to $n_{be0}=n_{b0}/2$, at five different times ($tc/L_x=0.3$, 0.6, 1.1, 2.3 and 4.6) in $x$–$y$ slices, at $z=0$ for the nearly 2-D simulation with $L_z=L_x/32$, and at $z=100\sigma \rho _0$ and $z=300\sigma \rho _0$ for $L_z=L_x$. A sample of magnetic field lines (traced while ignoring $B_z$) shows that $n_{e} \lesssim n_{be0}$ in the upstream region, where flux is unreconnected (cf. § 3.3), while $n_e > n_{be0}$ in most of the layer. Therefore, looking at density is a good if rough way to track the gross evolution of the shape of the layer.

Figure 16.

A time sequence of electron density $n_e$ in the $x$–$y$ reconnection plane, in nearly two dimensions ($L_z=L_x/32$, a,d,g,j,m) and three dimensions ($L_z=L_x$, b,e,h,k,n and c,f,i,l,o) simulations for $tc/L_x=0.3, 0.6, 1.1, 2.3, 4.6$. Blue lines follow magnetic field lines, ignoring $B_z$ (with a higher density of lines shown closer to the midplane). Panels (a,d,g,j,m) show $z=0$ (which looks fairly similar to all other $z$ values for this nearly 2-D case); (b,e,h,k,n) show $z=100 \sigma \rho _0\approx 0.2 L_z$ and (c,f,i,l,o) show $z=300 \sigma \rho _0 \approx 0.6L_z$.

For $L_z=L_x/32$, the picture is practically the same as in two dimensions: the initial current sheet tears, small plasmoids form and merge into larger plasmoids, until there is just one plasmoid (and a monster one at that) in the layer. Between plasmoids, the layer is thin at X-points (really, $X$-lines extended in the $z$ direction) where reconnection takes place. For $L_z=L_x$, the picture is initially very similar to that in two dimensions, when viewed in an $x$–$y$ slice. For example, at $t=0.3L_x/c$, a similar number of plasmoids is visible (although a slice at any value of $z$ looks qualitatively similar, these structures do not appear to extend uniformly in $z$; i.e. at early times we do not observe long flux ropes of length $L_z$). Soon (by $t=0.6L_x/c$), however, the 3-D layer shows general thickening without the clear plasmoid structures so familiar in two dimensions. At some rare places (e.g. $x\approx z \approx 100 \sigma \rho _0$, $tc/L_x=1.1$–$2.3$), the layer remains very thin, as at X-points in two dimensions. In the vicinity of these thin places we still see some ballooning of the layer, resembling a larger plasmoid, but smaller (in $y$) than in two dimensions, and with much less distinct density structure. After $2.3L_x/c$, most of the layer is thick, but not as thick as plasmoids in two dimensions.

The $L_z=L_x/32$ simulation is fairly uniform in $z$, consistent with its behaviour resembling 2-D reconnection. The 3-D simulation, however, is very non-uniform, as we will see in a series of plots showing various fields at time $t=2.3L_x/c$, at the centre of the layer, i.e. at $y_c(x,z,t)$ such that $B_x(x,y_c,z,t)=0$.

Figure 17 shows various quantities at the layer central surface, $y_c(x,z,t)$ (cf. § 3.2), vs $x$ and $z$, for the $L_z=L_x$ case, at $t=2.3L_x/c$. We focus particularly on region $A$, a patch of size (at this time) $\simeq 100\sigma \rho _0$ located around $x\approx z \approx 100\sigma \rho _0$, which we believe to be an actively reconnecting region; we saw in figure 16 that it is a site where the layer is near its thinnest. Panel (a) shows the current density $J_z$ in the $z$-direction, which is indeed most prominent in this reconnecting region, consistent with the layer being thin. (We note that current sheets are locally not necessarily parallel to $z$, but are nearly parallel to $z$ in many places.) Panel (b) shows strong outflows from this region in the $\pm x$ directions, as one would expect from reconnection; similarly, region $B$, around $x\approx 320\sigma \rho _0$ and $z\approx 400\sigma \rho _0$ with size (at this time) similar to patch $A$, has a high current density and strong outflows. There are also smaller reconnecting regions – e.g. at $x\approx 260\sigma \rho _0$ and $z\approx 280\sigma \rho _0$. The ‘line’ between the outflows thus appears to be a stagnation line, which is commonly associated with an X-line in the classic 2-D picture of (symmetric) reconnection. Panel (c) shows $\|\boldsymbol {E}\|$, which is sizeable ($\sim 0.1 B_0$) over a significant region, especially in the neighbourhood of large, strong outflows. However, (panel d) the parallel electric field $|E_\||$ (parallel to the local $\boldsymbol {B}$) is largest only in a region much closer to the stagnation line, as expected. Panel (e) shows that $\boldsymbol {E}\boldsymbol {\cdot } \boldsymbol {J}$ has more strongly positive areas than negative, indicating energy transfer from electromagnetic fields to the plasma; this transfer occurs around X-lines and extended outflows, but as panel (f) shows, the region where $E_\| J_\|$ is most positive is limited to areas closer to the largest stagnation lines (in regions $A$ and $B$). Figure 17 thus suggests that reconnection is occurring in patches (such as regions $A$ and $B$), perhaps not with such neat geometry as in two dimensions, but exhibiting basic signatures, such as a thin layer with high current density, strong outflows, particle energization and indeed particle energization by parallel electric fields or by electric fields in a region of weak magnetic field.

Figure 17.

Field quantities vs $x$ and $z$ at the layer central surface, i.e. at $y=y_c(x,z,t)$ such that $B_x(x,y_c,z,t)=0$ (cf. § 3.2); $A$ and $B$ indicate – roughly – examples of actively reconnecting regions. (a) The current density in the $z$-direction, $J_z$, is highest in reconnecting regions, where the layer is thin. (b) Shows $\overline {v_x}=(v_{e,x} + v_{i,x})/c$, the average of electron and positron $x$-velocities, shows strong outflows from reconnection regions. (c) The electric field, $\|\boldsymbol {E}\|$ is largest in reconnecting regions and also in strong outflows. (d) Although small, $E_\|$, the magnitude of the component of $\boldsymbol {E}$ parallel to the local $\boldsymbol {B}$, is stronger in reconnecting regions. (e) Here, $\boldsymbol {E}\boldsymbol {\cdot } \boldsymbol {J}$ indicates dominant power transfer from fields to plasma, mostly in reconnecting and outflow regions. (f) Shows that $\boldsymbol {E}_\|\boldsymbol {\cdot }\boldsymbol {J}_\|$ is positive and strongest in reconnecting regions, indicating significant plasma energization by parallel electric fields (or electric fields where $B\approx 0$). All plotted quantities were averaged over $y_c-\sigma \rho _0 < y < y_c+\sigma \rho _0$.

Very early in time, the layer exhibits many small patches of active reconnection. At any point in time, a patch tends to have roughly the same size in $x$ and $z$ (the $x$-extent can be seen in figure 16). For example, patches have size $\sim 20\sigma \rho _0$ at $t=0.3L_x/c$. It is perhaps therefore not surprising that the simulation with $L_z=L_x/32=16\sigma \rho _0$ – which is smaller than the patch size even at this early time – behaves like 2-D reconnection. Over time (for $L_z=L_x$), the number of reconnecting patches decreases rapidly, but the surviving patches grow in size, in both $x$ and $z$. By $2.3L_x/c$, the 3-D simulation has 2 major patches (regions $A$ and $B$), with perhaps a few much smaller, probably less active ones. By $4.6L_x/c$, only region $A$ remains, although it seems to have shrunk in the $z$ direction, perhaps spreading its current density a little wider in $x$. Importantly, region $A$ does continue to exhibit signs of active reconnection at these late-stage times, when the global reconnection rate is relatively slow in three dimensions. We leave a systematic investigation of the behaviour and statistics of these patches to future work.

Although the above-described patches in three dimensions locally resemble elementary reconnecting current sheets in two dimensions, important differences emerge away from these patches. As previously mentioned, we do not see the long-term storage of magnetic energy in highly structured plasmoids as in two dimensions; for example, in 3-D reconnection we do not see clear signatures of persistent, growing structures extended across the simulation in the $z$-direction (as we would expect flux ropes to be), e.g. in figure 17. In (nearly) two dimensions, the magnetic field is somewhat diminished far upstream of major X-lines, e.g. to $\approx 0.7 B_0$ more than $100\sigma \rho _0$ upstream at $t=4.6L_x/c$; at the same time, the magnetic field rarely falls below $0.3B_0$ anywhere (figure 18a). In two dimensions, plasmoids continue to grow roughly up to the size of the simulation, when the magnetic configuration prevents further reconnection. In three dimensions, however, the current layer becomes a somewhat thickened turbulent region with relatively sharp boundaries (where most current flows) inside which the magnetic field is generally diminished. The magnetic field is practically full strength outside of the sharp boundaries, while being less than $0.3 B_0$ in a large fraction of the volume within the boundaries (figure 18b,c).

Figure 18.

The magnetic field strength $|\boldsymbol {B}|$ at $t=4.6L_x/c$ in the same $z$-planes shown in figure 16; (a) nearly two dimensions, $z=0$, (b) three dimensions, $z=100\sigma \rho _0$ and (c) three dimensions, $z=300\sigma \rho _0$. Contours are drawn for $B/B_0=\sqrt {0.1}\approx 0.32$ (dotted white), $\sqrt {0.3}\approx 0.55$ (dashed green) and $\sqrt {0.5}\approx 0.71$ (solid blue).

The evolution of current sheets in three dimensions is complicated, and demands much more study. For now, however, we will move on to consider one of the most important consequences of current sheet evolution and resulting energy conversion, namely NTPA.

Reconnection-driven NTPA is a promising candidate to explain, e.g. high-energy X-ray and gamma-ray emission from a variety of astrophysical sources. It has been previously suggested that NTPA in three dimensions might be suppressed by the RDKI interfering with reconnection (Zenitani & Hoshino Reference Zenitani and Hoshino2007, Reference Zenitani and Hoshino2008). And since the reconnection electric field responsible for accelerating particles is proportional to the reconnection rate, and we have shown here that the global reconnection rate is slower in three dimensions than in two dimensions, we have another reason to suspect that NTPA might be considerably less efficient in three dimensions. However, that is emphatically not the case here. Figure 19 shows electron energy spectra $f(\gamma )$ (compensated by $\gamma ^2$) – compared both at the same time ($t=9 L_x/c$) and after the same amount of magnetic energy has been given to particles ($\Delta U_\textrm {plasma} \approx 0.17 \: U_{B0}$). Remarkably, all the spectra are pretty similar, exhibiting high-energy non-thermal behaviour with a power-law slope $\approx -4$ ( i.e. $f(\gamma )\sim \gamma ^{-4}$). However, though the difference may be small, 3-D reconnection yields more particles with very high energies than two dimensions; in addition, since ultimately 3-D reconnection can convert more magnetic energy to plasma energy than two dimensions, 3-D reconnection may accelerate more particles overall (not just at the very highest energies).

Figure 19.

Electron energy spectra $f(\gamma )$, compensated by $\gamma ^2$, where particle energy is $\gamma m_e c^2$, for $L_z/L_x=0$, 1/32, 1/8, 1, (a) at the same time, $t=9 L_x/c$ and (b) when the same amount of magnetic energy ($\simeq 17$ per cent) has been depleted. Dotted grey-on-yellow lines indicate a range of power-law slopes $f(\gamma )\sim \gamma ^{-3.4}$, $\gamma ^{-4}$ and $\gamma ^{-4.4}$.

NTPA continues to occur (in three dimensions) even during the second stage of reconnection, when the global rates of reconnection and magnetic energy conversion are smaller by an order of magnitude compared with the first stage. Figure 20 shows the spectra for $L_z=L_x$ at four different times, $tc/L_x=1.1, 2.3, 4.3, 8.6$, chosen so that between successive times, $U_\textrm {plasma}$ increases by equal amounts, namely $0.04 U_{B0}$ – i.e. $\Delta U_\textrm {plasma}/U_{B0} = 0.05, 0.09, 0.13, 0.17$. Even after $t\simeq 3L_x/c$, when the magnetic energy is at 89 per cent and after which the reconnection rate drops significantly, we see that $f(\gamma )$ continues to increase at high energies. Indeed, the increase in $f(\gamma )$ preserves the power-law slope at high energies; the high-energy slope of $f(\gamma )$ does not decrease after $3L_x/c$, nor does the high-energy cutoff decrease (if anything, the opposite occurs, possibly indicating re-acceleration of already-energetic particles that have returned coincidentally to actively reconnecting regions).

Figure 20.

Electron energy spectra $f(\gamma )$, compensated by $\gamma ^2$, at different times for the 3-D case $L_z/L_x=1$: the spectra are shown at times $tc/L_x=1.1$, 2.2, 4.3 and 8.6, when fractions 0.05, 0.09, 0.13 and 0.17 of the initial magnetic energy $U_{B0}$ have been depleted. Power-law indices of 3.4, 4 and 4.4 are shown for reference, as in figure 19. Even though reconnection is much slower after $4L_x/c$ than before, significant NTPA occurs between $tc/L_x=$4.3 and 8.6.

It is possible that NTPA in three dimensions is at least as efficient as in two dimensions, despite the slower (global) reconnection rate in three dimensions, because in three dimensions the local reconnection rate is as high as in two dimensions within actively reconnecting patches. Thus NTPA may be occurring in a similar manner in two dimensions and in three dimensions, but in a relatively smaller region in three dimensions. This result is consistent with Yin et al. (Reference Yin, Daughton, Karimabadi, Albright, Bowers and Margulies2008); the global reconnection rate drops because there are fewer/smaller patches, not because local reconnection rates drop.

Furthermore, the stark contrast in plasmoid evolution in two dimensions and three dimensions might promote more efficient NTPA in three dimensions. In two dimensions, ballooning plasmoids expand deeply into the upstream, and can contain higher magnetic field than the upstream region; the magnetic fields of large plasmoids can trap high-energy particles (with Larmor radii a large fraction of the plasmoid size; Sironi et al. Reference Sironi, Giannios and Petropoulou2016; Uzdensky Reference Uzdensky2020). In three dimensions, however, any structures (flux ropes) formed from reconnected magnetic field decay (see also Zhou, Loureiro & Uzdensky Reference Zhou, Loureiro and Uzdensky2020); high-energy particles may be confined to the layer by the upstream magnetic fields, but probably not by distinct sub-structures within the current layer. The lack of trapping may result in more particle acceleration in three dimensions; whereas a particle trapped in a plasmoid does not get accelerated by reconnection (but could be accelerated by a different, slower mechanism; Petropoulou & Sironi Reference Petropoulou and Sironi2018), untrapped particles could potentially pass multiple times through reconnecting layers, thus experiencing a sort of multi-stage acceleration (Dahlin et al. Reference Dahlin, Drake and Swisdak2015; Zhang et al. Reference Zhang, Sironi and Giannios2021).

It is worth estimating, from a Hillas criterion perspective (Hillas Reference Hillas1984), the energy gain of a particle in a typical reconnection electric field $E_\textrm {rec}\sim 0.04B_0$, traversing an entire patch of active reconnection. By $t=2.3L_x/c$, there is a patch (region $A$ in figure 17) in the simulation with $L_z=L_x$ that is about $\ell \sim 100\sigma \rho _0$ in size. Thus we estimate an energy gain $eE_\textrm {rec}\ell \sim 4 \sigma m_e c^2$, which is roughly consistent with the high-energy cutoff (a similar result was given for elementary current sheets in two dimensions, for high-$\sigma _h$ reconnection; Werner et al. Reference Werner, Uzdensky, Cerutti, Nalewajko and Begelman2016). We do not suggest this proves that all particles are being accelerated directly by the reconnection electric field; the acceleration mechanism is beyond the scope of this paper and remains under debate (see, e.g. Drake et al. Reference Drake, Swisdak, Che and Shay2006; Nalewajko et al. Reference Nalewajko, Uzdensky, Cerutti, Werner and Begelman2015; Ball et al. Reference Ball, Sironi and Özel2018, Reference Ball, Sironi and Özel2019; Guo et al. Reference Guo, Li, Daughton, Kilian, Li, Liu, Yan and Ma2019, Reference Guo, Li, Daughton, Li, Kilian, Liu, Zhang and Zhang2020; Kilian et al. Reference Kilian, Li, Guo and Li2020; Sironi & Beloborodov Reference Sironi and Beloborodov2020; Uzdensky Reference Uzdensky2020; Zhang et al. Reference Zhang, Sironi and Giannios2021), but at least for 2-D and 3-D reconnection with weak guide field, above-average-energy particles probably have gyroradii larger than the layer thickness and will experience more acceleration by the motional electric field of plasmoids, $\boldsymbol {E} = -(1/c)\boldsymbol {v}\times \boldsymbol {B}$ (Guo et al. Reference Guo, Li, Daughton, Kilian, Li, Liu, Yan and Ma2019, Reference Guo, Li, Daughton, Li, Kilian, Liu, Zhang and Zhang2020; Kilian et al. Reference Kilian, Li, Guo and Li2020; Uzdensky Reference Uzdensky2020). We offer this estimate simply as a way to compare system size with particle energies.

Thus we have seen both striking similarities and striking differences between 2-D and 3-D reconnection. Both convert significant amounts of magnetic energy to plasma energy, but the conversion rate is slower in three dimensions. Many of the signatures of reconnection at 2-D X-points (or X-lines) are apparent in 3-D simulations, but what happens in reconnection outflows is very different. While 2-D simulations beget a growing, highly structured hierarchy of distinct plasmoids containing energized particles and reconnected magnetic field, 3-D simulations only form small plasmoid-like features, which subsequently decay, allowing magnetic field energy in the layer to be converted to plasma energy. As a result, 3-D simulations ultimately convert more magnetic energy than in two dimensions. The decay of plasmoids in three dimensions means that the reconnected flux in those plasmoids is annihilated, while in two dimensions, flux is almost perfectly conserved, such that all the lost upstream flux ends up as reconnected flux inside plasmoids. In three dimensions, the lost upstream flux is mostly annihilated (and it is yet to be determined what fraction is directly annihilated and what fraction is first reconnected and later annihilated). Interestingly, perhaps the strongest similarity between 2-D and 3-D reconnection is NTPA; the spectrum of high-energy particles is very similar in both cases, especially if compared at times when the depleted magnetic energy is the same. With this overview finished, we will continue to examine differences between 2-D and 3-D reconnection in more depth.

5.2. Three dimensional reconnection: stochastic variability

An important difference between 2-D and 3-D current sheet simulations is the greater stochastic variability in three dimensions. When initializing the simulations, the initial electric and magnetic fields are precisely determined. However, the particle positions and velocities representing the plasma distribution in Monte Carlo fashion are chosen randomly (according to a precisely given distribution). Thus two ‘identical’ simulations will not be exactly the same, and the initial microscopic randomness can lead to variability at macroscopic scales.

In 2-D reconnection, two ‘identical’ simulations will yield different plasmoid formation, motion, and merging; this can affect the global reconnection and energy conversion rates. Smaller simulations often show greater stochastic variation in global quantities such as energy vs time, because there are fewer plasmoids overall, increasing the influence of erratic individual plasmoids. We have often observed that the process of plasmoid merging (especially in simulations without initial perturbation; cf. § 4.2) results in two monstrously large plasmoids in a layer spaced apart by nearly $L_x/2$. These two plasmoids will ultimately merge; but by symmetry due to periodic boundary conditions, one plasmoid can either move left or right to merge with the other. This symmetry is always eventually broken to yield a lower-energy magnetic configuration, but its breaking can be ultimately traced to the initial microscopic randomness. As these plasmoids decide which way to go for the final merger, the magnetic energy often stagnates at a roughly constant level for a short time before rapidly declining during the final plasmoid merger (see, e.g. figure 7, $L_x/\sigma \rho _0=640$, 1280 and 2560, around 3–5$L_x/c$). However, due to randomness, sometimes the final two plasmoids are not created so symmetrically, and the magnetic energy decline continues smoothly. Thus we often describe reconnection as ‘bursty’; in radiative reconnection, where high-energy particles cool rapidly, e.g. via inverse Compton radiation, this burstiness is visible even in the spectrum of high-energy particles (Werner et al. Reference Werner, Philippov and Uzdensky2019).

Thus initial random particle positions and velocities lead to macroscopic variability in two dimensions at system-size scales. However, although significant, this stochastic variability is a sequence of random short detours on a path to essentially the same final state. In fact, § 4.2 showed that varying the initial current sheet configuration (varying more than just different random velocities from the same distribution) still had little effect on the overall evolution of reconnection. Looking at coarse time scales of $O(L_x/c)$, the overall rate and amount of energy depletion are fairly independent of initial conditions (for the same upstream conditions).

In contrast, in a 3-D reconnection simulation, random initialization can lead to significant macroscopic differences that persist for long times – at least as long as we have been able to simulate. While a detailed study of random variation in large ensembles of ‘identical’ simulations is beyond the scope of this paper, we present here an anecdote of high variability in three dimensions in a case with $L_x=256 \sigma \rho _0$ and $L_z=L_x/4$, which we ran three times with ‘identical’ set-ups ($B_{gz}=0$, $\eta =5$, $\beta _d=0.3$, $a=0$). The three simulations ultimately exhibit very different magnetic energy and flux evolution (figure 21) even over scales of tens of $L_x/c$; this variation is much greater than we see in two dimensions (although not perfectly comparable, figures 3, 5 and 7 place a rough upper bound on variability in two dimensions; also see figure 22). Nevertheless, despite these differences, the particle energy spectra are similar despite different energy evolution (figure 21c), with the exception of one case exhibiting a bump above $\gamma \gtrsim 10\sigma$. At early times, the evolutions of magnetic energy and flux (figure 21a,b) are reasonably similar, but they diverge after $\simeq 5L_x/c$. Simulation (A) starts off with the highest magnetic conversion rate but ends up converting the least amount (30 per cent) of magnetic energy after $30L_x/c$. Simulation (B) rapidly converts approximately 40 per cent of the magnetic energy within $10L_x/c$, and then continues to convert energy at a much slower rate for at least the next $20L_x/c$. And simulation (C) exhibits a medium-slow stage of energy conversion between 5 and 15$L_x/c$, followed by a stagnant or very slow final stage similar to the others.

Figure 21.

Three-dimensional reconnection exhibits significant stochastic variability, as shown by three simulations (A, B, C) with $L_x=256\sigma \rho _0$ and $L_z=L_x/4$, with otherwise identical parameters, but different initial random particle positions and velocities. Panels show: transverse magnetic energy (a) and unreconnected flux (b) vs time, and (c) particle energy spectra at $20L_x/c$, compensated by $\gamma ^4$.