© ESO, 2014

1. Introduction

Magnetic reconnection has been the focus of extended studies since its first introduction by Giovanelli (1947, 1948) to explain the sudden release of energy in solar flares. The term itself was coined by Dungey (1958). It is now the key ingredient for theories of coronal heating, solar flares and jets, and coronal mass ejections in the Sun (Priest 1987), of magnetic storms and substorms in the Earth magnetosphere (Paschmann et al. 2013), and for the behavior of fusion plasmas with, for instance, the sawtooth oscillation in tokamaks (Biskamp 2000). Space physics proves that magnetic reconnection can quickly convert magnetic energy into kinetic energies (bulk flow, heat, non-thermal particles) with fast variability and high efficiency.

Such attributes made it attractive for high-energy astrophysics to explain, for example, radiation (Romanova & Lovelace 1992) and flares (Giannios et al. 2009) in active galactic nuclei (AGNs) jets or in gamma-ray bursts (Lyutikov 2006a; Lazar et al. 2009), the heating of AGN and microquasar coronae and associated flares (Di Matteo 1998; Merloni & Fabian 2001; Goodman & Uzdensky 2008; Reis & Miller 2013), the flat radio spectra from galactic nuclei and AGNs (Birk et al. 2001), the heating of the lobes of giant radio galaxies (Kronberg et al. 2004), the σ-paradox and particle acceleration at pulsar wind termination shocks (Kirk & Skjæraasen 2003; Pétri & Lyubarsky 2007; Sironi & Spitkovsky 2011), GeV flares from the Crab nebula (Cerutti et al. 2012a,b, 2013), transient outflow production in microquasars and quasars (de Gouveia dal Pino & Lazarian 2005; de Gouveia Dal Pino et al. 2010; Kowal et al. 2011; McKinney et al. 2012; Dexter et al. 2014), gamma-ray burst outflows and non-thermal emissions (Drenkhahn & Spruit 2002; McKinney & Uzdensky 2012), X-ray flashes (Drenkhahn & Spruit 2002), soft gamma-ray repeaters (Lyutikov 2006b; Uzdensky 2011), flares in double pulsar systems (Lyutikov & Lazarian 2013), or energy extraction in the ergosphere of black holes (Koide & Arai 2008).

As pointed out by Uzdensky (2006), magnetic reconnection is of dynamical importance in any environment both where magnetic fields dominate the energy budget, so that the energy transfer can have dynamical and observable consequences, and where the rates of reconnection are fast. The latter is known to hold both in collisionless plasmas (Birn et al. 2001) and in collisional plasmas, either via turbulence (Lazarian & Vishniac 1999; Lapenta & Lazarian 2012; Lazarian et al. 2012), or via plasmoid induced reconnection (Uzdensky et al. 2010; Zanotti & Dumbser 2011; Loureiro et al. 2007, 2012).

Many of these environments are collisionless (Ji & Daughton 2011), so that fast reconnection must be triggered and sustained by non-ideal terms other than collisional ones, which implies kinetic processes on scales close to the electron inertial length or Larmor radius, with particles largely out of equilibrium and possibly comprising high-energy tails. These non-ideal terms can be linked to particle inertia and wave-particle resonant interactions or to finite Larmor radius effects in magnetic field gradients. Simulation studies thus require full kinetic codes, such as Vlasov solvers or particle-in-cell algorithms.

Most of the above environments are also relativistic, either because of relativistic velocities (bulk flows or currents) or because the thermal kinetic energy and/or the magnetic energy density exceeds the restmass energy of the particles. The latter translates into an inflow magnetization larger than unity, ${\mathit{\sigma }}_{\mathrm{in}\mathit{,s}}\mathrm{=}{\mathit{B}}_{\mathrm{in}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}{\mathit{\mu }}_{\mathrm{0}}{\mathit{n}}_{\mathrm{in}}{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}\mathrm{\right)}\mathit{>}\mathrm{1}$, with s denoting ions or electrons. This magnetic energy can be transferred to the particles, and because it is greater than the particles restmass, relativistic particles are expected. The relation h_0,out,sΓ_out,s = 1 + σ_in,s, with h_0,out,s the enthalpy and Γ_out,s the bulk Lorentz factor of the reconnection outflow (see Sect. 3.4, Eq. (11)), indeed shows that either relativistic temperatures (h_0,out,s> 1) or relativistic bulk velocities (Γ_out,s> 1) are obtained for the outflows. The relevant magnetization is thus not that of the plasma, which is low because of the ion mass, but that of each species taken individually.

Studies of relativistic reconnection are scarcer than their non-relativistic counterparts (for the latter, see the reviews by Birn & Priest 2007; Treumann & Baumjohann 2013), and they mainly deal with pair plasmas. For relativistic pair plasmas, they include analytical works (Tenbarge et al. 2010; Kojima et al. 2011) and Sweet-Parker-like analysis (Blackman & Field 1994; Lyutikov & Uzdensky 2003; Lyubarsky 2005), 2D MHD simulations (Watanabe & Yokoyama 2006; Zenitani et al. 2011a; Takahashi et al. 2011; Zanotti & Dumbser 2011; Takamoto 2013; Baty et al. 2013), two-fluid simulations (Zenitani et al. 2009a,b), test-particle simulations (Bulanov & Sasorov 1976; Romanova & Lovelace 1992; Larrabee et al. 2003; Cerutti et al. 2012a), 1D PIC simulations (Pétri & Lyubarsky 2007), 2D PIC simulations (Jaroschek et al. 2008; Sironi & Spitkovsky 2011; Bessho & Bhattacharjee 2012; Cerutti et al. 2012b, 2013; Zenitani & Hoshino 2001, 2005a, 2007, 2008; Zenitani & Hesse 2008), and 3D PIC simulations (Zenitani & Hoshino 2008, 2005b; Sironi & Spitkovsky 2011, 2014; Kagan et al. 2013; Cerutti et al. 2014). Relativistic reconnection in ion-electron plasmas has been studied less. We find a test-particle simulation (Romanova & Lovelace 1992), a resolution of the diffusion equation (Birk et al. 2001), and a discussion by Sakai et al. (2002) in a 2D PIC simulation of laser fusion beams.

The focus of the present work is on relativistic reconnection, as compared to non-relativistic studies, and on ion-electron plasmas, as compared to pair plasmas. Our goals are to carve out aspects that are particular to this regime, to shed light on the underlying physical causes, and to ultimately put our findings in the, admittedly speculative, larger astrophysical context of microquasar and AGN disk coronae and magnetospheres, and of other possible environments with ion-electron relativistic plasmas. Part of our results are also of interest for pair plasmas and for non-relativistic cases.

In Sect. 2 we describe the simulation setup and parameters. Section 3 presents the results of simulations with antiparallel asymptotic magnetic fields. We investigate the structure of the two-scale diffusion region in Sect. 3.2, and explain why we see sharp transitions at the entrance of this region. Next, we turn to the relativistic Ohm’s law. In non-relativistic reconnection, non-ideal terms are dominated by thermal inertia, i.e., by the divergence of off-diagonal elements of the pressure tensor. There are, however, PIC studies (see references of Sect. 3.3) suggesting that for relativistic reconnection, thermal inertia can be dominated by bulk inertia. In Sect. 3.3 we show that this is the case in our simulations with large inflow magnetization. We demonstrate in Sect. 3.8 that this is to be expected on the basis of an analytical model. Concerning the reconnection outflows, mass and energy conservation imply that relativistic inflow magnetization results in relativistic temperatures and/or relativistic bulk velocities in the outflows, but say nothing about the balance between the two. In Sect. 3.4 we show that in our simulations, thermal energy largely dominates the bulk kinetic energy. We demonstrate analytically in Sect. 3.8 that this is to be expected for large inflow magnetization, under the assumption that thermal inertia significantly contributes in Ohm’s law. This is an important question that has observational consequences. In Sect. 3.5 we detail the structure of the magnetic islands and of their central density dips and isolated centers. Section 3.6 studies the reconnection electric field. The relevant normalization is non-trivial for relativistic setups, and we propose to use the relativistic Alfvén speed in the inflow, which leads to rates in a close range. Section 4 highlights differences resulting from the presence of a guide magnetic field. We summarize and conclude our work in Sect. 5 and discuss applications to astrophysical objects.

2. Problem setup

2.1. Description of the relativistic Harris equilibrium

We use the explicit particle-in-cell code Apar-T, presented and tested in Melzani et al. (2013). Broadly speaking, it is a parallel electromagnetic relativistic 3D PIC code with a staggered grid, where the fields are integrated via Faraday and Maxwell-Ampère equations, currents computed by a charge conserving volume weighting (CIC), and fields interpolated accordingly.

The simulations start from a Harris equilibrium, which is a solution of the Vlasov-Maxwell system. The magnetic field is ${{B}}_{\mathrm{rec}}\mathrm{=}\hat{{z}}\hspace{0.17em}{\mathit{B}}_{\mathrm{0}}\tanh {}^{\mathrm{(}}\mathit{x}\mathit{/}{\mathit{L}}^{\mathrm{)}}$ (see Figs. 1 or 2 for axis orientation), and is sustained by a population of electrons and ions of equal number density n_cs(x) = n_cs(0) / cosh²(x/L) (cs stands for current sheet), flowing with bulk velocities U_e and U_i = −U_e in the ± y directions. We denote the associated Lorentz factors by Γ_e and Γ_i. Each species follows a Maxwell-Jüttner distribution (Eq. (A.1)) of normalized temperature Θ_s = 1 /μ_s = T_s/m_sc².

We derived the equilibrium relations for relativistic temperatures and current drift speeds, as well as for arbitrary ion to electron mass ratios and temperature ratios, in Melzani et al. (2013). Details specific to the present application can be found in Appendix A. The equilibrium depends on the ratio ω_ce/ω_pe, with ω_pe = (n_cs(0)e²/ (ϵ₀m_e))^{1 / 2} the electron plasma pulsation defined with the lab-frame number density n_cs(0) (not including background particles), and ω_ce = eB₀/m_e the electron cyclotron pulsation in the asymptotic magnetic field (not including the guide field, e> 0). Our simulations are, however, not loaded exactly with the equilibrium values, but with a temperature and current speed uniformly in excess of ~ 10% in order to shorten the otherwise rather long stable phase.

We also set a background plasma of number density (for electrons or for ions) n_bg and of temperature T_bg,e for electrons and T_bg,i for ions. Finally, a guide magnetic field is sometimes considered, i.e., a uniform component ${{B}}_{\mathrm{G}}\mathrm{=}{\mathit{B}}_{\mathrm{G}}\hat{{y}}$. Adding the background plasma or the guide field does not change the Harris equilibrium.

2.2. Magnetization and energy fluxes

There are several ways to characterize the magnetization of the configuration. The ratio ω_ce/ω_pe has no direct physical meaning and is mostly used as a simulation label. The magnetization ${\mathit{\sigma }}_{\mathit{s}}^{\mathrm{hot}}$ of the background plasma species s is the ratio of the energy flux in the reconnecting magnetic field to that in the particles (restmass, thermal, bulk). The electromagnetic energy flux is the Poynting flux. Far from the current sheet, it reads as $\frac{{E}\mathrm{\wedge }{B}}{{\mathit{\mu }}_{\mathrm{0}}}\mathrm{=}\frac{{E}\mathrm{\wedge }{{B}}_{\mathrm{rec}}}{{\mathit{\mu }}_{\mathrm{0}}}\mathrm{+}\frac{{E}\mathrm{\wedge }{{B}}_{\mathrm{G}}}{{\mathit{\mu }}_{\mathrm{0}}}\mathrm{=}\frac{{\mathit{B}}_{\mathrm{rec}}^{\mathrm{2}}}{{\mathit{\mu }}_{\mathrm{0}}}{\mathit{v}}_{{\mathit{E}}_{\mathit{y}}\mathrm{\times }{\mathit{B}}_{\mathrm{rec}}}\hat{{x}}\mathrm{+}\frac{{\mathit{B}}_{\mathrm{G}}^{\mathrm{2}}}{{\mathit{\mu }}_{\mathrm{0}}}{\mathit{v}}_{{\mathit{E}}_{\mathit{y}}\mathrm{\times }{\mathit{B}}_{\mathrm{rec}}}\hat{{x}}\mathit{,}$(1)where ${B}\mathrm{=}{\mathit{B}}_{\mathrm{rec}}\hat{{z}}\mathrm{+}{\mathit{B}}_{\mathrm{G}}\hat{{y}}$ and v_{Ey × Brec} = E_y/B_rec. This splitting of the energy flux into two contributions, one from the magnetic field that will reconnect, the other from the guide field that will mostly be compressed, is possible only if the electric field is normal to B, which is indeed the case in the ideal outer area because of the tendency of the plasma to screen parallel electric fields. For the particles, the energy flux of species s is (Eq. (B.9)) ${\mathit{n}}_{\mathrm{lab}\mathit{,s}}\mathrm{⟨}\mathit{v\gamma }{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{s}}\mathrm{=}{\mathit{n}}_{\mathrm{lab}\mathit{,s}}{\mathrm{\Gamma }}_{\mathit{s}}{\mathit{h}}_{\mathrm{0}\mathit{,s}}\mathit{v̅}\mathit{s}{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}\mathit{,}$(2)with n_lab,s the particle number density in the lab frame (=Γ_s times that in the comoving frame), ⟨ ·⟩ _s denoting an average over momentum of the distribution function, m_s the particle mass, their bulk velocity, Γ_s the associated Lorentz factor, and h_0,s their comoving enthalpy (drawn in Fig. B.1 for a thermal distribution). All in all, the magnetization of species s is $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{s}}^{\mathrm{hot}}\mathrm{\left(}{\mathit{B}}_{\mathrm{rec}}\mathrm{\right)}& \mathrm{=}& \frac{\mathit{E}\mathrm{\times }{\mathit{B}}_{\mathrm{rec}}\mathit{/}{\mathit{\mu }}_{\mathrm{0}}}{{\mathit{n}}_{\mathrm{lab}\mathit{,s}}\mathrm{⟨}\mathit{v\gamma }{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{s}}}\mathrm{=}\frac{{\mathit{B}}_{\mathrm{rec}}^{\mathrm{2}}\mathit{/}{\mathit{\mu }}_{\mathrm{0}}}{{\mathit{n}}_{\mathrm{lab}\mathit{,s}}{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}{\mathrm{\Gamma }}_{\mathit{s}}{\mathit{h}}_{\mathrm{0}\mathit{,s}}}\\ & \mathrm{=}& \end{array}$(3)with ${\mathit{\sigma }}_{\mathit{s}}^{\mathrm{cold}}$ the magnetization of the plasma without taking temperature effects and relativistic bulk motion into account: ${\mathit{\sigma }}_{\mathit{s}}^{\mathrm{cold}}\mathrm{\left(}\mathit{B}\mathrm{\right)}\mathrm{=}\frac{{\mathit{B}}^{\mathrm{2}}}{{\mathit{\mu }}_{\mathrm{0}}{\mathit{n}}_{\mathrm{lab}\mathit{,s}}{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}}\mathrm{·}$(4)If ${\mathit{\sigma }}_{\mathit{s}}^{\mathrm{cold}}\mathrm{\left(}{\mathit{B}}_{\mathrm{rec}}\mathrm{\right)}$ exceeds unity, then it is possible to pass to the particles an amount of energy from the reconnecting field that exceeds their restmass; i.e., it is possible to obtain relativistic particles. We do not include the guide field ${\mathit{B}}_{\mathrm{G}}\hat{{y}}$ in the definition of the magnetization because it is mostly compressed and does not transfer energy to the particles.

Finally, the total magnetization of the plasma is ${\mathit{\sigma }}_{\mathrm{i}\mathrm{+}\mathrm{e}}\mathrm{\left(}{\mathit{B}}_{\mathrm{rec}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{B}}_{\mathrm{rec}}^{\mathrm{2}}\mathit{/}{\mathit{\mu }}_{\mathrm{0}}}{{\sum }_{\mathit{s}}{\mathit{n}}_{\mathrm{lab}\mathit{,s}}{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}{\mathrm{\Gamma }}_{\mathit{s}}{\mathit{h}}_{\mathrm{0}\mathit{,s}}}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{i}}^{\mathrm{cold}}\mathrm{\left(}{\mathit{B}}_{\mathrm{rec}}\mathrm{\right)}}{{\sum }_{\mathit{s}}{\mathrm{\Gamma }}_{\mathit{s}}{\mathit{h}}_{\mathrm{0}\mathit{,s}}\mathrm{(}{\mathit{m}}_{\mathit{s}}\mathit{/}{\mathit{m}}_{\mathrm{i}}\mathrm{)}}\mathrm{·}$(5)In the inflow part of our simulations, we have h_0,i ~ 1, h_0,e < 2.6, and m_i = 25m_e for the range of background temperatures considered here, so that the particle energy flux is largely dominated by the restmass energy flux of the ions, which is ${\mathit{n}}_{\mathrm{lab}\mathit{,}\mathrm{i}}\mathit{v̅}{}_{\mathrm{i}}\mathit{m}_{\mathrm{i}}{\mathit{c}}^{\mathrm{2}}$, and has no temperature dependence. We thus have ${\mathit{\sigma }}_{\mathrm{i}\mathrm{+}\mathrm{e}}\mathrm{~}{\mathit{\sigma }}_{\mathrm{i}}^{\mathrm{cold}}$, and σ_{i + e} is not a good representative of the electron physics and of the possibility that they are relativistically magnetized. The inflow magnetizations in our simulations are presented in Table 2.

2.3. Alfvén velocities

We give the definitions of the Alfvén speeds that will be used to discuss the normalization of the reconnection electric field in Sect. 3.6. They are reported in Table 2. The Alfvén velocity in the inflow plasma, far from the current sheet, is expressed in the comoving plasma frame and is, respectively in the non-relativistic and relativistic cases, $\begin{array}{ccc}{\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{in}}^{\mathrm{NR}}& \mathrm{=}& \\ {\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{in}}^{\mathrm{R}}& \mathrm{=}& \end{array}$where σ_{i + e}(B_tot) is to be expressed in the comoving frame (Eq. (5) with n_lab,s = n_0,bg the comoving density and Γ_s = 1), and where ${\mathit{B}}_{\mathrm{tot}}\mathrm{=}\mathrm{(}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{+}{\mathit{B}}_{\mathrm{G}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$. For the relativistic expression 6b, the first equality is general and derived from the relativistic ideal MHD description (Gedalin 1993), while the second holds only because m_i ≫ m_e so that the total enthalpy is dominated by the ion contribution. When there is a guide magnetic field, we show in Sect. 4.3 that it is relevant to project the Alfvén speed into the direction of the reconnecting magnetic field ($\hat{{z}}$), i.e., to consider ${\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{in}}^{\mathrm{R}}\cos \mathit{\theta }$ with tanθ = B_G/B₀.

A hybrid Alfvén speed is often defined in the literature as depending on the asymptotic magnetic field (without the guide field) and on the comoving density at the center of the current sheet: ${\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{0}}^{\mathrm{NR}}\mathrm{=}\frac{{\mathit{B}}_{\mathrm{0}}}{\sqrt{{\mathit{\mu }}_{\mathrm{0}}{\mathit{n}}_{\mathrm{0}\mathit{,}\mathrm{cs}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{(}{\mathit{m}}_{\mathrm{e}}\mathrm{+}{\mathit{m}}_{\mathrm{i}}\mathrm{)}}}\mathrm{=}{\left(\frac{{\mathit{m}}_{\mathrm{e}}}{{\mathit{m}}_{\mathrm{e}}\mathrm{+}{\mathit{m}}_{\mathrm{i}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\frac{{\mathit{\omega }}_{\mathrm{ce}}}{{\mathit{\omega }}_{\mathrm{0}\mathit{,}\mathrm{pe}}}\hspace{0.17em}\mathit{c,}$(7)where a subscript 0 indicates a comoving quantity. Its relativistic generalization is denoted by ${\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{0}}^{\mathrm{R}}$ and is obtained with Eq. (6b) but with parameters of the plasma at the center of the current sheet in the magnetization.

2.4. Simulation parameters and resolution tests

The physical parameters of the main simulations are given in Tables 1 and 2. We consider a mass ratio m_i/m_e = 25, except for one simulation with pairs. The background plasma number density is n_bg = 0.1n_cs(0) or 0.3n_cs(0). Its temperature is varied between T_bg = 1.5 × 10⁷K (2.5 × 10^-3m_ec²) and 3 × 10⁹K (0.5 m_ec²). The magnetization depends on the ratio ω_ce/ω_pe = 1, 3, or 6, leading to inflow magnetizations ${\mathit{\sigma }}_{\mathit{s}}^{\mathrm{hot}}$ between 10 and 260 for electrons, or 0.4 and 14 for ions. The current sheet is either of initial halfwidth L = 0.5d_i = 2.5d_e or L = 1d_i in the ω_ce/ω_pe = 6 case, with d_e, d_i the inertial lengths defined at current sheet center. We stress that for relativistic temperatures the sheet width in terms of Larmor radii will not be the same for simulations with different ω_ce/ω_pe (Eq. (A.5)); see L/r_ce in Table 1.

The numerical resolution is set by the number of cells n_x per electron inertial length d_e, by the number of timesteps n_t per electron plasma period 2π/ω_pe, and by the number of computer particles (the so-called superparticles) per cell ρ_sp. The quantities d_e, ω_pe, and ρ_sp are defined at t = 0 at the center of the current sheet, where the particle density is highest. For the simulations of Table 2, we take n_x = 9 and n_t = 150 (250 for ω_ce/ω_pe = 6). We checked with a simulation with twice this resolution (ω_ce/ω_pe = 3, n_x = 18, n_t = 250) that all of the presented results are not affected¹.

Concerning the number of superparticles per cell, the simulations of Table 2 use ρ_sp = 1820 (1090 for n_bg/n_cs(0) = 0.3). For the case n_bg/n_cs(0) = 0.1, this corresponds to 1650 electron and ion superparticles per cell for the plasma of the current sheet and to 170 for the background plasma. The density profile of the current sheet plasma is set by changing the number of superparticles per cell when going away from the center. We stressed in Melzani et al. (2013, 2014a) that because of their low numbers of superparticles per cell when compared to real plasmas, PIC simulations present high levels of collisionality. One should thus ensure that collisionless kinetic processes remain faster than collisional effects (e.g., for thermalization), essentially by taking a large enough number Λ^PIC of superparticles per Debye sphere. For example, with Θ_e = 2.4 the electron Debye length is 20 cells large, and we have initially at the center of the current sheet: Λ^PIC ~ 364 × 20 × 20 = 7.3 × 10⁵ superparticles. For a background plasma with T_bg = 2 × 10⁸ K: Λ^PIC = 133. We performed a simulation identical to the one with ω_ce/ω_pe = 3, n_bg/n_cs(0) = 0.1, T_bg = 2 × 10⁸ K, but with half the superparticles, and found no change in the results. It shows that the main simulations use large enough ρ_sp.

Boundaries are periodic along z and y. At the top and bottom x boundaries, we use reflective boundaries; i.e., we place a perfectly conducting wall that reflects waves and particles. The number of cells for the standard simulations is 4100 × 6144. The length along y is of no dynamical importance, and the dimensions correspond to a 2D simulation with 455 initial electron inertial lengths along x and 683 along z, typically with 4 × 10⁹ superparticles. It takes 70 T_pe for light waves to start from the current sheet, reflect at the ± x boundaries, and come back to the sheet. for, respectively, simulations with ω_ce/ω_pe = 1, 3, 6. The light travel time in the z direction is longer. Except for run ω_ce/ω_pe = 1, all the analyses presented here are for shorter times and are thus not affected by boundaries. To check this, we performed a larger simulation with 8000 × 10 240 cells (i.e., 888 × 1138 initial electron inertial length) for the case ω_ce/ω_pe = 3, n_bg/n_cs(0) = 0.1, T_bg = 2 × 10⁸ K, with the same n_t, n_x, ρ_sp. The corresponding light-crossing time is now $\mathrm{136}\hspace{0.17em}{\mathit{T}}_{\mathrm{pe}}\mathrm{=}\mathrm{101}{\mathit{\omega }}_{\mathrm{ci}}^{-1}\mathrm{=}\mathrm{2535}{\mathit{\omega }}_{\mathrm{ce}}^{-1}$. All the results are the same, which shows that we do not suffer from boundary effects.

3. Results with no guide field

This section explores results for simulations with no guide field, where the magnetic field above and below the current sheet is antiparallel.

3.1. Overall structure and evolution

We first summarize the general picture, raising important points that are explained in the next sections. The initial kinetic equilibrium is unstable to the collisionless tearing mode, which in all presented simulations is triggered by the noise level and by the slightly out-of-equilibrium initial state. We have studied in detail the linear phase of this instability with PIC simulations (Melzani et al. 2013) for pair plasmas and found growth rates within 5% of the analytical derivations of Pétri & Kirk (2007) made on the basis of a linearization of the Vlasov-Maxwell system. Physically, the instability is driven by the particles freely bouncing in the layer of magnetic field reversal (Coppi et al. 1966), with a mechanism similar to a filamentation instability: perturbations in B_x and B_z lead to a bunching of the particles, which in turn increase the magnetic field perturbation. It leads to the formation of alternating X and O points, here with no privileged location because we impose no localized initial perturbation (Fig. 1).

With the appearance of X and O points, the magnetic flux variations across fixed contours induce an out-of-plane electric field ${\mathit{E}}_{\mathit{y}}\hat{{y}}\mathrm{\propto }\mathrm{-}\hat{{y}}$, which amplifies the initial current along −y, which in turn increases the magnetic field in order to cancel out the former magnetic flux variations and prevent reconnection. However, non-ideal processes forbid an ideal plasma response and allow the triggering of reconnection and the existence of a finite electric field at the current sheet center, where ideal Ohm’s law would otherwise read E_y = 0. This electric field E_y is at the heart of the reconnection process, because it is responsible for the transfer of energy from the magnetic field to particles in the diffusion region. We detail Ohm’s law and the contribution of each non-ideal term in Sect. 3.3.

Plasma and magnetic fields decouple in the non-ideal region, and flux tubes can “reconnect”, producing new flux tubes strongly bent that accelerate the plasma outward in the $\mathrm{\pm }\hat{{z}}$ directions, thus producing the exhaust outflows. This depletion of particles and/or the spreading of the electric field E_y outside the current sheet create an inflow from the ideal zone toward the current sheet: particles E × B drift at a speed ${{v}}_{\mathrm{in}}\mathrm{=}{E}\mathrm{\times }{B}\mathit{/}{\mathit{B}}^{\mathrm{2}}\mathrm{=}{\mathit{E}}_{\mathit{y}}\mathit{/}{\mathit{B}}_{\mathit{z}}\hat{{x}}$. The incoming particles are then accelerated along $\hat{{y}}$ by E_y once they enter the non-ideal region where they are unmagnetized (because there E>B). The structure of this central region is investigated in Sect. 3.2.

The exhaust outflows, which in the MHD view are driven by the magnetic field tension force, are produced by particles accelerated by E_y along $\hat{{y}}$ and then slowly rotating owing to the increasingly strong magnetic field B_x as one goes away from the X point (Fig. 3). As they do so, particles still gain energy as long as ${\mathit{q}}_{\mathit{s}}{\mathit{E}}_{\mathit{y}}\hat{{y}}\mathrm{·}{v}\mathit{>}\mathrm{0}$. The energy content of these outflows comprises a Poynting flux, bulk kinetic energy, and thermal energy, with respective weights that depend on the background plasma parameters, as studied in Sect. 3.4. The balance between inflow and outflow can lead to a steady state Sweet-Parker-like configuration, depending on the conditions.

In the configuration of the simulations the initial perturbation is not localized in space, so that several X points appear, with islands inbetween that collect the flux of particles and of the reconnected magnetic field. The islands are trapped between two exhausts, and the bulk energy of the outflows is converted into heat by random scatterings in the complex electric and magnetic structures at the island entrance and inside the islands (Fig. 10), however with particle distributions that are not necessarily thermal (see Sect. 3.5). The islands grow and, since they are threaded by parallel currents (along $\mathrm{-}\hat{{y}}$), attract each other via the Lorentz force and merge, thus growing even more. As time goes on, the island number dwindles and the space inbetween them increases, forming elongated current sheets composed of a X point surrounded by two elongated exhausts (see Fig. 2). We stop the simulations when only two or three islands remain.

3.2. Inflow: two-scale diffusion region and sharp transitions

We now examine the inflow of plasma into the diffusion region. In the literature, for antiparallel reconnection (i.e., no guide field), its width is found to scale with the particles inertial length, a result that we show to also hold for relativistic reconnection in Sect. 3.2.1. The originality of the following results is the formation of a very sharp transition at the entrance of the diffusion regions, which we explore in Sect. 3.2.2.

3.2.1. Width of the diffusion region

The diffusion region is, by definition, the area where impeding mechanisms (which can be collisions, inertia and collective interactions, or finite Larmor effects) prevent the particles from responding in an ideal way to the electric fields induced by magnetic flux variations. The magnetic field and the plasma are then not coupled anymore, the former can freely diffuse, and reconnection can start or be sustained. Defining the diffusion region is thus a matter of finding the area where the non-ideal processes dominate ideal behavior.

Defining unambiguous criteria to identify this region is a complex subject (Ishizawa & Horiuchi 2005; Klimas et al. 2010), especially in the presence of a guide field (Hesse et al. 2002, 2004; Liu et al. 2014) in 3D simulations (Pritchett 2013) or in asymmetric configurations (Zenitani et al. 2011b). In our case, as we show in Sect. 3.2.2, there is a sharp increase in particle density when the inflow plasma reaches the central part, where particles are retained by bouncing motions around the reversing magnetic field. It is associated with a sharp drop in inflow velocity, a rise in temperature, and a violation of the frozen-in relation ${E}\mathrm{+}\overline{{v}}{}_{\mathit{s}}\mathrm{\wedge }{B}\mathrm{=}\mathrm{0}$. We identify the diffusion region with this area of increased density.

Because of their lighter mass and fastest response, electrons remain frozen to the magnetic field longer than ions. The ion diffusion region is thus larger than the electron one. In all the antiparallel simulations, we find the expected two-scale (ion and electron) diffusion region. It can be seen in Fig. 4, where we present a cut along x through an X point at a given time for different simulations. The width δ_s of the diffusion region is roughly given by the inertial length d_s of the corresponding species, defined with the particle density at the center of the diffusion region, i.e., ${\mathit{\delta }}_{\mathit{s}}\mathrm{~}{\mathit{d}}_{\mathit{s}}\mathrm{=}\mathit{c}\mathit{/}\sqrt{{\mathit{n}}_{\mathit{s}}\mathrm{(}\mathit{x}\mathrm{=}\mathrm{0}\mathit{,t}\mathrm{)}{\mathit{q}}_{\mathit{s}}^{\mathrm{2}}\mathit{/}{\mathit{ϵ}}_{\mathrm{0}}{\mathit{m}}_{\mathit{s}}}$, with throughout all simulations 0.5 ≤ δ_i/d_i ≤ 1 and 1 ≤ δ_e/d_e ≤ 1.5 (Fig. 5). This is also the scaling found in PIC simulations of non-relativistic ion-electron magnetic reconnection (e.g., Pritchett 2001; Klimas et al. 2010).

In the case of hot background electrons (T_bg,e = 3 × 10⁹ K, with a corresponding background plasma β_e = 1.1 × 10^-2), the transition is less sharp, and the width is larger than the inertial length. The sharpness is discussed in Sect. 3.2.2, and the larger extent is expected because inflowing particles have faster speeds, and thus larger Larmor radii and larger bouncing motions. The width δ_s should thus also depend on the β_s parameter of the inflow.

3.2.2. Sharpness of the inflow boundaries

We see from Fig. 4 that the boundary of the diffusion region, defined by the increase in particle number density, is very sharp in some cases (especially for the case ω_ce/ω_pe = 3 with a cold background plasma). These sharp transitions are not the trace of a shock between the incoming plasma and the over-dense diffusion region. First, because the inflow bulk velocity is not supersonic, in the sense that we have , with C_fms the phase speed of waves propagating perpendicular to the magnetic field B₀, i.e., the fast magnetosonic velocity in MHD (Alfvén and slow waves do not propagate perpendicularly to B₀). Second, because the width of the transition between inflow plasma and diffusion region plasma is, in the sharpest case, less than an electron inertial length, while we know that the thermalization of a cold inflow by collisionless kinetic processes occurs on a width of some inertial length, also with the formation of a precursor upstream (e.g., Plotnikov et al. 2013) that is not seen here. Third, because there is no variations in the magnetic field across the sharp particle density variation.

Instead, we explain these sharp transitions by quasi-ballistic motions of the particles entering the diffusion region. Figure 6 shows some particle trajectories: they travel toward the diffusion region because of the electric field E with a guiding center velocity given by v_{E × B} = E × B/B², then they reach the magnetic field reversal, overshoot to the other side because of their inertia, turn around the increasingly strong magnetic field, and then oscillate around the B = 0 line before being ejected toward the islands (because they have a high v_y in the B_x field, see Sect. 3.1 and Fig. 3 for a trajectory in 3D). This explains the fluid point of view of Fig. 4 well: averaging over particles to obtain the fluid particle number density or fluid velocity, the sharpness of the transitions comes from all particles of the same species turning around at roughly the same location.

It also explains the M shape of the particle densities of Fig. 4. Particles spend more time at the edge of the diffusion region, when turning back, than at the center, hence the peaks at the edges and the depletion at the center. The fluid velocity profiles are also well explained by the particle view. We note that an M shape is also reported in the context of the inversion electric field layer (Chen et al. 2011, and references therein). Here we do not find any trace of the inversion layer in the electric field. A possible explanation may be that our electrons are relativistic.

Concerning the sharp transitions, the question is thus to know why and under which conditions all the particles of a species turn back at the same location. They will do so if they all enter the diffusion region with the same velocity, i.e., if their thermal velocity is lower than their bulk velocity: ${\mathit{v}}_{\mathrm{in}\mathit{,s}}^{\mathrm{th}}\mathrm{\ll }\mathit{v̅}{}_{\mathrm{in}\mathit{,s}}\mathrm{~}{\mathit{v}}_{\mathit{E}\mathrm{\times }\mathit{B}}$. Table 3 presents the ratio ${\mathit{v}}_{\mathrm{in}\mathit{,s}}^{\mathrm{th}}\mathit{/}\mathit{v̅}{}_{\mathrm{in}\mathit{,s}}$ for the different simulations of Fig. 4. For a given magnetization, here ω_ce/ω_pe = 3, we clearly see a correlation between a small ${\mathit{v}}_{\mathrm{in}\mathit{,s}}^{\mathrm{th}}\mathit{/}\mathit{v̅}{}_{\mathrm{in}\mathit{,s}}$ and a sharp transition. In particular, in the case with very hot electrons (${\mathit{v}}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{th}}\mathit{/}\mathit{v̅}{}_{\mathrm{th}\mathit{,}\mathrm{e}}\mathrm{=}\mathrm{1.8}$) and relatively cold ions (${\mathit{v}}_{\mathrm{in}\mathit{,}\mathrm{i}}^{\mathrm{th}}\mathit{/}\mathit{v̅}{}_{\mathrm{in}\mathit{,}\mathrm{i}}\mathrm{=}\mathrm{0.1}$), the electron particle number and velocities present smooth variations, while the same curves for ions do exhibit sharp transitions. This is in accordance with the explanation of sharp transitions by the cold nature, in terms of the ratio ${\mathit{v}}_{\mathrm{in}}^{\mathrm{th}}\mathit{/}{\mathit{v}}_{\mathit{E}\mathrm{\times }\mathit{B}}$, of the inflowing plasma.

The influence of the magnetization on the sharpness of the transitions is seen by comparing the two first simulations of Table 3, with ω_ce/ω_pe = 1 and 3 and same ratio ${\mathit{v}}_{\mathrm{in}\mathit{,s}}^{\mathrm{th}}\mathit{/}\mathit{v̅}{}_{\mathrm{in}\mathit{,s}}$ (plotted in the top part of Fig. 4). A smaller magnetization implies a smoother shape. We understand this as a consequence of particles turning back on a scale given by their Larmor radii in the magnetic field at the sheet entrance, which is smaller when ω_ce/ω_pe is higher.

As a final remark concerning resolution, we note that in the coldest cases, the thermal Larmor radius of the electrons is smaller than a cell length. The resolution of the Larmor radius on the grid is, however, not important for integrating particle trajectories in constant fields, because interpolation of field quantities to particle position then provides the same (constant) values, regardless of the cell size. What matters is the temporal resolution, ω_csΔt < 1 (see Melzani et al. 2013, Sect. 3.1, for a discussion). Here, we have ω_csΔt = 0.04 and 0.12 for ω_ce/ω_pe = 1 or 3. Also, a test run with a spatial and temporal resolution increased by a factor two leads to the same structures.

In summary, the sharp transitions are explained by collective bouncing motions allowed by the cold nature of the inflow: ${\mathit{v}}_{\mathrm{in}\mathit{,s}}^{\mathrm{th}}\mathrm{\ll }\mathit{v̅}{}_{\mathrm{in}\mathit{,s}}$. Such transitions may also occur in non-relativistic reconnection, but then the inflow speed is low and the background plasma should be very cold. They are likely to be more common in relativistic reconnection, where is larger. This may be why, to our knowledge, they have never been reported in other simulations.

3.3. The relativistic Ohm’s law

We explore the balance of terms in Ohm’s law. The literature concerning 2D non-relativistic reconnection largely shows that for antiparallel reconnection, the dominant term is thermal inertia either in ion-electron (Klimas et al. 2010; Shay et al. 2007; Fujimoto 2009) or pair (Bessho & Bhattacharjee 2005) plasmas, and this is the key element of various analytical models (e.g., Hesse et al. 2011). On the other hand, reconnection with a guide field is sustained by electron bulk inertia on skin-depth scales, and thermal inertia on Larmor radius scales (Hesse et al. 2002, 2004; Liu et al. 2014). Existing studies with relativistic pair plasmas confirm the non-relativistic trend: with no guide field, Hesse & Zenitani (2007), Bessho & Bhattacharjee (2012), and Zenitani et al. (2009a) find a dominating contribution of thermal inertia, while with a guide field Hesse & Zenitani (2007) find a significant contribution of the time derivative of bulk momentum . Here we explore the situation in ion-electron relativistic plasmas for antiparallel reconnection. We show for the first time in the antiparallel case that bulk inertia contributes at least as much as thermal inertia for high inflow magnetizations. We explain that this is due to a relativistic inflow magnetization in Sect. 3.8.

Ohm’s law is the fluid equation of motion for the lighter particles (the electrons) and is a means to evaluating the relative weight of the different non-ideal terms, allowing reconnection with a simple fluid picture. Understanding the weight and localization of each term is a first step toward building an effective resistivity for fluid models (resistive MHD, two-fluid codes, hybrid codes), where concrete resistive parametrizations can lead to very different behaviors, for instance changing from steady to unsteady states in Zenitani et al. (2009b) or triggering (or not) a Petsheck-like configuration depending on the localization and gradients of the resistivity (Baty et al. 2006).

The relativistic Ohm’s law for electrons is directly derived from the equation of conservation of momentum for the electron fluid, Eq. (B.4), which is itself obtained from Vlasov equation in Appendix B. It reads as $\begin{array}{ccc}\frac{{\mathit{q}}_{\mathrm{e}}}{{\mathit{m}}_{\mathrm{e}}}\mathrm{(}{E}\mathrm{+}\overline{{v}}\mathrm{e}\mathrm{\wedge }{B}\mathrm{)}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{n}}_{\mathrm{e}}}\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{\left(}{\mathit{n}}_{\mathrm{e}}\overline{{p}}\mathrm{e}\mathrm{\right)}\mathrm{+}\frac{\mathrm{1}}{{\mathit{n}}_{\mathrm{e}}}\frac{\mathit{\partial }}{\mathit{\partial }{x}}\mathrm{·}\mathrm{(}{\mathit{n}}_{\mathrm{e}}\mathrm{⟨}{{p}}_{\mathrm{e}}{{v}}_{\mathrm{e}}\mathrm{⟩}\mathrm{)}\\ & \mathrm{=}& \begin{array}{c}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\frac{\mathrm{1}}{{\mathit{n}}_{\mathrm{e}}}\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{\left(}{\mathit{n}}_{\mathrm{e}}\overline{{p}}\mathrm{e}\mathrm{\right)}\mathrm{+}\frac{\mathrm{1}}{{\mathit{n}}_{\mathrm{e}}}\frac{\mathit{\partial }}{\mathit{\partial }{x}}\mathrm{·}\mathrm{\left(}{\mathit{n}}_{\mathrm{e}}\overline{{p}}\mathrm{e}\overline{{v}}\mathrm{e}\mathrm{\right)}}\\ \mathrm{bulkinertia}\end{array}\\ & & \end{array}$(8)Here we use for simplicity quantities computed in the simulation (or lab) frame; e.g., n_e is the electron number density in the labframe (denoted by n_lab,e in Appendix B). Brackets ⟨ ·⟩ _s or an overbar denote an average in momentum space over the particles distribution function. We also used the definition $\mathit{\delta }{p}\mathrm{=}{p}\mathrm{-}\overline{{p}}$ where p = γv is the momentum. The left hand side of Eq. (8) is the ideal part and is set equal to 0 in ideal MHD. On the right hand side are terms linked to finite inertia (i.e., particles do not respond perfectly to the electric field variations because of their inertia):

• The first term is a part of bulk inertia. However, it vanishes in steady state so we neglect it (validated a posteriori).

• The second term is inertia linked to the bulk flow and is denoted as bulk inertia. It comes from the overall structure of the flow around the sheet (the profiles of the mean quantities: increase in , ).

• The third term is inertia linked to microscopic thermal motion and is denoted as thermal inertia. It comes from the divergence of the off-diagonal terms of the pressure tensor and can be anticipated by a study of the temperature curves.

We analyze Ohm’s law in the direction of the reconnection electric field ($\hat{{y}}$ here). Given the invariance along y, the divergence terms have two contributions: ∑ _k∂_k(p_kv_y) = ∂_x(p_xv_y) + ∂_z(p_zv_y). A computation of the divergence requires a smoothing of the quantities, especially for the thermal inertia term, which is very noisy.

We show in Fig. 7 the results for a cut through the X point for the simulation with ω_ce/ω_pe = 3 and T_bg = 1.5 × 10⁷ K (${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{hot}}\mathrm{=}\mathrm{89}$). Ohm’s law is satisfied everywhere, except at the sharp transitions at the entrance of the electron diffusion region, where the derivatives diverge. Different areas emerge:

• The electrons are ideal outside of the ion diffusion region.

• In the ion diffusion region, $\mathrm{|}\overline{{v}}{}_{\mathrm{e}}\mathrm{\wedge }{B}\mathrm{|}$ decreases linearly. The bulk inertia term rises linearly to compensate. The term wins out over . The contribution of the former is understandable when looking at and , which increase when we get closer to the sheet in this region (see Fig. 4 for ). The thermal inertia term ∑ _k∂_k(n_eδp_kδv_y) is slightly positive, and partly cancels the contribution of bulk inertia. This cancellation is also reported in Fujimoto (2009) and Klimas et al. (2010) for non-relativistic ion-electron plasmas, and in Bessho & Bhattacharjee (2012) for relativistic pairs. Only ∂_x(n_eδp_xδv_y) contributes and is negative, which is easily seen when looking at the temperature curves T_xy,e and T_zy,e (Fig. 8).

• In the electron diffusion region, the ${\overline{{v}}}_{\mathrm{e}}\mathrm{\wedge }{B}$ term vanishes (because B is very weak and ). The bulk inertia term is constant and due only to the term , which is expected to contribute, given that and in this region (Fig. 9). The other term, , vanishes because in this area. The thermal inertia term ∑ _k∂_k(n_eδp_kδv_y) contributes as much as the bulk inertia term. Only ∂_x(n_eδp_xδv_y) contributes, and it is negative, which is easily seen when looking at the temperature curves T_xy,e and T_zy,e in Fig. 8.

A cut along x through other X points in the simulation leads to the same results. Also, a cut along z through the X point reveals that the results of the electron region hold throughout the center of the current sheet, with a slow increase in the ${\overline{{v}}}_{\mathrm{e}}\mathrm{\wedge }{B}$ term as we get near the islands.

In summary, non-ideal terms in the ion regions are due to bulk inertia and, in the electron diffusion region, to an interestingly equal contribution of bulk and thermal inertia. For other runs with ω_ce/ω_pe = 3 (${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{hot}}\mathrm{=}\mathrm{27}$ to 89), we also find an equal contribution from thermal and bulk inertia. For the most magnetized run, ω_ce/ω_pe = 6 (${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{hot}}\mathrm{=}\mathrm{260}$), the contribution of bulk inertia exceeds that of thermal inertia by a factor 1.5 to 3. We show in Sect. 3.8 with analytical estimations that the high magnetization for electrons indeed allows bulk inertia to overreach thermal inertia, with the former scaling as ${\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}$ and the latter as ${{}^{\mathrm{\right(}}{\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}^{\mathrm{\left)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}$. This effect is present in our simulations and not in the references previously mentioned with antiparallel fields, because our background electron magnetization is greater. It is thus a new regime that challenges the thermal inertia paradigm at high electron magnetizations. We discuss the possible consequences in Sect. 5.2.

3.4. Outflow: energy content of the exhaust jets

It can easily be shown (Sect. 3.4.1) from analytical considerations that the outflows from the diffusion region should have relativistic bulk velocities and/or relativistic temperatures. In our simulation data, the thermal part always clearly dominates the bulk kinetic energy part, and more strongly for more relativistic cases (Sect. 3.4.2). A refined analytical estimate explains why in Sect. 3.8.

3.4.1. A simple analytical estimation

As explained in Sect. 3.1, bipolar outflow jets are naturally produced from each side of the X points. They are clearly visible in Fig. 2. An estimation of the energy content of these outflows can be easily obtained in steady state by using the conservation of particle number and of energy. To do so, we consider that the diffusion region for species s has a length D_s (along z) and a width δ_s (along x). We generalize the situation to cases where there is a guide field B_G. We denote quantities entering (leaving) this region by a subscript “in” (“out”, respectively).

Conservation of particle number (Eq. (B.3)) gives n_in,sv_in,sD_s = n_out,sv_out,sδ_s. Regarding energy conservation (Eq. (B.9)), the inflow energy flux includes the particle energies, and the reconnecting and guide field Poynting fluxes (Eqs. (1) and (2)): ${\mathit{D}}_{\mathit{s}}{\mathit{n}}_{\mathrm{in}\mathit{,s}}{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}\mathit{p̅}{}_{\mathrm{in}\mathit{,s}}\mathrm{+}{\mathit{D}}_{\mathit{s}}{\mathit{v}}_{\mathrm{in}\mathit{,s}}\mathrm{·}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}{\mathit{\mu }}_{\mathrm{0}}\mathrm{+}{\mathit{D}}_{\mathit{s}}{\mathit{v}}_{\mathrm{in}\mathit{,s}}\mathrm{·}{\mathit{B}}_{\mathrm{G}\mathit{,}\mathrm{in}}^{\mathrm{2}}\mathit{/}{\mathit{\mu }}_{\mathrm{0}}$. We assume that in the outflow the energy in the reconnected magnetic field B₀ is negligible compared to particle energy, so that the energy flux is ${\mathit{\delta }}_{\mathit{s}}{\mathit{n}}_{\mathrm{out}\mathit{,s}}{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}\mathit{p̅}{}_{\mathrm{out}\mathit{,s}}\mathrm{+}{\mathit{\delta }}_{\mathit{s}}{\mathit{v}}_{\mathrm{out}\mathit{,s}}\mathrm{·}{\mathit{B}}_{\mathrm{G}\mathit{,}\mathrm{out}}^{\mathrm{2}}\mathit{/}{\mathit{\mu }}_{\mathrm{0}}$. Equating the two fluxes and combining this with the conservation of particle number, we obtain $\frac{\mathit{p̅}\mathrm{out}\mathit{,s}}{\mathit{v̅}\mathrm{out}\mathit{,s}}\mathrm{=}\frac{\mathit{p̅}\mathrm{in}\mathit{,s}}{\mathit{v̅}\mathrm{in}\mathit{,s}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{in}\mathit{,s}}^{\mathrm{cold}}\mathrm{\left(}{\mathit{B}}_{\mathrm{0}}\mathrm{\right)}\mathrm{+}{\mathit{\sigma }}_{\mathrm{in}\mathit{,s}}^{\mathrm{cold}}\mathrm{\left(}{\mathit{B}}_{\mathrm{G}\mathit{,}\mathrm{in}}\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{)}\mathit{,}$(9)with $\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{)}\mathrm{=}\left(\mathrm{1}\mathrm{-}\frac{{\mathit{n}}_{\mathrm{in}\mathit{,s}}{\mathit{B}}_{\mathrm{G}\mathit{,}\mathrm{out}}^{\mathrm{2}}}{{\mathit{n}}_{\mathrm{out}\mathit{,s}}{\mathit{B}}_{\mathrm{G}\mathit{,}\mathrm{in}}^{\mathrm{2}}}\right)$. The guide field is usually merely compressed, so that 1 − α ≃ 0.

Equation (9) is independent of the p dependence of the distribution function f_s. However, some insight can be gained by considering a distribution that is isotropic in the comoving frame, for which we have the result $\overline{{p}}\mathrm{=}{\mathit{h}}_{\mathrm{0}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{\Gamma }\overline{{v}}$, with h₀ the comoving enthalpy (as defined and pictured for a Maxwell-Jüttner distribution in Fig. B.1), and $\mathrm{\Gamma }\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\overline{{v}}{}^{\mathrm{2}}\mathit{/}{\mathit{c}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$. If, in addition, we neglect the contribution of the guide field and assume an inflow plasma with non-relativistic temperatures and non-relativistic velocities, Eq. (9) becomes ²

${\mathit{h}}_{\mathrm{0}\mathit{,}\mathrm{out}\mathit{,s}}{\mathrm{\Gamma }}_{\mathrm{out}\mathit{,s}}\mathrm{=}\mathrm{1}\mathrm{+}{\mathit{\sigma }}_{\mathrm{in}\mathit{,s}}^{\mathrm{cold}}\mathrm{\left(}{\mathit{B}}_{\mathrm{0}}\mathrm{\right)}\mathit{.}$(11)We underline that the magnetization ${\mathit{\sigma }}_{\mathrm{in}\mathit{,s}}^{\mathrm{cold}}\mathrm{\left(}{\mathit{B}}_{\mathrm{0}}\mathrm{\right)}$ is to be taken at the entrance of the diffusion region of species s, where it can differ from its asymptotic value because of a decrease in magnetic field and particle number density (as in Fig. 4).

From Eq. (11), we clearly see that for a relativistic inflow plasma where B²/μ₀>nmc² and hence ${\mathit{\sigma }}_{\mathrm{in}\mathit{,s}}^{\mathrm{cold}}\mathrm{\left(}{\mathit{B}}_{\mathrm{0}}\mathrm{\right)}\mathit{>}\mathrm{1}$, magnetic reconnection is expected to produce outflows with either relativistic bulk velocities (Γ_out,s> 1) or relativistic temperatures (h_0,out,s> 1), or both. We also see that since ${\mathit{\sigma }}_{\mathit{s}}^{\mathrm{cold}}\mathrm{\propto }\mathrm{1}\mathit{/}{\mathit{m}}_{\mathit{s}}$, electrons will be more accelerated/heated than ions and that relativistic electrons (${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{cold}}\mathit{>}\mathrm{1}$) can be expected even at low ion magnetizations (${\mathit{\sigma }}_{\mathrm{i}}^{\mathrm{cold}}\mathrm{=}{\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{cold}}\mathrm{\times }{\mathit{m}}_{\mathrm{e}}\mathit{/}{\mathit{m}}_{\mathrm{i}}\mathrm{\ll }\mathrm{1}$).

3.4.2. Results from simulations

We first check whether the energy estimate of Eqs. (9) and (11) holds in all simulations. This is the case up to a factor ≲ 2. An only approximate correspondence is to be expected because this relation assumes a simple geometry with, in particular, a constant inflow velocity along the diffusion region and no energy exchange between the species. For example, in Fig. 4, for ω_ce/ω_pe = 3 and T_bg = 1.5 × 10⁷ K, we measure $\mathrm{1}\mathrm{+}{\mathit{\sigma }}_{\mathrm{in}\mathit{,s}}^{\mathrm{cold}}\mathrm{\left(}{\mathit{B}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{1.8}$ in the inflow for ions and 6.2 for electrons, while at the outflow maximal velocity (Fig. 9) we have for ions and 13 for electrons. These orders of magnitude hold for all cases.

A more refined analysis of the energy content of the outflow, split into its thermal and bulk contributions, can be performed. To do so, we decompose the particle energy flux as (see Appendix B) $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{n}}_{\mathit{s}}\mathrm{⟨}\mathit{\gamma }{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}{v}{\mathrm{⟩}}_{\mathit{s}}\mathrm{=}{\mathit{n}}_{\mathit{s}}{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}{\mathit{h}}_{\mathrm{0}\mathit{,}\mathrm{out}\mathit{,s}}{\mathrm{\Gamma }}_{\mathrm{out}\mathit{,s}}\overline{{v}}\mathrm{out}\mathit{,s}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{=}{\mathit{n}}_{\mathit{s}}{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}\overline{{v}}\mathrm{out}\mathit{,s}\left[\mathrm{1}\mathrm{+}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{out}\mathit{,s}}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{+}{\mathrm{\Gamma }}_{\mathrm{out}\mathit{,s}}\mathrm{(}{\mathit{h}}_{\mathrm{0}\mathit{,}\mathrm{out}\mathit{,s}}\mathrm{-}\mathrm{1}\mathrm{)}\right]\mathit{.}\end{array}$(12)On the right hand side, the first term is the restmass energy flux and is the same as from the inflow. The second is the kinetic energy of a cold bulk flow of velocity $\overline{{v}}{}_{\mathrm{out}\mathit{,s}}$. The third is the energy transported by thermal motions in the plasma restframe and by pressure work, and we denote it as the enthalpy flux. We note that these definitions match those of Zenitani et al. (2009a), who performed a similar analysis with two-fluid simulations of relativistic reconnection in pair plasmas.

We measure the maximum outflow velocity , deduce the Lorentz factor Γ_out,s, measure the maximum in momentum , and compute the enthalpy . From these values, we estimate in Table 4 the balance of particle energy between each of the terms of Eq. (12).

In all cases, a large fraction of the particle energy flux is in thermal kinetic energy, not in bulk flow kinetic energy. For electrons, we see that the thermal part clearly dominates more as one increases the relativistic nature of the inflow (e.g., 69% in the thermal part for the less relativistic case, 99% for the most relativistic). This is also the case for ions: from 80% to 92% in the thermal part as the ion magnetization increases. The T_bg,i = 2 × 10⁸ K case is exceptional with 91% in the thermal energy, but this high percentage is likely explained by interactions with the hot electrons T_bg,e = 3 × 10⁹ K. We explain why thermally dominated outflows are expected at high inflow magnetization with a refined analytical model in Sect. 3.8.

3.5. Islands structure

Turning to the magnetic islands, we emphasize that they are magnetically isolated, have an M-shaped density distribution, and are hot with anisotropic temperatures. After being expelled along ± z in the outflow, the particles meet the magnetic islands that separate each pair of X points. The islands are initially formed by the tearing of the current sheet. They then consist only of particles of the current sheet, plus those of the background plasma that were inside the current sheet location. Small at the beginning, they grow by collecting particles from the outflows at their periphery and by merging with other islands. A remarkable property is that particles from the background plasma cannot enter the islands, but are scattered by the strong magnetic field structure surrounding the island and circle around it. Consequently, the particles at the island centers remain the same throughout the whole simulation, even after many island-merging events. This matter is explored in more detail in Melzani et al. (2014b). We stress here two main points.

First, the trapped particles are heated by the island contraction (which occurs when two islands merge), so that the central temperatures are very high for both species, highly anisotropic (Figs. 9 and 10), with a dominance of T_zz. The island centers are also where the currents are the strongest.

Second, as said above, most of the inflowing background particles populate a region around the center, while the central part of the island mainly consists of particles originally from the current sheet. As a result, the surrounding part becomes denser than the isolated central part after some time. A cut along z through an island center indeed (Fig. 9) reveals, for the particle density, an M shape with a central dip and two shoulders. This may explain observed density dips at the center of magnetic islands during magnetotail reconnection events (Khotyaintsev et al. 2010), without invoking island merging or particle escape along the flux tube.

3.6. Reconnection electric field and reconnection rate

The rate of variation in magnetic field flux across an X point, $\mathrm{d}{\mathrm{\Phi }}_{{\mathit{B}}_{\mathit{z}}}\mathit{/}\mathrm{d}\mathit{t}\mathrm{=}\mathrm{(}\mathrm{d}\mathit{/}\mathrm{d}\mathit{t}\mathrm{)}{{}^{\mathrm{\int }}}_{\mathit{x}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{0}}^{\mathrm{Xpoint}}{\mathit{B}}_{\mathit{z}}\mathrm{d}\mathit{x}$, is equal in 2D configurations to the y component E_y of the electric field at the X-point location. In addition, dΦ_Bz/ dt is in part determined by the outflow velocity, because the latter sets the rate at which magnetic field is extracted from around the X point (see, e.g., in a resistive MHD context, Borovsky & Hesse 2007; Cassak & Shay 2007). Since one expects in non-relativistic setups, the reconnection rate E_y is usually normalized either to ${\mathit{B}}_{\mathrm{0}}{\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{0}}^{\mathrm{NR}}$, with ${\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{0}}^{\mathrm{NR}}$ the hybrid Alfvén speed of Eq. (7), or to ${\mathit{B}}_{\mathrm{0}}{\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{in}}^{\mathrm{NR}}$, with ${\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{in}}^{\mathrm{NR}}$ the Alfvén speed in the inflow of Eq. (6a). These normalizations are chosen so that the normalized rate, E^∗ = E_y/B₀V_A, stays close to the same set of values. For example, it has been shown that it gives identical results when varying the mass ratio (e.g., Hesse et al. 1999; or Ricci et al. 2002, 2003, for m_i/m_e = 25,180,1836 with an implicit PIC code)³.

In the following we turn to our relativistic case and ask whether a normalization can be found that confines the range of values for E^∗ in a narrow range and relaxes to the above normalization in the non-relativistic case. We argue here that the normalization by the hybrid Alfvén speed is not relevant, because it does not depend on the particle number density of the inflow, while the ratio E_y/B₀ clearly does. This is seen for the simulation with n_bg = 0.3n_cs(0), for which E_y/B₀ peaks at 0.13c, compared to the otherwise identical simulation with n_bg = 0.1n_cs(0), where E_y/B₀ peaks at 0.20c. On the other hand, the inflow Alfvén speed ${\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{in}}^{\mathrm{NR}}\mathrm{\propto }\mathrm{1}\mathit{/}\sqrt{{\mathit{n}}_{\mathrm{bg}}}$, and thus leads to closer normalized rates. We consequently exclude hybrid quantities for normalization.

In a relativistic configuration, the non-relativistic Alfvén speed can increase to infinity. However, the ratio E_y/B₀ is also the E × B velocity of the incoming plasma and cannot exceed the speed of light. The normalizing Alfvén velocity should thus also saturate to some value, which is why we choose to normalize the electric field by ${\mathit{E}}^{\mathrm{\ast }}\mathrm{=}\frac{{\mathit{E}}_{\mathit{y}}}{{\mathit{B}}_{\mathrm{0}}{\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{in}}^{\mathrm{R}}}\mathit{,}$(13)with ${\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{in}}^{\mathrm{R}}$ the relativistic Alfvén speed in the inflow (Eq. (6b)), which cannot exceed c. The time evolution of E^∗ is shown in Fig. 11. Several comments can be made.

First, the rate E^∗ is not sensitive to the background plasma temperatures, as can be seen for the simulations ω_ce/ω_pe = 3, n_bg = 0.1n_cs(0), no guide field, and T_bg = 1.5 × 10⁷, 2 × 10⁸ and 3 × 10⁹ K. This contrasts with the interpretation of Hesse & Zenitani (2007), who attribute a lower rate to a higher inflow temperature. In addition to the temperatures, the magnetization of their simulation also changes and may also affect the rates. Coming back to our simulations, we note that we use very low background plasma β (< 10^-2, Table 2) and that a weak plasma β dependence is expected for higher values (e.g., TenBarge et al. 2014, have rates E^∗ multiplied by ~ 2 when β passes from 0.01 to 1).

Second, the reconnection rate for the simulation with a higher background particle density (n_bg = 0.3n_cs(0), E^∗ = 0.18) remains lower than its counterpart with n_bg = 0.1n_cs(0) (E^∗ = 0.23). This is in line with the pair plasma simulations of Bessho & Bhattacharjee (2012) who found a similar rate for n_bg = 0.1n_cs(0) (E^∗ = 0.19) and a higher rate for n_bg = 0.01n_cs(0) (E^∗ = 0.36). The reconnection rate thus increases with decreasing background plasma density, which is also coherent with the β dependence mentioned above.

Finally, the normalization leads to very similar values of E^∗ for the relativistic cases (ω_ce/ω_pe = 3 or 6), with E^∗ = 0.17−0.24, but to a significantly lower rate for the less relativistic case (ω_ce/ω_pe = 1), with E^∗ = 0.14. More generally, the values for the relativistic cases are higher than those reported in the literature for undriven, symmetric reconnection with zero guide field in non-relativistic ion-electron plasmas. For the peak values of E^∗ (once normalized in the same way as here) we can cite: Birn et al. (2001); Pritchett (2001): 0.09, Fujimoto (2006, 2009): 0.15, Daughton et al. (2006): 0.08, Klimas et al. (2010): 0.07−0.09, and the theoretical work of Hesse et al. (2009a,b) predicting a maximum rate of 0.28. Our results thus suggest higher rates for relativistic reconnection, which is already seen in relativistic simulations of pair plasmas with, e.g., Zenitani & Hoshino (2007) (E^∗ = 0.2), Cerutti et al. (2012b) (E^∗ = 0.17), or Bessho & Bhattacharjee (2012) (E^∗ = 0.19 and 0.36).

In conclusion, the relativistic Alfvén speed of the inflow provides the best normalization for the reconnection electric field, in that it is robust from non-relativistic to relativistic flows. Corresponding rates are in a close range (0.14–0.25), which is higher than the rates found in non-relativistic simulations with the same normalization (0.07–0.15). The rate does not depend on the inflow temperature at low β, but is nevertheless not universal: it decreases with increasing background particle number density. Generalization to the presence of a guide field is discussed in Sect. 4.3.

3.7. Hall field and dispersive waves

We can see in Fig. 4 that inside the ion-diffusion region, but outside of the electron-diffusion region, ions have a low fluid velocity, while electrons still E × B drift toward their diffusion region. This results in a net current roughly given by ${\mathit{q}}_{\mathrm{e}}{\mathit{n}}_{\mathrm{e}}\overline{{v}}{}_{\mathrm{e}}\mathrm{=}{\mathit{q}}_{\mathrm{e}}{\mathit{n}}_{\mathrm{e}}{E}\mathrm{\wedge }{B}\mathit{/}{\mathit{B}}^{\mathrm{2}}$, which is the Hall current. This current continues along the magnetic separatrices in the outflow direction, and is at the origin of a quadripolar magnetic field directed along$\mathrm{\pm }\hat{{y}}$. This Hall magnetic field, with the quadripolar structure, is present in our simulations. It has a weak intensity (between 1% and 10% of B₀). The charge separation between electrons and ions (Fig. 4) also leads to the creation of a Hall electric field directed along $\mathrm{+}\hat{{x}}$ in the x < 0 region and $\mathrm{-}\hat{{x}}$ in the x < 0 region. Both the magnetic and electric Hall fields are absent in a simulation with pairs.

The difference in the dynamical response of ions and electrons also allows the existence of waves with a quadratic dispersion relation, ω ∝ k², below ion scales (either whistler waves or kinetic Alfvén waves, see Rogers et al. 2001). Observations of the same reconnection rate for any simulation model allowing these waves (PIC, electron-MHD, Hall-MHD, two-fluid with and without electron inertia, hybrid simulations, see Birn et al. 2001; Shay et al. 2007; Rogers et al. 2001), as well as theoretical considerations, have led to the thesis that these waves are essential to allow for fast reconnection rates. However, this view is questioned by a number of simulations that do not support quadratic dispersive waves, but still support fast rates (hybrid simulations with no Hall term, pair plasmas, or strong guide field regime, see Karimabadi et al. 2004; Bessho & Bhattacharjee 2005; Daughton et al. 2006; Daughton & Karimabadi 2007; Liu et al. 2014). It is thus interesting to see whether our simulation data can provide any further insight into the matter.

A prediction of the dispersive wave physics is that the reconnection rate is controlled solely by the ion physics and not by the electrons. According to Daughton et al. (2006), it should be independent of the electron diffusion region length. We could, however, not reproduce their analysis because the electron diffusion zone length is, in our case, limited by the standing islands. It cannot stretch to high values, so we are unable to conclude for or against of the dispersive wave paradigm.

However, we underline that the simulation with m_i/m_e = 1 that we performed features an identical (and even slightly greater, Fig. 11) reconnection rate than simulations with m_i/m_e = 25. It thus points toward a negligible influence of the dispersive waves or to another mechanism that allows fast rates in pair plasmas.

3.8. Simulation-based scaling analysis

The energy content of the outflows and the balance between thermal and bulk inertia in Ohm’s law were explored through the simulations in Sects. 3.3 and 3.4. The aim of the present section is to investigate these points with a simple analytical model in order to gain physical insight into these phenomena and to extrapolate our simulation results to a larger parameter space.

We extend the analytical results of Sect. 3.4.1, where particle number and energy conservation allowed an estimation of the quantity h_0,out,sΓ_out,s (Eqs. (9) or (11)) by now also using the equation of conservation of momentum (Eq. (B.4)).

3.8.1. Thermal versus bulk electron inertia

We first investigate the relative weight of thermal and bulk electron inertia. At the center of the electron diffusion region, we learn from Sect. 3.3 that the reconnection electric field is sustained by electron thermal and bulk inertia, with only the terms ∂_x(n_e ⟨ δp_xδv_y ⟩) and contributing to either one of them, respectively.

• Concerning thermal inertia, the temperature tensor is defined via Eq. (B.10), so that ⟨ δp_xδv_y ⟩ = c²Θ_xy,e/ Γ_e. We see in Fig. 8 that Θ_xy,e is linear in the electron diffusion region. It vanishes at the center because there the distribution function f_e is symmetric with respect to v_x. It is maximal at the diffusion region edge with a value ${\mathrm{\Theta }}_{\mathit{xy,}\mathrm{e}}^{\mathrm{edge}}$. Consequently, we approximate the thermal inertia contribution by $\mathrm{(}{\mathit{c}}^{\mathrm{2}}{\mathrm{\Theta }}_{\mathit{xy,}\mathrm{e}}^{\mathrm{edge}}\mathit{/}{\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{in}}\mathrm{)}\mathit{/}{\mathit{\delta }}_{\mathrm{e}}$, where δ_e is the width of the electron diffusion region.

• For the bulk inertia term, we use the fact that rises linearly from the center to its maximum value denoted by $\mathit{p̅}\begin{array}{c}\mathrm{out}\\ \mathit{z}\end{array}$ over a distance D_e/ 2 and that has a vanishing derivative at the center (Fig. 9). Consequently, it can be estimated as $\mathit{v̅}\begin{array}{c}\mathrm{center}\\ \mathit{y}\end{array}\mathit{p̅}\begin{array}{c}\mathrm{out}\\ \mathit{z}\end{array}\mathit{/}{\mathit{D}}_{\mathrm{e}}$.

All in all, from Ohm’s law (Eq. (8)), the electric field at the center of the diffusion region is $\begin{array}{ccc}{\mathit{E}}_{\mathit{y}}^{\mathrm{center}}& \mathrm{=}& \frac{{\mathit{m}}_{\mathrm{e}}}{{\mathit{q}}_{\mathrm{e}}{\mathit{n}}_{\mathrm{e}}}{\left(\frac{\mathit{\partial }}{\mathit{\partial }{x}}\mathrm{·}\mathrm{(}{\mathit{n}}_{\mathrm{e}}\mathrm{⟨}\mathit{\delta }{{p}}_{\mathrm{e}}\mathit{\delta }{{v}}_{\mathrm{e}}\mathrm{⟩}\mathrm{)}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial }{x}}\mathrm{·}\mathrm{\left(}{\mathit{n}}_{\mathrm{e}}\overline{{p}}\mathrm{e}\overline{{v}}\mathrm{e}\mathrm{\right)}\right)}_{\mathit{y}}\\ & \mathrm{~}& \end{array}$(14)The next step is to use the constancy of E_y, which is respected well in the simulations: ${\mathit{E}}_{\mathit{y}}^{\mathrm{center}}\mathrm{=}{\mathit{E}}_{\mathit{y}}^{\mathrm{in}}\mathrm{=}\mathit{v̅}{}_{\mathrm{in}}\mathit{B}_{\mathrm{0}}$. We note, however, that at the entrance of the diffusion regions B is different than the asymptotic value B₀ and that the E × B drift does not strictly hold (see Fig. 4). If we introduce the inertial length in the inflow, ${\mathit{d}}_{\mathrm{e}}^{\mathrm{in}}\mathrm{=}\mathit{c}\sqrt{{\mathit{ϵ}}_{\mathrm{0}}{\mathit{m}}_{\mathrm{e}}\mathit{/}\mathrm{\left(}{\mathit{n}}_{\mathrm{e}}^{\mathrm{in}}{\mathit{e}}^{\mathrm{2}}\mathrm{\right)}}$, and the inflow magnetization ${\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}\mathrm{=}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}{\mathit{\mu }}_{\mathrm{0}}{\mathit{m}}_{\mathrm{e}}{\mathit{n}}_{\mathrm{e}}^{\mathrm{in}}{\mathit{c}}^{\mathrm{2}}\mathrm{\right)}$, we ultimately obtain $\frac{{\mathit{\delta }}_{\mathrm{e}}}{{\mathit{d}}_{\mathrm{e}}^{\mathrm{in}}}{\mathrm{(}{{\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}\frac{\mathit{v̅}\mathrm{in}}{\mathit{c}}\mathrm{=}\frac{{\mathrm{\Theta }}_{\mathit{xy,}\mathrm{e}}^{\mathrm{edge}}}{{\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{in}}}\mathrm{+}\frac{{\mathit{\delta }}_{\mathrm{e}}}{{\mathit{D}}_{\mathrm{e}}}\frac{\mathit{v̅}\begin{array}{c}\mathrm{center}\\ \mathit{y}\end{array}\mathit{p̅}\begin{array}{c}\mathrm{out}\\ \mathit{z}\end{array}}{{\mathit{c}}^{\mathrm{2}}}\mathrm{·}$(15)We now proceed to derive approximate scaling relations for cases where either thermal or bulk inertia dominate the reconnection electric field.

• First, if thermal inertia dominates over bulk inertia, then Eq. (15) gives $\frac{{\mathrm{\Theta }}_{\mathit{xy,}\mathrm{e}}^{\mathrm{edge}}}{{\mathrm{\Gamma }}_{\mathrm{e}}^{\mathrm{in}}}\mathrm{=}\frac{{\mathit{\delta }}_{\mathrm{e}}}{{\mathit{d}}_{\mathrm{e}}^{\mathrm{in}}}{\mathrm{(}{{\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}\frac{\mathit{v̅}\mathrm{in}\mathit{,}\mathrm{e}}{\mathit{c}}\mathrm{\propto }{\mathrm{(}{{\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(16)There are thus several factors contributing to ${\mathrm{\Theta }}_{\mathit{xy,}\mathrm{e}}^{\mathrm{edge}}$. The diffusion zone width δ_e is dynamically set during the reconnection process. It can be close to the particles gyroradius at the center of the current sheet or to the plasma inertial length at the center of the current sheet. In all our simulations, we find that the latter assumption holds throughout time to within a factor 2 (Sect. 3.2.1), and in any case, ${\mathit{\delta }}_{\mathrm{e}}\mathit{/}{\mathit{d}}_{\mathrm{e}}^{\mathrm{in}}$ is expected to be of order unity. The inflow speed is set by the reconnection electric field, $\mathit{v̅}{}_{\mathrm{in}}\mathrm{=}{\mathit{E}}_{\mathit{y}}\mathit{/}{\mathit{B}}_{\mathrm{0}}\mathrm{=}{\mathit{E}}^{\mathrm{\ast }}{\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{in}}^{\mathrm{R}}$ with E^∗ the normalized reconnection rate (which lies in the range 0.1−0.25, Sect. 3.6) and ${\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{in}}^{\mathrm{R}}$ the relativistic Alfvén speed in the inflow. For relativistic setups we thus have . The inflow magnetization can be arbitrarily large. It is thus the main actor in producing relativistic temperatures, and thermal inertia scales as ${\mathrm{\Theta }}_{\mathit{xy,}\mathrm{e}}^{\mathrm{edge}}\mathrm{\propto }{{}^{\mathrm{\right(}}{\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}^{\mathrm{\left)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}$.

• Second, the term corresponding to bulk inertia in Eq. (15) can be estimated with the help of Eq. (9) (with $\mathit{p̅}{}_{\mathrm{out}\mathit{,s}}\mathrm{=}\mathit{p̅}\begin{array}{c}\mathrm{out}\\ \mathit{z}\end{array}$, $\mathit{v̅}{}_{\mathrm{out}\mathit{,s}}\mathrm{=}\mathit{v̅}\begin{array}{c}\mathrm{out}\\ \mathit{z}\end{array}$, and neglecting the guide field): $\frac{{\mathit{\delta }}_{\mathrm{e}}}{{\mathit{D}}_{\mathrm{e}}}\frac{\mathit{v̅}\begin{array}{c}\mathrm{center}\\ \mathit{y}\end{array}\mathit{p̅}\begin{array}{c}\mathrm{out}\\ \mathit{z}\end{array}}{{\mathit{c}}^{\mathrm{2}}}\mathrm{=}\frac{{\mathit{\delta }}_{\mathrm{e}}}{{\mathit{D}}_{\mathrm{e}}}\frac{\mathit{v̅}\begin{array}{c}\mathrm{center}\\ \mathit{y}\end{array}\mathit{v̅}\begin{array}{c}\mathrm{out}\\ \mathit{z}\end{array}}{{\mathit{c}}^{\mathrm{2}}}\left(\frac{\mathit{p̅}\mathrm{in}}{\mathit{v̅}\mathrm{in}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}\right)\mathit{.}$(17)The ratio δ_e/D_e is of order 1 / 10 in our simulations. If we neglect the term , which is of order unity for non-relativistic inflow temperatures, we see that bulk inertia scales with ${\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}$.

In conclusion, thermal inertia scales at most as ${{}^{\mathrm{\right(}}{\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}^{\mathrm{\left)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}$, and bulk inertia as ${\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}$. Consequently, regarding the non-ideal terms in Ohm’s law in the electron diffusion region, we expect bulk inertia to outweight thermal inertia at high inflow electron magnetization.

3.8.2. Energy content of the outflows

We now turn to the energy content of the outflows to see whether we can explain their thermally dominated character for relativistic runs. The temperature in the outflows is dominated by Θ_xx,e or Θ_yy,e, which we denote by ${\mathrm{\Theta }}_{\mathrm{e}}^{\mathrm{out}}$. We first have to link ${\mathrm{\Theta }}_{\mathrm{e}}^{\mathrm{out}}$ to ${\mathrm{\Theta }}_{\mathit{xy,}\mathrm{e}}^{\mathrm{edge}}$. The outflow temperature at the center of the diffusion region is roughly constant along z throughout the area of linear increase of (Fig. 9), because particles on their way from the X point to the exhaust mainly turn into the reconnected magnetic field and thus do not really gain thermal agitation, but convert it from one component of Θ to another. We can thus assume ${\mathrm{\Theta }}_{\mathrm{e}}^{\mathrm{out}}\mathrm{=}{\mathrm{\Theta }}_{\mathrm{e}}^{\mathrm{center}}$. We now would like to assume ${\mathrm{\Theta }}_{\mathit{xy,}\mathrm{e}}^{\mathrm{edge}}\mathrm{~}{\mathrm{\Theta }}_{\mathit{xx,}\mathrm{e}}^{\mathrm{center}}$. This indeed holds for electrons in the case of Fig. 8. However, this does not hold in all simulations, and ${\mathrm{\Theta }}_{\mathit{xy,}\mathrm{e}}^{\mathrm{edge}}$ is between 1 / 10 to 10 times ${\mathrm{\Theta }}_{\mathit{xx,}\mathrm{e}}^{\mathrm{center}}$. This is due to the different origin of these components: ${\mathrm{\Theta }}_{\mathit{xx,}\mathrm{e}}^{\mathrm{center}}$ reflects particles in Speiser orbits going up and down along x with a zero bulk x velocity, while ${\mathrm{\Theta }}_{\mathit{xy,}\mathrm{e}}^{\mathrm{edge}}$ reflects the asymmetry of the distribution function with respect to v_x due to the newly entering particles at the edge of the diffusion zone.

With the previous remark in mind, we still make the hypothesis ${\mathrm{\Theta }}_{\mathit{xy,}\mathrm{e}}^{\mathrm{edge}}\mathrm{~}{\mathrm{\Theta }}_{\mathit{xx,}\mathrm{e}}^{\mathrm{center}}$. Next, if we assume that thermal inertia contributes significantly in Ohm’s law, we obtain, with the scaling of Eq. (16), ${\mathrm{\Theta }}_{\mathrm{e}}^{\mathrm{out}}\mathrm{\propto }{\mathrm{(}{{\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(18)For relativistic temperatures we have ${\mathit{h}}_{\mathrm{0}\mathit{,}\mathrm{out}\mathit{,}\mathrm{e}}\mathrm{\simeq }\mathrm{4}{\mathrm{\Theta }}_{\mathrm{e}}^{\mathrm{out}}$ (Fig. B.1), so that with Eq. (18) we see that a relativistic inflow magnetization implies ${\mathit{h}}_{\mathrm{0}\mathit{,}\mathrm{out}\mathit{,}\mathrm{e}}\mathrm{\propto }{{}^{\mathrm{\right(}}{\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}^{\mathrm{\left)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}$. On the other hand, energy conservation gives, in its simplest form (Eq. (11) with ${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{cold}}\mathrm{\left(}{\mathit{B}}_{\mathrm{0}}\mathrm{\right)}\mathrm{\gg }\mathrm{1}$): ${\mathit{h}}_{\mathrm{0}\mathit{,}\mathrm{out}\mathit{,}\mathrm{e}}{\mathrm{\Gamma }}_{\mathrm{out}\mathit{,}\mathrm{e}}\mathrm{~}{\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}\mathit{.}$(19)Thus, ${\mathrm{\Gamma }}_{\mathrm{out}\mathit{,}\mathrm{e}}\mathrm{\propto }{\mathrm{(}{{\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(20)We finally turn to the ratio of energy fluxes in the outflow. We see with Eq. (12) that the flux associated with kinetic bulk energy is Γ_out,e − 1. With Eq. (20) (and for Γ_out,e ≫ 1), this flux is ${\mathrm{\Gamma }}_{\mathrm{out}\mathit{,}\mathrm{e}}\mathrm{\propto }{{}^{\mathrm{\right(}}{\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}^{\mathrm{\left)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}$. The flux associated with thermal kinetic energy and pressure work is h_0,out,eΓ_out,e − 1, and with Eq. (19) we have ${\mathit{h}}_{\mathrm{0}\mathit{,}\mathrm{out}\mathit{,}\mathrm{e}}{\mathrm{\Gamma }}_{\mathrm{out}\mathit{,}\mathrm{e}}\mathrm{-}\mathrm{1}\mathrm{~}{\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}$. The ratio of thermal-to-bulk energy fluxes is thus $\mathrm{\propto }{{}^{\mathrm{\right(}}{\mathit{\sigma }}_{\mathrm{in}\mathit{,}\mathrm{e}}^{\mathrm{cold}}^{\mathrm{\left)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}$, meaning that relativistic inflow magnetization inevitably implies reconnection exhausts dominated by thermal energy, which is what we see in our simulations (Table 4), even if the scalings derived here do not hold exactly because of the many assumptions involved.

In conclusion, we have shown that under the hypothesis of non-ideal effects sustained by thermal inertia, relativistic inflow magnetizations produce thermally dominated outflows. The physical reason is that the reconnection electric field E_y is large in the inflow region, so that thermal inertia must be high in order to sustain E_y in the central region, which implies high temperatures.

However, we also demonstrated that thermal inertia is not expected to dominate for very relativistic inflows. When this is the case, there is no constraints from Ohm’s law on the temperature, and we cannot decide on the ratio of thermal to bulk energy fluxes. Since this ratio is (h_0,out,eΓ_out,e − 1) / (Γ_out,e − 1) ~ h_0,out,e, a relativistic outflow temperature on the order of m_ec² suffices to give thermally dominated outflows.

For our simulations, thermal inertia contributes equally or less than bulk inertia (Sect. 3.3), but still significantly, so that the outflows are thermally dominated.

4. Effects of a guide field

Except in special configurations, the generic reconnection geometry involves asymptotic fields that are not antiparallel. An angle different from 180° can be described by the addition of a uniform guide magnetic field ${{B}}_{\mathrm{G}}\mathrm{=}{\mathit{B}}_{\mathrm{G}}\hat{{y}}$ to the antiparallel configuration. Such configurations have been largely studied in the non-relativistic case (e.g., Pritchett & Coroniti 2004; Drake et al. 2005; Goldman et al. 2011; Le et al. 2013; Graf von der Pahlen & Tsiklauri 2014), and feature significant differences with the antiparallel case. Here we only focus on the reconnection rates and on the island structure, and postpone a study of other points to a future publication. We present results from two simulations, with B_G = 0.5B₀ and B_G = B₀.

4.1. Overall structure

We first stress that because of the presence of the guide field, in both cases B_G = 0.5B₀ and B₀, the relation B>E holds everywhere through time and space, hence also the relation B>E_perp (with E_perp the component perpendicular to B). Consequently, particle acceleration is not possible in directions perpendicular to B and is only allowed along the field lines at places where E·B ≠ 0. Such parallel electric fields are allowed by the non-idealness of the plasma response, ${E}\mathrm{+}\overline{{v}}{}_{\mathit{s}}\mathrm{\wedge }{B}\ne \mathrm{0}$, and are indeed found at and around the X points.

Just as in the zero guide field case, the plasma accelerated by E_y is slowly deviated by the reconnected B_x component, which produces outflows directed along ± z. However, particles from these outflows feel the Lorentz force from the guide field, and their trajectories are tilted against the z axis, as is evident in Fig. 12. Reversing the guide field from $\mathrm{+}\hat{{y}}$ to $\mathrm{-}\hat{{y}}$ tilts them in the opposite direction. We underline that while the tilt angle is indeed smaller for a smaller guide field, it also depends on the background plasma pressure, as shown by TenBarge et al. (2014).

In the present case, ${{E}}_{\mathit{y}}\mathrm{\propto }\mathrm{-}\hat{{y}}$. Electrons are thus accelerated along the field lines in the $\mathrm{+}\hat{{y}}$ direction. Their motion along the field lines results in a projected (x-z plane) motion directed toward positive z where ${B}\mathrm{·}\hat{{z}}\mathit{>}\mathrm{0}$ (i.e., in the x> 0 area), or toward negative z where ${B}\mathrm{·}\hat{{z}}\mathit{<}\mathrm{0}$ (i.e., in the x < 0 area). It results in large and fast electron flows above and below the current sheet, and to a rotation around the islands in a direction opposed to cyclotron gyration. Ions are accelerated toward $\mathrm{-}\hat{{y}}$ and counterstream with respect to the electrons.

Particles reaching the central part are accelerated along y by E_y and, because they always feel a magnetic field B>E, they are guided by the magnetic field and spend more time in the acceleration region for strong B_G (Fig. 13). Consequently, is greater than with no guide field under similar conditions, and reaches high values on a larger area (compare in Figs. 2 and 12).

4.2. Islands’ structure

The structure of the islands resembles the no guide field case: their centers are composed of particles initially in the current sheet, with background particles only circling at the periphery. They are the hottest and strongest current-carrying part of the simulations. There are, however, important differences.

First, the inclination of the outflows makes the island asymmetric, with electrons rotating around in a direction opposite to that of ions (when looking at the motion projection in the x-z plane). Second, as islands form and contract, the guide magnetic field is compressed and increases in strength. In the simulation with B_G = 0.5B₀, it passes from 0.5B₀ initially to 1.8−2.1B₀ in the islands, while in the simulation with B_G = B₀, it passes from B₀ initially to 2.1−2.4B₀ in the islands. Third, because of the strong magnetic field along y, temperatures are isotropized along x and z. (In the no guide field case we had Θ_zz ~ 2Θ_xx.) Here Θ_zz ≃ Θ_yy ≃ Θ_xx, equal to up to 4 for electrons and 0.2 for ions. This is a value close to the average (Θ_zz + Θ_yy + Θ_xx) / 3 of the zero guide field case. Off-diagonal terms are an order of magnitude smaller.

4.3. Reconnection electric field and reconnection rate

In the normalization of Sect. 3.6, ${\mathit{E}}^{\mathrm{\ast }}\mathrm{=}{\mathit{E}}_{\mathit{y}}\mathit{/}\mathrm{\left(}{\mathit{B}}_{\mathrm{0}}{\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{in}}^{\mathrm{R}}\mathrm{\right)}$, the Alfvén velocity includes the total magnetic field ${B}\mathrm{=}{\mathit{B}}_{\mathrm{0}}\mathrm{(}\hat{{z}}\mathrm{+}\mathit{\alpha }\hat{{y}}$). In the presence of a guide field (α ≠ 0), Alfvén waves propagating along the magnetic field do so in a direction that makes an angle θ = arctanα with the outflow direction $\hat{{z}}$. As we explained in Sect. 3.6, the reconnection electric field is in part set by the velocity at which the field lines are extracted from the X point, i.e., by the outflow velocity. This outflow velocity is mostly set by the projection of the Alfvén speed onto the outflow direction (Eq. (10) in the non-relativistic limit, Eq. (9) in the general case). Consequently, it seems justified to normalize the electric field with the projection of the Alfvén speed onto the $\hat{{z}}$ direction, i.e., to use ${\mathit{E}}^{\mathrm{\ast }}\mathrm{=}{\mathit{E}}_{\mathit{y}}\mathit{/}\mathrm{(}{\mathit{B}}_{\mathrm{0}}{\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{in}}^{\mathrm{R}}\cos \mathit{\theta }$).

Figure 14 shows the time evolution of E^∗ for three simulations with B_G/B₀ = 0, 0.5, 1, with otherwise identical parameters. The peak reconnection rate decreases when the guide field increases: 0.24, 0.22, and 0.20. This is also the kind of dependence found in PIC simulations of ion-electron, non-relativistic reconnection (Ricci et al. 2003; Huba 2005; TenBarge et al. 2014; Liu et al. 2014), relativistic pair reconnection (Hesse & Zenitani 2007), or two-fluid simulations of relativistic reconnection (Zenitani et al. 2009a)⁴. Consequently, relativistic reconnection in ion-electron plasmas does not bring new effects in this respect. However, we underline that the normalization used here allows minimizing the scatter in E^∗ for the various simulations.

5. Summary and discussion

5.1. Summary

We investigate magnetic reconnection in low β ion-electron plasmas (mass ratio of 25) with 2D PIC simulations, under relativistic conditions, i.e., the magnetic energy of the inflowing plasma exceeds its restmass energy. The simulations start from a Harris kinetic equilibrium with no localized perturbation. For diagnostics and analytical models, we use momentum and energy fluid equations based on lab-frame quantities (Appendix B). They have the advantage of being valid regardless of the particle distribution function, while the usual relativistic fluid equations using comoving quantities are restricted to isotropic comoving distribution functions.

For antiparallel reconnection, the structure of the diffusion region has a width (in the inflow direction) δ_s given by the respective inertial length d_s of the species s, measured at the center of the diffusion region. A high inflow temperature increases this width. At the entrance of the diffusion regions for simulations at low background β ≤ 2.5 × 10^-3, we find sharp transitions in the fluid quantities that were not reported elsewhere. We argue that they are not shocks. Instead, they occur when the inflowing particles have a thermal velocity that is far lower than their bulk E × B velocity, so that they enter the diffusion region with the same velocity and bounce back at the same location. We stress that these sharp transitions are not a feature of relativistic reconnection, because they depend only on the cold nature of the inflow. However, the phenomenon should be more common in relativistic reconnection because then the inflow bulk velocity v_{E × B} ~ E/B is large.

We detail the balance of terms in the relativistic Ohm’s law for antiparallel reconnection. The ion diffusion region is dominated by bulk inertia (as defined in Eq. (8)). In the electron diffusion region, bulk inertia contributes equally or more than thermal inertia. The latter result challenges the thermal-inertia-dominated paradigm that holds for non-relativistic or mildly relativistic antiparallel reconnection. We show analytically that a significant contribution of bulk inertia is to be expected whenever the inflow magnetization ${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{cold}}$ (cold, meaning that the temperature is not taken into account, see Eq. (4)) of the electrons is high, because then bulk inertia ${\mathit{\partial }}_{\mathit{z}}\mathit{p̅}{}_{\mathit{z}}\mathit{v̅}{}_{\mathit{y}}\mathrm{~}\mathit{p̅}{}_{\mathit{z}}\mathit{c}\mathit{/}\mathit{D}\mathrm{\propto }{\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{cold}}\mathit{/}\mathit{D}$ can exceed thermal inertia ${\mathit{\partial }}_{\mathit{x}}\mathit{\delta }{\mathit{p}}_{\mathit{x}}\mathit{\delta }{\mathit{v}}_{\mathit{y}}\mathrm{\propto }\mathrm{(}{\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{cold}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{/}\mathit{\delta }$. This is a new result that should hold for any large electron magnetization.

For the reconnection outflows we analytically show from mass and energy conservation that reconnection is expected to produce relativistic outflow temperatures and/or relativistic outflow bulk velocities. From simulations we find that outflow thermal energy dominates over bulk kinetic energy (from 70% to 99%, for simulations with increasing background magnetization). A more refined analytical analysis shows that this is expected if the reconnection electric field is sustained by thermal inertia. If bulk inertia dominates the thermal inertia, as expected at very high inflow magnetization, then our simple analytical model does not allow a conclusion about the cold or hot nature of the outflows. Also, our simulations do not probe high enough electron magnetizations to reach this regime: at ${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{cold}}\mathrm{=}\mathrm{90}$, thermal inertia contributes 50% of the reconnection electric field, and this fraction goes down to 25−40% at ${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{cold}}\mathrm{=}\mathrm{360}$, which is significant enough for the hypothesis of E_rec sustained by thermal inertia to hold.

For the islands we show that, with or without a guide field, their centers consist mainly of particles initially in the current sheet that were gathered inside the island during the tearing instability, that do not mix with the background plasma even after many island merging events. Particles of the background plasma cannot reach the inner parts because of the strong magnetic field surrounding the islands, and thus circle around the central part. As a result, the central part is less dense than its immediate surrounding. This may explain observed density dips at the center of magnetic islands during magnetotail reconnection events (Khotyaintsev et al. 2010), without invoking island merging or particle escape along the flux tube. Islands are also the hottest parts of the flow, with fully anisotropic temperatures in the antiparallel case, and distributions close to gyrotropic with a guide field.

We argue that the reconnection rates are to be normalized by the asymptotic magnetic field and relativistic Alfvén speed in the inflow, projected onto the outflow direction if there is a guide field: ${\mathit{E}}^{\mathrm{\ast }}\mathrm{=}{\mathit{E}}_{\mathit{y}}\mathit{/}\mathrm{\left(}{\mathit{B}}_{\mathrm{0}}{\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{in}}^{\mathrm{R}}\cos \mathit{\theta }\mathrm{\right)}$. This leads to rates in a narrow range: E^∗ peaks between 0.14–0.25. However, we stress that there is no universal value for E^∗ as defined here or elsewhere. First, because other studies show that it depends on the inflow plasma β (increasing with decreasing β). Here we find no dependence on the background plasma temperature, but smaller rates for larger particle number densities. Second, we find larger rates for the relativistic setups (0.18−0.25) than for the mildly relativistic case (0.15). These rates are also larger than those reported in the literature for ion-electron non-relativistic reconnection (0.07−0.15 for Birn et al. 2001; Pritchett 2001; Fujimoto 2006, 2009; Daughton et al. 2006; Klimas et al. 2010). This points toward relativistic reconnection being slightly faster than non-relativistic reconnection. This trend is reinforced by simulations in relativistic pair plasmas (E^∗ = 0.3 in our case or, e.g., 0.17−0.36 for Zenitani & Hoshino 2007; Bessho & Bhattacharjee 2012; Cerutti et al. 2012b). We note that this is against the interpretation of Hesse & Zenitani (2007) of a smaller rate for more relativistic setups. Third, we confirm that a guide field leads to a smaller normalized rate.

We explore the consequences of adding a guide magnetic field. The flow structure is strongly disturbed for two reasons: the Lorentz force associated with the guide field, and the relation E < B everywhere. The acceleration region is now defined by the condition E·B ≠ 0. Inflowing ion and electron Larmor radii are smaller than the island scales or magnetic gradient scales, and remain so even after the acceleration phase by E_rec, because this phase conserves the perpendicular-to-B momentum. Particles thus remain tied to the field lines everywhere, including in the acceleration region where they spend more time before being deviated in the outflows.

5.2. Discussion and astrophysical outlook

This study may serve as microphysics input for analyses on larger scales concerning magnetic energy conversion in relativistic ion-electron plasmas, as should be encountered in the coronae of AGN and microquasar accretion flows, in the lobe of radio galaxies or inside relativistic jets from AGNs or GRBs. We now discuss such applications and give estimates for key parameters in these objects: magnetic field B, electron number density n_e, magnetizations ${\mathit{\sigma }}_{\mathit{s}}^{\mathrm{cold}}$ (where cold means that only the rest mass energy is taken into account, Eq. (4)), with s = i, e for ions or electrons, and Alfvén speeds ${\mathit{V}}_{\mathrm{A}}^{\mathrm{R}}$. The properties of magnetic reconnection as studied here only depend on the inflow magnetization and temperatures, regardless of the real values of B and n_e. This is true at least as long as effects such as pair creation and annihilation, radiative braking, or Compton drag on the electrons, can be neglected (see Melzani et al. 2014b, for a discussion on these effects).

5.2.1. Toward a new regime: non-dissipative reconnection?

Our finding of a reconnection electric field sustained equally or more by bulk inertia than by thermal inertia for large inflow electron magnetization (${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{cold}}\mathrm{\ge }\mathrm{90}$), and the extrapolation of Sect. 3.8, indicate that bulk inertia might end up dominating at even higher inflow electron magnetizations. This was also envisioned by Hesse & Zenitani (2007), and reconnection in such a regime would bear significant differences with the standard picture. However, we nuance the assertion of Hesse & Zenitani (2007) that reconnection would then be a reversible process: as underlined in Sect. 3.8, the reconnection outflows may be thermally dominated even when bulk inertia dominates Ohm’s law. A definite answer to these questions requires very high magnetizations, which we can hardly reach with a PIC code and which may require relativistic gyrokinetic codes.

Highly magnetized environments, such as pulsar winds near the termination shock (${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{cold}}\mathrm{=}{\mathrm{10}}^{\mathrm{13}}$, Bucciantini et al. 2011; Sironi & Spitkovsky 2011), or other objects with ${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{cold}}⋙\mathrm{1}$ discussed in what follows, are likely to support this reconnection regime. We note that magnetar magnetospheres, even if they feature magnetizations exceeding ${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{cold}}\mathrm{=}{\mathrm{10}}^{\mathrm{13}}$ (Lyutikov & Lazarian 2013), support a very different reconnection regime due to pair and photon creations (Uzdensky 2011).

5.2.2. Large scale transient outflow production, the example of microquasars

We have shown that the reconnection outflows are thermally dominated, with a bulk Lorentz factor not necessarily increasing with the inflow magnetization and featuring low values (1.63 at most, Table 4). However, applications to large scale outflows triggered by reconnection events require some care. The outflows studied in the present manuscript originate in the electron diffusion region, and feature ion/electron decoupling. At larger distances, if not bounded by the islands and by our periodic setup, electrons and ions are expected to couple and to follow the ideal MHD dynamic. The scale on which they can propagate is fixed by larger scales than simulated here.

On the other hand, it is expected and observed (Khotyaintsev et al. 2006) that magnetic energy conversion also takes place along the magnetic separatrices far away from the dissipation region, on length scales of hundreds of ion inertial lengths. This conversion occurs through instabilities that produce thermal and non-thermal electrons (Drake et al. 2005; Egedal et al. 2009, 2012), and through the complex structure of collisionless nonlinear waves (slow shock, compound wave, rotational wave) by which the magnetized inflowing plasma transits to the hot and unmagnetized outflow on MHD scales (Liu et al. 2012; Higashimori & Hoshino 2012). It is this large-scale outflow that should be identified as the transient reconnection-driven jets in microquasar models (de Gouveia dal Pino & Lazarian 2005; de Gouveia Dal Pino et al. 2010; Kowal et al. 2011; McKinney et al. 2012; Dexter et al. 2014). In the magnetosphere close to the black hole, de Gouveia dal Pino & Lazarian (2005) estimate on the basis of an analytical model, n_e ~ 5 × 10¹⁵ cm^-3, B ~ 7 × 10⁷ G, which gives electron and ion magnetizations ${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{cold}}\mathrm{~}{\mathrm{10}}^{\mathrm{5}}$ and ${\mathit{\sigma }}_{\mathrm{i}}^{\mathrm{cold}}\mathrm{~}\mathrm{60}$, and an Alfvén speed ${\mathit{V}}_{\mathrm{A}}^{\mathrm{R}}\mathrm{~}\mathit{c}$. This is clearly in the relativistic regime. The energy content of the large scale outflows in this case has not been studied, but we can expect from the collisionless slow shocks or rotational discontinuities at the separatrices, to produce a thermally dominated outflow. It may not be so for other jet production mechanisms, and could help distinguishing in favor or against reconnection scenarios.

Another unknown is what becomes of the ambient plasma that is expelled by the first reconnected field lines ahead of the dipolarization front. In our study, it would correspond to half of a magnetic island, ejected out of the simulation box. The ambient plasma would be the plasma from the current sheet trapped in the island. As we demonstrate, this plasma does not mix with the reconnected plasma and is simply compressed and heated (see Vapirev et al. 2013, for a 3D case where instabilities imply magnetic to kinetic energy conversion). In an open configuration, it would be at the head of the large scale outflow.

5.2.3. Plasma heating in AGN and microquasar coronae and in galaxy radio lobes

Photon emission in the hard state of microquasars and AGNs is believed to come from inverse-Compton scattering of seed photons by the electrons of a corona. To achieve this, these electrons must reach temperatures on the order of 10⁹ K, i.e., Θ_e = T_e/m_ec² ~ 0.2. A non-thermal population of electrons is also required by the observation of MeV photons (Poutanen & Veledina 2014). A proposed mechanism for plasma heating is by magnetic reconnection (Di Matteo 1998; Merloni & Fabian 2001; Reis & Miller 2013). The plasma Alfvén speed estimated by these authors lies in the range 0.03c–0.3c. Associated electron magnetizations are ${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{cold}}\mathrm{~}\mathrm{1.7}\mathrm{-}\mathrm{180}$, which is in the range of the present study. A crucial question is the energy distribution between ions and electrons: if most of the magnetic energy goes to ions, and because of the low collisionality of these dilute environments, a large temperature difference can be sustained (Di Matteo et al. 1997). Our study shows that ions are slightly more heated than electrons. This can be seen with the temperatures of Figs. 8–10, and 12. More generally, the kinetic energy of the particles trapped in the magnetic islands is distributed as 55% for ions and 45% for electrons, and the kinetic energy of the particles from the background plasma that are accelerated when reaching the current sheet is also distributed as 60% for ions and 40% for electrons (for details see Melzani et al. 2014b). The energy distribution by acceleration processes far downstream of the diffusion region requires another study.

Similar questions arise concerning the heating of the lobes of radio galaxies (Kronberg et al. 2004). There, n ~ 3 × 10^-6 cm^-3 for the number densities, and B ~ 5 μG for the equipartition magnetic field with values that can be locally ten times higher, which gives magnetizations ${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{cold}}\mathrm{~}\mathrm{0.8}\mathrm{-}\mathrm{80}$ and Alfvén speeds ~ 0.02c − 0.2c. Our conclusion for the energy repartition between ions and electrons also holds.

5.2.4. Flares and “mini”-jets in extragalactic jets and in GRBs

Flare-like activity in the GeV–TeV range is observed from extragalactic jets and may possibly be explained by local reconnection events inside the jet, which produce smaller jets (the reconnection exhausts) which in turn radiate the expected photons (Giannios et al. 2009). This γ-ray emission region may be located close to the black hole ( < 0.05 pc, Giroletti et al. 2004), where B ~ 0.02−0.2 G and the plasma magnetization is high. For example, Giannios et al. (2009) take ${\mathit{\sigma }}_{\mathrm{i}}^{\mathrm{cold}}\mathrm{=}\mathrm{100}$, which leads to ${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{cold}}\mathrm{=}\mathrm{2}\hspace{0.17em}\mathrm{\times }{\mathrm{10}}^{\mathrm{5}}$ and V_A ~ c. This is in the regime where bulk inertia should dominate in Ohm’s law. Also, Giannios et al. (2009) use energy considerations to estimate that the blobs emitted from the reconnection exhausts should be ~ 10¹⁴ cm, i.e., based on its estimated particle density, 10¹⁰ ion inertial lengths. Here again, the physics far from the dissipation region should play an important role in producing such large scale structures.

5.2.5. Radio emission from extragalactic jets

Another case for relativistic magnetic reconnection is inside jets from AGNs on scales of 10−100 kpc. Radio spectra may be explained by radiation linked to reconnection events (Romanova & Lovelace 1992). Observations of AGN jets indicate B ~ 10−30 μG, n ~ 0.8−5 × 10^-8 cm^-3, and electron magnetizations in the range ${\mathit{\sigma }}_{\mathrm{e}}^{\mathrm{cold}}\mathrm{~}\mathrm{500}\mathrm{-}\mathrm{2500}$ (Schwartz et al. 2006), which implies ion magnetizations ${\mathit{\sigma }}_{\mathrm{i}}^{\mathrm{cold}}\mathrm{~}\mathrm{0.3}\mathrm{-}\mathrm{1.3}$ and Alfvén speeds ~ 0.5−0.8c. Again, our results apply in these cases, and in particular the electron magnetizations are in the very relativistic range where bulk inertia should dominate in Ohm’s law.

5.2.6. High-energy particle production

The proposed normalization of the reconnection rate, ${\mathit{E}}^{\mathrm{\ast }}\mathrm{=}{\mathit{E}}_{\mathrm{rec}}\mathit{/}\mathrm{(}{\mathit{B}}_{\mathrm{0}}{\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{in}}^{\mathrm{R}}\cos \mathit{\theta }$) with ${\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{in}}^{\mathrm{R}}$ the relativistic inflow Alfvén speed, leads to E^∗ in a close range (0.14−0.25), and because it only relies on inflow quantities, it allows for an easy prediction of the reconnection electric field. In particular, the ratio E_rec/B₀ is a key quantity for estimating the time scale of energy dissipation or the maximum energy gain for particles. It is interesting to notice that for very relativistic plasmas, ${\mathit{V}}_{\mathrm{A}\mathit{,}\mathrm{in}}^{\mathrm{R}}$ saturates at c, so that E_rec/B₀ saturates at ~ 0.2c. It may imply that the hardness of the high-energy tails saturates. We explore these matters in a forthcoming paper (Melzani et al. 2014b). Briefly, we find for a given species a power-law tail whenever its background magnetization is relativistic (above a few), with an index depending mainly on the inflow magnetization.

5.2.7. Other complications

We finally point out that the present study is oversimplified in many respects. Magnetic reconnection in magnetized coronae and jets probably often implies asymmetric plasmas from each side of the current sheet, guide fields (Aunai et al. 2013; Hesse et al. 2013; Eastwood et al. 2013), and also normal fields (along $\hat{{x}}$ here) reminiscent of the ambient magnetic field. The last point has been studied in the context of the Earth magnetotail (Pritchett 2005a, 2010; Sitnov & Swisdak 2011). Also, the initial conditions chosen in the simulations are arbitrary and do not necessarily reflect the real environments. Explored alternatives to the Harris sheet include X-point collapse (e.g., Graf von der Pahlen & Tsiklauri 2014) or force-free equilibrium (e.g., Liu et al. 2014) and show little differences with the Harris case. However, 3D initial configurations should also be considered, because in a real environment most of the energy dissipation may occur at 3D nulls, involving for example spine-fan reconnection, or at quasi-separatrix layers (Birn & Priest 2007; Pontin 2011). Few kinetic simulations exist for such setups (Baumann & Nordlund 2012; Olshevsky et al. 2013). A related point is the external forcing, i.e., the large scale plasma flow that can increase the magnetic field gradients and trigger reconnection. Studies (Pei et al. 2001; Pritchett 2005b; Ohtani & Horiuchi 2009; Klimas et al. 2010) show that the reconnection rate E^∗ is then fixed by the boundary conditions, and is thus larger than the spontaneous rate. The timescale of the forcing also proves to be important (Pei et al. 2001). These considerations, as well as some of the points made earlier on, highlight the multiscale nature of reconnection in the context of concrete astrophysical objects – and demonstrate the need for corresponding multiscale simulation studies, a field still in its infancy (e.g., Horiuchi et al. 2010; Innocenti et al. 2013)

Another central question is the validity of the 2D findings in three dimensions. Magnetic islands then become extended filaments, modulated or broken by instabilities in the third dimension or by a lack of coherence of the tearing instability (Daughton et al. 2011; Kagan et al. 2013; Markidis et al. 2013). It may imply more mixing of the current sheet particles with those of the background plasma. Concerning the validity of our claims on Ohm’s law, 3D PIC simulations in non-relativistic plasmas have shown that anomalous resistivity due to microinstabilities remains a negligible dissipation mechanism in the diffusion region (Liu et al. 2013; Karimabadi et al. 2013), where the reconnection electric field is still sustained by thermal electron inertia. The scaling analysis of Sect. 3.8 should thus remain valid, as should the conclusion that bulk inertia dominates at high inflow magnetization.

---

Acknowledgments

We thank the referee, D. Uzdensky, for comments that allowed to significantly improve the manuscript. M. Melzani would like to thank H. Baty, B. Cerutti, and M. E. Innocenti for stimulating discussions. The simulations were performed using HPC resources from GENCI (Grand Équipement National de Calcul Intensif) at CINES, CCRT, and IDRIS under the allocation x2013046960. Tests were conducted at pôle scientifique de modélisation numérique, PSMN, at the École Normale Supérieure de Lyon, whose staff we thank for their steady technical support. This work has been financially supported by the Programme National Hautes Energies (PNHE).

References

Appendix A: Relativistic Harris equilibrium

We derived the equilibrium relations for relativistic temperatures and current drift speeds, as well as for arbitrary ion to electron mass ratio and temperature ratio, in Melzani et al. (2013).

Each species follows a Maxwell-Jüttner distribution: ${\mathit{f}}_{\mathit{s}}\mathrm{\left(}{x}\mathit{,}\begin{array}{c}\mathrm{􏽥}\\ {p}\end{array}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\mu }}_{\mathit{s}}\hspace{0.17em}{\mathit{n}}_{\mathrm{cs}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}{\mathrm{4}\mathit{\pi }{\mathrm{\Gamma }}_{\mathit{s}}{\mathit{K}}_{\mathrm{2}}\mathrm{\left(}{\mathit{\mu }}_{\mathit{s}}\mathrm{\right)}}\exp \left\{\mathrm{-}{\mathit{\mu }}_{\mathit{s}}{\mathrm{\Gamma }}_{\mathit{s}}\left(\sqrt{\mathrm{1}\mathrm{+}\stackrel{}{{{\mathit{p}}^{\mathrm{2}}}_{\mathrm{˜}}}}\mathrm{-}{\mathit{U}}_{\mathit{s}}\begin{array}{c}\mathrm{˜}\\ {\mathit{p}}_{\mathit{y}}\end{array}\mathit{/}\mathit{c}\right)\right\}\mathit{,}$(A.1)with s = i for ions or e for electrons, μ_s = 1 / Θ_s = m_sc²/T_s, $\begin{array}{c}{}^{\mathrm{􏽥}}\\ {p}\end{array}\mathrm{=}{p}\mathit{/}\mathit{c}\mathrm{=}\mathit{\gamma }{v}\mathit{/}\mathit{c}$, and K₂ the modified Bessel function of the second kind. We note that f_s is indeed normalized with respect to $\begin{array}{c}{\mathrm{􏽥}}^{}\\ {p}\end{array}$ to n_cs(x) (Melzani et al. 2013), so that the particle density in the simulation frame is n_cs(x), and that in the comoving plasma frame of each species is n_0,cs(x) = n_cs(x) / Γ_s. Loading the distribution A.1 when both Θ_s and U_s are relativistic is non-trivial, and we use the method explained in Melzani et al. (2013).

For the special case where ions and electrons have the same temperatures, fulfilling Vlasov and Maxwell’s equations leads to $\begin{array}{ccc}{\mathit{n}}_{\mathrm{cs}}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \frac{{\mathit{n}}_{\mathrm{cs}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}}{{\cosh }^{\mathrm{2}}\mathrm{(}\mathit{x}\mathit{/}\mathit{L}\mathrm{)}}\mathit{,}\\ {\mathrm{\Theta }}_{\mathrm{i}}& \mathrm{=}& \mathrm{(}{\mathit{m}}_{\mathrm{e}}\mathit{/}{\mathit{m}}_{\mathrm{i}}\mathrm{)}{\mathrm{\Theta }}_{\mathrm{e}}\mathit{,}\\ {\mathrm{\Theta }}_{\mathrm{e}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{4}}{\left(\frac{{\mathit{\omega }}_{\mathrm{ce}}}{{\mathit{\omega }}_{\mathrm{0}\mathit{,}\mathrm{pe}}}\right)}^{\mathrm{2}}\mathit{,}\\ \frac{{\mathrm{\Gamma }}_{\mathit{s}}{\mathit{U}}_{\mathit{s}}}{\mathit{c}}& \mathrm{=}& \mathrm{-}\mathrm{2}{\mathrm{\Theta }}_{\mathit{s}}\frac{{\mathit{d}}_{\mathrm{0}\mathit{,}\mathrm{e}}}{\mathit{L}}\frac{{\mathit{\omega }}_{\mathrm{0}\mathit{,}\mathrm{pe}}}{{\mathit{\omega }}_{\mathrm{c}\mathit{s}}}\mathrm{sgn}\mathrm{\left(}{\mathit{q}}_{\mathit{s}}\mathrm{\right)}\mathit{,}\end{array}$with sgn(q_s) the sign of the charge q_s, ω_cs = eB₀/m_s the cyclotron pulsation defined in the asymptotic magnetic field B₀ (e> 0 here), ${\mathit{\omega }}_{\mathrm{0}\mathit{,}\mathrm{pe}}\mathrm{=}\sqrt{{\mathit{n}}_{\mathrm{0}\mathit{,}\mathrm{cs}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}{\mathit{e}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}{\mathit{ϵ}}_{\mathrm{0}}{\mathit{m}}_{\mathrm{e}}\mathrm{\right)}}$ the electron plasma pulsation defined by the comoving number density n_0,cs(0) = n_cs(0) / Γ_s at the center of the current sheet, and d_0,e = c/ω_0,pe the associated inertial length.

Inserting Eqs. (A.2c) into (A.2d), the latter becomes $\frac{{\mathrm{\Gamma }}_{\mathit{s}}{\mathit{U}}_{\mathit{s}}}{\mathit{c}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{d}}_{\mathrm{0}\mathit{,}\mathrm{e}}}{\mathit{L}}\frac{{\mathit{\omega }}_{\mathrm{ce}}}{{\mathit{\omega }}_{\mathrm{0}\mathit{,}\mathrm{pe}}}\mathrm{sgn}\mathrm{\left(}{\mathit{q}}_{\mathit{s}}\mathrm{\right)}\mathit{,}$(A.3)so that we see that Eqs. (A.2a)−(A.2c) and (A.3) are well suited to expressing the equilibrium relations in terms of the comoving quantities ω_ce/ω_0,pe and L/d_0,e alone. Some manipulations are needed to express everything in terms of the lab-frame quantities ω_ce/ω_pe and L/d_e, where ω_pe is the counterpart of ω_0,pe in the simulation frame (${\mathit{\omega }}_{\mathrm{pe}}\mathrm{=}\sqrt{{\mathrm{\Gamma }}_{\mathrm{e}}}{\mathit{\omega }}_{\mathrm{0}\mathit{,}\mathrm{pe}}$), and d_e = c/ω_pe. To do so, we note in Eq. (A.2d) that ω_0,ped_0,e = c = ω_ped_e, and we express Eq. (A.2c) as Θ_e = Γ_e(ω_ce/ω_pe)²/ 4. With this, we obtain $\begin{array}{ccc}\frac{{\mathit{U}}_{\mathit{s}}}{\mathit{c}}& \mathrm{=}& \\ {\mathrm{\Theta }}_{\mathrm{e}}& \mathrm{=}& \end{array}$We see from Eq. (A.4a) that we have the condition ω_ce/ω_pe < 2L/d_e. If this is not the case, the equilibrium cannot be achieved. Why it is so can be seen by rewriting Eq. (A.2d) as Θ_e = (L/ 2d_e)(ω_ce/ω_pe)Γ_eU_e/c: satisfying Vlasov equation is possible only if Θ_e < (L/ 2d_e)(ω_ce/ω_pe)Γ_e, but this is not possible if the pressure balance condition (which is Eq. (A.4b)) requires a temperature exceeding this limit to balance the magnetic field pressure.

As a final note, we express the thermal Larmor radius of the particles, defined as ⟨ (γv_perp)² ⟩ ^{1 / 2}/ω_ce where ⟨·⟩ is an average over the distribution function, at current sheet center: $\frac{\mathrm{⟨}{\mathit{r}}_{\mathrm{ce}}\mathrm{⟩}}{{\mathit{d}}_{\mathrm{e}}}\mathrm{=}\frac{{\mathit{\omega }}_{\mathrm{pe}}}{{\mathit{\omega }}_{\mathrm{ce}}}\sqrt{{\mathrm{\Theta }}_{\mathrm{e}}}\sqrt{{\mathit{\kappa }}_{\mathrm{32}}\mathrm{\left(}{\mathit{\mu }}_{\mathrm{e}}\mathrm{\right)}}\mathrm{=}\sqrt{{\mathrm{\Gamma }}_{\mathrm{e}}{\mathit{\kappa }}_{\mathrm{32}}\mathrm{\left(}{\mathit{\mu }}_{\mathrm{e}}\mathrm{\right)}\mathit{/}\mathrm{2}}\mathit{,}$(A.5)where the first part of the equality is general (Melzani et al. 2013), and the second is obtained for the Harris equilibrium using Eq. (A.4b) for the ratio ω_ce/ω_pe. The function κ₃₂ is the plasma comoving enthalpy and is plotted in Fig. B.1. The thermal Larmor radius is consequently temperature dependent via κ₃₂(1 / Θ_e).

Appendix B: From Vlasov to fluid equations

Appendix B.1: Fluid equations

Fluid equations employed in numerical codes are usually expressed in terms of comoving quantities, such as, for species s, the comoving particle number density n_0s, the comoving enthalpy h_0s, the comoving pressure P_0s, and in terms of the fluid velocity $\overline{{v}}{}_{\mathit{s}}$ and its associated Lorentz factor Γ_s. The conservation of particle number and of momentum for each species, and of total energy, then read as (see e.g., Mihalas & Mihalas 1984; Barkov et al. 2014) $\begin{array}{cc}& \\ \hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{(}{\mathrm{\Gamma }}_{\mathit{s}}^{\mathrm{2}}{\mathit{n}}_{\mathrm{0}\mathit{s}}{\mathit{h}}_{\mathrm{0}\mathit{s}}{\overline{{v}}}_{\mathit{s}}\mathrm{)}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial }{x}}\mathrm{·}\mathrm{(}{\mathit{n}}_{\mathrm{0}\mathit{s}}{\mathit{h}}_{\mathrm{0}\mathit{s}}{\mathrm{\Gamma }}_{\mathit{s}}\overline{{v}}{}_{\mathit{s}}\mathrm{\Gamma }_{\mathit{s}}\overline{{v}}{\mathit{s}}_{}\mathrm{)}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{{\mathit{m}}_{\mathit{s}}}\frac{\mathit{\partial }{\mathit{P}}_{\mathrm{0}\mathit{s}}}{\mathit{\partial }{x}}\mathrm{+}\frac{{\mathit{q}}_{\mathit{s}}}{{\mathit{m}}_{\mathit{s}}}{\mathrm{\Gamma }}_{\mathit{s}}{\mathit{n}}_{\mathrm{0}\mathit{s}}\mathrm{(}{E}\mathrm{+}{\overline{{v}}}_{\mathit{s}}\mathrm{\wedge }{B}\mathrm{)}\mathit{,}\end{array}& \end{array}$These equations are, however, not well suited to analysis of particle simulations. First, because accessible quantities are those defined in the simulation (or lab) frame, while those in the comoving frame of the plasma must be obtained by a boost at the local mean speed $\overline{{v}}{}_{\mathit{s}}$. Second, because they assume a comoving particle distribution that is isotropic in momentum space in order to use a scalar pressure P_0s instead of the full pressure tensor and, as we show below, to write relations such as $\overline{{p}}{}_{\mathit{s}}\mathrm{=}\mathrm{⟨}\mathit{\gamma }{v}{\mathrm{⟩}}_{\mathit{s}}\mathrm{=}{\mathit{h}}_{\mathrm{0}\mathit{s}}{\mathrm{\Gamma }}_{\mathit{s}}\overline{{v}}{}_{\mathit{s}}$ for the mean momentum. This is not the case in the out-of-equilibrium plasmas that we study.

Instead, we derive the fluid equations directly from the collisionless Vlasov equation. The latter reads as $\frac{\mathit{\partial }{\mathit{f}}_{\mathit{s}}\mathrm{\left(}{x}\mathit{,}{p}\mathit{,t}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{+}{v}\mathrm{·}\frac{\mathit{\partial }{\mathit{f}}_{\mathit{s}}}{\mathit{\partial }{x}}\mathrm{+}\frac{{\mathit{q}}_{\mathit{s}}}{{\mathit{m}}_{\mathit{s}}}\mathrm{(}{E}\mathrm{+}{v}\mathrm{\wedge }{{B}}^{\mathrm{)}}\mathrm{·}\frac{\mathit{\partial }{\mathit{f}}_{\mathit{s}}}{\mathit{\partial }{p}}\mathrm{=}\mathrm{0}\mathit{,}$(B.2)where p = γv, and f_s is the distribution function in the simulation or lab frame. We denote its counterpart in the comoving frame by f_0s.

The first moment (with 1) of Eq. (B.2) gives the equation of conservation of the number of particles: $\frac{\mathit{\partial }}{\mathit{\partial t}}{\mathit{n}}_{\mathrm{lab}\mathit{,s}}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial }{x}}\mathrm{·}\mathrm{\left(}{\mathit{n}}_{\mathrm{lab}\mathit{,s}}\overline{{v}}\mathit{s}\mathrm{\right)}\mathrm{=}\mathrm{0.}$(B.3)We note that n_lab,s = Γ_sn_0s, so that we indeed recover Eq. (B.1).

The second moment (with p) gives the equation of conservation of momentum: $\begin{array}{ccc}\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{\left(}{\mathit{n}}_{\mathrm{lab}\mathit{,s}}\overline{{p}}\mathit{s}\mathrm{\right)}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial }{x}}\mathrm{·}\mathrm{\left(}{\mathit{n}}_{\mathrm{lab}\mathit{,s}}\overline{{p}}\mathit{s}\overline{{v}}\mathit{s}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{\partial }}{\mathit{\partial }{x}}\mathrm{·}\mathrm{(}{\mathit{n}}_{\mathrm{lab}\mathit{,s}}\mathrm{⟨}\mathit{\delta }{{p}}_{\mathit{s}}\mathit{\delta }{{v}}_{\mathit{s}}{\mathrm{⟩}}_{\mathit{s}}\mathrm{)}\mathrm{+}\frac{{\mathit{q}}_{\mathit{s}}}{{\mathit{m}}_{\mathit{s}}}{\mathit{n}}_{\mathrm{lab}\mathit{,s}}\mathrm{(}{E}\mathrm{+}\overline{{v}}\mathit{s}\mathrm{\wedge }{B}\mathrm{)}\mathit{.}& & \end{array}$(B.4)We used the definition $\mathit{\delta }{p}\mathrm{=}{p}\mathrm{-}\overline{{p}}{}_{\mathit{s}}$, where p = γv is the momentum, and similarly for δv. Also, ⟨ ·⟩ _s denotes an average in p over the distribution function f_s. To recover Eq. (B.1a), we use the relation $\overline{{p}}{}_{\mathit{s}}\mathrm{=}{\mathit{h}}_{\mathrm{0}\mathit{s}}{\mathrm{\Gamma }}_{\mathit{s}}\overline{{v}}{}_{\mathit{s}}$ for the mean momentum, which is only valid if the comoving distribution f_0s is isotropic in p (Melzani et al. 2013). Also, we note that the stress tensor is defined as , and is also equal (again if f_0s is isotropic) to ${\mathrm{\Pi }}_{\mathit{ij}}\mathrm{=}{\mathit{P}}_{\mathrm{0}\mathit{s}}{\mathit{\delta }}_{\mathit{ij}}\mathrm{+}{\mathrm{\Gamma }}_{\overline{{v}}{}_{\mathit{s}}}^{\mathrm{2}}{\mathit{h}}_{\mathrm{0}\mathit{s}}{\mathit{n}}_{\mathrm{0}\mathit{s}}{\mathit{m}}_{\mathit{s}}\mathit{v̅}{}_{\mathit{s,i}}\mathit{v̅}{}_{\mathit{s,j}}$. Inserting these expressions into Eq. (B.4) does lead to Eq. (B.1a). We also note that Eq. (B.4) can be put into a conservative form by using the conservation of momentum for the electromagnetic field, which reads as $\frac{\mathit{\partial }}{\mathit{\partial t}}{\mathit{ϵ}}_{\mathrm{0}}\mathrm{(}{E}\mathrm{\wedge }{B}\mathrm{)}\mathrm{-}\frac{\mathit{\partial }}{\mathit{\partial }{x}}\mathrm{·}{T}\mathrm{=}\mathrm{-}\mathrm{[}\mathit{\rho }{E}\mathrm{+}{j}\mathrm{\wedge }{B}\mathrm{]}\mathit{,}$(B.5)where ${\mathit{T}}_{\mathit{ij}}\mathrm{=}{\mathit{ϵ}}_{\mathrm{0}}\mathrm{(}{\mathit{E}}_{\mathit{i}}{\mathit{E}}_{\mathit{j}}\mathrm{-}{\mathit{E}}^{\mathrm{2}}{\mathit{\delta }}_{\mathit{ij}}\mathit{/}\mathrm{2}\mathrm{)}\mathrm{+}{\mathit{\mu }}_{\mathrm{0}}^{-1}\mathrm{(}{\mathit{B}}_{\mathit{i}}{\mathit{B}}_{\mathit{j}}\mathrm{-}{\mathit{B}}^{\mathrm{2}}{\mathit{\delta }}_{\mathit{ij}}\mathit{/}\mathrm{2}\mathrm{)}$ is Maxwell stress tensor, ρ = ∑ _sq_sn_lab,s is the charge density, and ${j}\mathrm{=}{\sum }_{\mathit{s}}{\mathit{q}}_{\mathit{s}}{\mathit{n}}_{\mathrm{lab}\mathit{,s}}\overline{{v}}{}_{\mathit{s}}$ is the current density. One thus has to sum Eq. (B.4) over all species and then use Eq. (B.5), to obtain $\begin{array}{ccc}\frac{\mathit{\partial }}{\mathit{\partial t}}\left({\mathit{ϵ}}_{\mathrm{0}}{E}\mathrm{\wedge }{B}\mathrm{+}\sum _{\mathit{s}}{\mathit{m}}_{\mathit{s}}{\mathit{n}}_{\mathrm{lab}\mathit{,s}}\overline{{p}}\mathit{s}\right)\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial }{x}}\mathrm{·}\left(\mathrm{-}{T}\mathrm{+}\sum _{\mathit{s}}{\mathit{m}}_{\mathit{s}}{\mathit{n}}_{\mathrm{lab}\mathit{,s}}\left[\overline{{p}}\mathit{s}\overline{{v}}\mathit{s}\mathrm{+}\mathrm{⟨}\mathit{\delta }{{p}}_{\mathit{s}}\mathit{\delta }{{v}}_{\mathit{s}}{\mathrm{⟩}}_{\mathit{s}}\right]\right)\mathrm{=}\mathrm{0.}& & \end{array}$(B.6)Finally, multiplying Vlasov Eq. (B.2) by γm_sc² and integrating over p gives the equation of conservation of energy: $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{(}{\mathit{n}}_{\mathrm{lab}\mathit{,s}}\mathrm{⟨}\mathit{\gamma }{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{s}}\mathrm{)}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial }{x}}\mathrm{·}\mathrm{(}{\mathit{n}}_{\mathrm{lab}\mathit{,s}}\mathrm{⟨}{v}\mathit{\gamma }{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{s}}\mathrm{)}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{=}{\mathit{q}}_{\mathit{s}}{\mathit{n}}_{\mathrm{lab}\mathit{,s}}\mathrm{⟨}{E}\mathrm{·}{v}{\mathrm{⟩}}_{\mathit{s}}\mathit{.}\end{array}$(B.7)The righthand side accounts for the coupling between the species and the electromagnetic fields, and thus possibly with other species via collective interactions. The non-relativistic limit of this equation is easily obtained by making the difference between Eqs. (B.7) and (B.3). Also, Eq. (B.7) can be put into a useful conservative form by expressing its righthand side through the energy equation for the fields, which is $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\frac{\mathit{\partial }}{\mathit{\partial t}}\left(\frac{{\mathit{E}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}{\mathit{c}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{B}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}}\right)\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial }{x}}\mathrm{·}\frac{{E}\mathrm{\wedge }{B}}{{\mathit{\mu }}_{\mathrm{0}}}\mathrm{=}\mathrm{-}{E}\mathrm{·}{j}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{=}\mathrm{-}{E}\mathrm{·}\sum _{\mathit{s}}\mathrm{\int }{\mathrm{d}}^{\mathrm{3}}{p}{\mathit{f}}_{\mathit{s}}\mathrm{\left(}{x}\mathit{,}{p}\mathrm{\right)}{\mathit{q}}_{\mathit{s}}{v}\\ \begin{array}{c}\end{array}\mathrm{=}\mathrm{-}\sum _{\mathit{s}}{\mathit{n}}_{\mathrm{lab}\mathit{,}\mathrm{s}}{\mathit{q}}_{\mathit{s}}\mathrm{⟨}{v}\mathrm{·}{E}{\mathrm{⟩}}_{\mathit{s}}\mathit{.}\end{array}$(B.8)We thus have to sum Eq. (B.7) for all species, to obtain $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\frac{\mathit{\partial }}{\mathit{\partial t}}\left\{\sum _{\mathit{s}}\mathrm{(}{\mathit{n}}_{\mathrm{lab}\mathit{,}\mathrm{s}}\mathrm{⟨}\mathit{\gamma }{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{s}}\mathrm{)}\mathrm{+}\frac{{\mathit{E}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}{\mathit{c}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{B}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\mu }}_{\mathrm{0}}}\right\}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial }{x}}\mathrm{·}\left\{\sum _{\mathit{s}}\mathrm{(}{\mathit{n}}_{\mathrm{lab}\mathit{,}\mathrm{s}}\mathrm{⟨}{v}\mathit{\gamma }{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{s}}\mathrm{)}\mathrm{+}\frac{{E}\mathrm{\wedge }{B}}{{\mathit{\mu }}_{\mathrm{0}}}\right\}\mathrm{=}\mathrm{0.}\end{array}$(B.9)To recover Eq. (B.1b), we use the relation ${\mathit{n}}_{\mathrm{lab}\mathit{,s}}\mathrm{⟨}\mathit{\gamma }{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}{\mathrm{⟩}}_{\mathit{s}}\mathrm{=}{\mathrm{\Gamma }}_{\mathit{s}}^{\mathrm{2}}{\mathit{n}}_{\mathrm{0}\mathit{s}}{\mathit{h}}_{\mathrm{0}\mathit{s}}{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}\mathrm{-}{\mathit{P}}_{\mathrm{0}\mathit{s}}$ (Melzani et al. 2013, Table 5), and the relations previously used for $\overline{{p}}{}_{\mathit{s}}$ and n_lab,s.

Appendix B.2: Measure of the temperature tensor

The kinetic temperature is a second-order moment of the particle distribution function and is, as such, not always suited to characterizing the velocity distribution of a population strongly out-of-equilibrium. We nevertheless use it as an indication of the thermal energy content of the population, and of the momentum flux transfers, the latter being especially relevant to Ohm’s law. We define the temperature tensor as the ratio of the pressure tensor n_lab,sm_s ⟨ δp_sδv_s ⟩ _s (appearing in the equation of conservation of momentum, Eq. (B.4)), to the comoving particle number density n_0s. Since n_lab,s = Γ_sn_0,s, the temperature tensor is $\frac{{\mathit{T}}_{\mathit{ij,s}}}{{\mathit{m}}_{\mathit{s}}{\mathit{c}}^{\mathrm{2}}}\mathrm{=}{\mathrm{\Gamma }}_{\mathit{s}}\frac{\mathrm{⟨}\mathit{\delta }{\mathit{p}}_{\mathit{i}}\mathit{\delta }{\mathit{v}}_{\mathit{j}}{\mathrm{⟩}}_{\mathit{s}}}{{\mathit{c}}^{\mathrm{2}}}\mathrm{·}$(B.10)For the special case of a Maxwell-Jüttner distribution function of temperature Θ_s = 1 /μ_s = T_s/ (m_sc²), of drift velocity β_0,s and associated Lorentz factor Γ_0,s, given in the simulation frame by ${\mathit{f}}_{\mathit{s}}\mathrm{\left(}{x}\mathit{,}\begin{array}{c}{}_{\mathrm{˜}}\\ {p}\end{array}\mathrm{\right)}\mathrm{=}{\mathit{n}}_{\mathrm{lab}\mathit{,s}}\mathrm{\left(}{x}\mathrm{\right)}{\mathit{g}}_{\mathit{s}}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {p}\end{array}\mathrm{\right)}$, $\begin{array}{c}{}_{\mathrm{˜}}\\ {p}\end{array}\mathrm{=}\mathit{\gamma }{v}\mathit{/}\mathit{c}$, and ${\mathit{g}}_{\mathit{s}}\mathrm{\left(}\begin{array}{c}\mathrm{˜}\\ {p}\end{array}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\mu }}_{\mathit{s}}}{\mathrm{4}\mathit{\pi }{\mathit{K}}_{\mathrm{2}}\mathrm{\left(}{\mathit{\mu }}_{\mathit{s}}\mathrm{\right)}{\mathrm{\Gamma }}_{\mathrm{0}\mathit{,s}}}\exp \left\{\mathrm{-}\mathit{\mu }{\mathrm{\Gamma }}_{\mathrm{0}\mathit{,s}}\left(\sqrt{\mathrm{1}\mathrm{+}\stackrel{}{{{\mathit{p}}^{\mathrm{2}}}_{\mathrm{˜}}}}\mathrm{-}\overline{{v}}\mathit{s}\mathit{/}\mathit{c}\mathrm{·}\begin{array}{c}\mathrm{˜}\\ {p}\end{array}\right)\right\}\mathit{,}$(B.11)we do have (Melzani et al. 2013) ⟨ δp_iδv_j ⟩ _s/c² = Θ_sδ_ij/ Γ_s where δ_ij is Kronecker’s delta. We also note that this distribution is normalized to unity: ${}^{\#}\mathrm{d}^{\mathrm{3}}\begin{array}{c}{}_{\mathrm{˜}}\\ {p}\end{array}\hspace{0.17em}{\mathit{g}}_{\mathit{s}}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {p}\end{array}\mathrm{\right)}\mathrm{=}\mathrm{1}$. Based upon these considerations, we measure the temperatures with expression (B.10).

All Tables

All Figures