Relativistic Magnetic Reconnection in Electron-Positron-Proton Plasmas: Implications for Jets of Active Galactic Nuclei

TL;DR

This study uses particle-in-cell simulations to explore magnetic reconnection in electron-positron-proton plasmas, revealing how energy is distributed among particles and magnetic fields, with implications for understanding AGN jet emissions.

Contribution

It investigates the largely unexplored regime of pair-proton plasma reconnection, providing empirical relations for energy partition and particle spectra relevant to AGN jets.

Findings

Energy is shared roughly equally among magnetic fields, pairs, and protons at high magnetization.

Lepton Lorentz factors depend on magnetization, temperature, and plasma composition, following a specific empirical relation.

Reconnection produces power-law lepton energy distributions consistent with AGN jet observations.

Abstract

Magnetic reconnection is often invoked to explain the non-thermal radiation of relativistic outflows, including jets of active galactic nuclei (AGN). Motivated by the largely unknown plasma composition of AGN jets, we study reconnection in the unexplored regime of electron-positron-proton (pair-proton) plasmas with large-scale two-dimensional particle-in-cell simulations. We cover a wide range of pair multiplicities (lepton-to-proton number ratio $\kappa=1-199$) for different values of the all-species plasma magnetization ($\sigma=1,3$ and 10) and electron temperature ($\Theta_e\equiv kT_e/m_ec^2=0.1-100$). We focus on the dependence of the post-reconnection energy partition and lepton energy spectra on the hot pair plasma magnetization $\sigma_{e,h}$ (i.e., the ratio of magnetic to pair enthalpy densities). We find that the post-reconnection energy is shared roughly equally between magnetic fields, pairs, and protons for $\sigma_{e,h}\gtrsim 3$. We empirically find that the mean lepton Lorentz factor in the post-reconnection region depends on $\sigma, \Theta_e$, and $\sigma_{e,h}$ as $\langle \gamma_e-1\rangle \approx \sqrt{\sigma}(1+4\Theta_e)\left(1+\sigma_{e,h}/30\right)$, for $\sigma\ge1$. The high-energy part of the post-reconnection lepton energy distributions can be described by a power law, whose slope is mainly controlled by $\sigma_{e,h}$ for $\kappa \gtrsim 3-6$, with harder power laws obtained for higher magnetizations. We finally show that reconnection in pair-proton plasmas with multiplicities $\kappa \sim 1-20$, magnetizations $\sigma \sim 1-10$, and temperatures $\Theta_e \sim 1-10$ results in particle power law slopes and average electron Lorentz factors that are consistent with those inferred in leptonic models of AGN jet emission.

Figures

Tables

Table 1. Table 1: Simulation parameters.

Run	$\sigma$	${\Theta }_e$	$\kappa$	${\sigma }_{e,h}$aaHot pair plasma magnetization defined in Eq. (1).	$L/{\rho }_{Le}$	$L/{\rho }_{Li}$	$T_{\max }$bbDuration of the simulation in units of $L/c$.
A0*	1	1	199	2.9	2706.5	13.6	1.8
A1	1	1	66	6.9	1769.0	26.9	1.8
A2	1	1	19	21.4	1004.7	52.9	1.8
A3	1	1	6	69.3	273.2	48.2	6.1ccThe reconnection rate decreases after $\sim 3L/c$ due to the formation of a large boundary island.
A4*	1	1	6	69.3	557.9	98.4	1.7
A5	1	1	6	69.3	1195.4	210.9	1.3
A6	1	1	3	130.0	398.9	133.0	1.7
A7	1	1	1.2	317.6	257.7	206.6	1.7
A8	1	1	1	387.9	245.2	245.2	5.5
A9	1	1	199	2.9	5799.8	29.1	1.5
B0	1	10	199	1.2	4099.7	20.6	2.0
B1	1	10	66	1.7	3492.7	53.2	1.9
B2	1	10	19	3.4	2466.1	129.8	1.9
B3	1	10	6	9.0	1510.7	266.6	1.8
B4	1	10	3	16.1	1104.6	368.2	1.8
B5	1	10	1	46.4	688.2	688.2	3.1
C1	3	1	199	8.8	1562.6	7.85	1.7
C2	3	1	66	20.7	1021.3	15.5	1.6
C3	3	1	19	64.1	580.1	30.5	1.5
C4*	3	1	6	207.8	322.1	56.8	3.1
C5	3	1	6	207.8	690.2	121.8	1.3
C6	3	1	1	1163.8	141.5	141.5	3.9ccThe reconnection rate decreases after $\sim 3L/c$ due to the formation of a large boundary island.
D1	3	10	66	5.1	1975.4	30.1	1.2
D2	3	10	19	10.2	1423.8	74.9	1.7
D3	3	10	6	27.1	872.2	153.9	1.5
D4	3	10	1	139.3	397.3	397.3	3.9
E1	10	1	199	29.4	838.4	4.2	1.5
E2	10	1	19	213.6	311.2	16.4	1.4
E3*	10	1	1	3879.2	77.5	77.5	3.9ccThe reconnection rate decreases after $\sim 3L/c$ due to the formation of a large boundary island.
E4	10	1	1	3879.2	229.8	229.8	1.0
F1	10	10	199	12.9	1296.4	6.5	1.6
F2	10	10	6	90.2	468.0	82.6	1.5
F3	10	10	1	464.5	217.6	217.6	3.9
G1	1	0.1	19	76.4	1230.6	64.8	1.5
G2	1	0.1	3	478.3	491.7	163.9	1.3
H1	1	100	199	1.0	4376.9	22.0	1.9
H2	1	100	19	1.3	3911.1	205.8	2.0
H3	1	100	3	2.7	2610.6	870.2	2.0
H4	1	100	1	6.2	1758.0	1758.0	6.0

Table 2. Table 2: Summary table of the power-law index of the pair energy distribution p 𝑝 p , mean electron Lorentz factor ⟨ γ e − 1 ⟩ delimited-⟨⟩ subscript 𝛾 𝑒 1 \langle\gamma_{e}-1\rangle , hot pair plasma magnetization σ e , h subscript 𝜎 𝑒 ℎ \sigma_{e,h} (Eq. ( 1 )), and electron plasma parameter β e subscript 𝛽 𝑒 \beta_{e} (Eq. ( A7 )) from our simulations of reconnection in pair-proton plasmas (see also Figs.  12 and 13 ). A systematic error of ± 0.2 plus-or-minus 0.2 \pm 0.2 applies to all p 𝑝 p values. Results from the box-size scaling simulations are not included here.

Run	$\sigma$	${\Theta }_e$	$\kappa$	${\sigma }_{e,h}$	${\beta }_e$	$p$	$⟨{\gamma }_e-1⟩$
A0	1	1	199	2.9	0.072	4.8	6.0
A1	1	1	66	6.9	0.031	3.5	7.0
A2	1	1	19	21.4	0.010	2.5	9.8
A4	1	1	6	69.3	0.003	2.0	15.8
A6	1	1	3	130.0	0.002	2.1	25.8
A7	1	1	1.2	317.6	0.001	2.4	75.0
A8	1	1	1	387.9	0.001	2.5	147.5
B0	1	10	199	1.2	0.199	4.4	44.1
B1	1	10	66	1.7	0.146	4.6	48.2
B2	1	10	19	3.4	0.076	4.6	53.9
B3	1	10	6	9.0	0.032	3.5	65.9
B4	1	10	3	16.1	0.020	3.4	76.3
B5	1	10	1	46.4	0.010	4.0	201.3
C1	3	1	199	8.8	0.024	3.6	13.6
C2	3	1	66	20.7	0.010	2.5	20.0
C3	3	1	19	64.1	0.003	1.9	36.3
C4	3	1	6	207.8	0.001	1.7	79.8
C6	3	1	1	1163.8	0.0004	2.0	602.9
D1	3	10	66	5.1	0.049	4.0	86.0
D2	3	10	19	10.2	0.025	3.1	107.9
D3	3	10	6	27.1	0.011	2.5	153.7
D4	3	10	1	139.3	0.003	3.0	691.2
E1	10	1	199	29.4	0.007	2.4	40.6
E2	10	1	19	213.6	0.001	1.6	166.3
E3	10	1	1	3879.2	0.0001	1.4aafootnotemark:	2114.7
F1	10	10	199	12.9	0.020	3.0	203.1
F2	10	10	6	90.2	0.003	2.0	592.4
F3	10	10	1	464.5	0.001	1.8	2371.3
G1	1	0.1	19	76.4	0.001	1.9	0.31
G2	1	0.1	3	478.3	0.0002	1.6	0.25
H1	1	100	199	1.0	0.244	3.2	454.0
H2	1	100	19	1.3	0.205	3.3	463.1
H3	1	100	3	2.7	0.121	3.6	512.2

Table 3. Table 3: Description, symbol, and definition of parameters used in this study.

Parameter	Symbol	Definition
Pair multiplicity	$\kappa$	$n_{e^{\pm }}/n_i$
Proton fraction	$q$	$n_i/n_{e^{-}}=2/\left(\kappa +1\right)$
Lepton adiabatic index	${\hat{\gamma }}_e$	Synge (1957)
Proton adiabatic index	${\hat{\gamma }}_i$	Synge (1957)
Total plasma magnetization	$\sigma$	Eq. (A1)
Hot pair plasma magnetization	${\sigma }_{e,h}$	Eq. (1)
Cold pair plasma magnetization	${\sigma }_{e,c}$	Eq. (A4)
Hot proton plasma magnetization	${\sigma }_{i,h}$	Eq. (A5)
Cold proton plasma magnetization	${\sigma }_{i,c}$	Eq. (A6)
Electron plasma $\beta$	${\beta }_e$	Eq. (A7)
Plasma electron frequency	${\omega }_{pe-}$	Eq. (A9)
Plasma proton frequency	${\omega }_{pi}$	Eq. (A11)
Electron Larmor radius	${\rho }_{Le}$	Eq. (A12)
Proton Larmor radius	${\rho }_{Li}$	Eq. (A13)

Table 4. Table 4: Parameter values (with their 1 σ 𝜎 \sigma statistical errors) obtained from a χ 2 superscript 𝜒 2 \chi^{2} fit of Eq. ( D1 ) to the mean lepton Lorentz factor derived from our simulations for different σ 𝜎 \sigma and Θ e subscript Θ 𝑒 \Theta_{e} values. We exclude cases with less data points than the number of free model parameters.

$\sigma$	${\Theta }_e$	$a$	$\chi$	$b$
1	1	$114.9\pm 18.1$	$1.5\pm 0.2$	$6.5\pm 0.7$
1	10	$149.7\pm 17.6$	$1.3\pm 0.2$	$46.8\pm 2.2$
1	100	$328.9\pm 7.1$	$1.6\pm 0.1$	$456.7\pm 3.4$
3	1	$565.7\pm 73.6$	$1.1\pm 0.1$	$13.5\pm 2.1$
3	10	$603.7\pm 36.7$	$1.2\pm 0.1$	$85.7\pm 7.0$

Equations

$$\frac{{\sigma }_{e,h}}{\sigma }=\frac{q\left(\frac{m_i}{m_e}+{\Theta }_e\frac{{\hat{\gamma }}_i}{{\hat{\gamma }}_i-1}\right)+(2-q)\left(1+{\Theta }_e\frac{{\hat{\gamma }}_e}{{\hat{\gamma }}_e-1}\right)}{(2-q)\left(1+{\Theta }_e\frac{{\hat{\gamma }}_e}{{\hat{\gamma }}_e-1}\right)},$$

$$⟨{\gamma }_e-1⟩\approx \sqrt{\sigma }\left(1+4{\Theta }_e\right)\left(1+\frac{{\sigma }_{e,h}}{30}\right)\cdot$$

$$\displaystyle f(\gamma)\propto\left\{\begin{array}[]{cc}\gamma^{-s_{1}},&\gamma\leq\gamma_{inj}\\ \gamma^{-s_{2}},&\gamma>\gamma_{inj}\end{array}\right.$$

$$\sigma$$

$${\sigma }_c=\frac{B_0^2}{4\pi \left(n_i{m}_i{c}^2+n_{e^{\pm }}{m}_e{c}^2\right)}\cdot$$

$${\sigma }_{e,h}=\sigma \frac{q\left(\frac{m_i}{m_e}+\frac{{\hat{\gamma }}_i{\Theta }_e}{{\hat{\gamma }}_i-1}\right)+(2-q)\left(1+\frac{{\hat{\gamma }}_e{\Theta }_e}{{\hat{\gamma }}_e-1}\right)}{(2-q)\left(1+\frac{{\hat{\gamma }}_e{\Theta }_e}{{\hat{\gamma }}_e-1}\right)}\cdot$$

$${\sigma }_{e,c}=\sigma \frac{q\left(\frac{m_i}{m_e}+\frac{{\hat{\gamma }}_i{\Theta }_e}{{\hat{\gamma }}_i-1}\right)+(2-q)\left(1+\frac{{\hat{\gamma }}_e{\Theta }_e}{{\hat{\gamma }}_e-1}\right)}{2-q}\cdot$$

$${\sigma }_{i,h}=\sigma \frac{q\left(\frac{m_i}{m_e}+\frac{{\hat{\gamma }}_i{\Theta }_e}{{\hat{\gamma }}_i-1}\right)+(2-q)\left(1+\frac{{\hat{\gamma }}_e{\Theta }_e}{{\hat{\gamma }}_e-1}\right)}{q\left(\frac{m_i}{m_e}+\frac{{\hat{\gamma }}_i{\Theta }_e}{{\hat{\gamma }}_i-1}\right)},$$

$${\sigma }_{i,c}=\frac{\sigma m_e}{q{m}_i}\left[q\left(\frac{m_i}{m_e}+\frac{{\hat{\gamma }}_i{\Theta }_e}{{\hat{\gamma }}_i-1}\right)+(2-q)\left(1+\frac{{\hat{\gamma }}_e{\Theta }_e}{{\hat{\gamma }}_e-1}\right)\right],$$

$${\beta }_e\equiv \frac{8\pi n_{e^{-}}k{T}_e}{B_0^2}=\frac{2{\Theta }_e}{{\sigma }_{e,h}(2-q)\left(1+\frac{{\hat{\gamma }}_e}{{\hat{\gamma }}_e-1}{\Theta }_e\right)}\cdot$$

$${\omega }_p^2={\omega }_{p{e}^{-}}^2+{\omega }_{p{e}^{+}}^2+{\omega }_{pi}^2.$$

$${\omega }_{p{e}^{-}}^2$$

$${\omega }_{p{e}^{+}}^2$$

$${\omega }_{pi}^2={\omega }_{p{e}^{-}}^2\frac{m_e}{m_i}q\frac{1+\frac{{\Theta }_e}{{\hat{\gamma }}_e-1}}{1+\frac{{\Theta }_i}{{\hat{\gamma }}_i-1}}.$$

$${\rho }_{Le}\equiv \frac{{\sigma }_{e,c}{m}_e{c}^2}{e{B}_0}=\frac{c}{{\omega }_{p{e}^{-}}}{\left(\frac{{\sigma }_{e,c}}{2-q}\right)}^{1/2}{\left(1+\frac{{\Theta }_e}{{\hat{\gamma }}_e-1}\right)}^{-1/2}$$

$${\rho }_{Li}\equiv \frac{{\sigma }_{i,c}{m}_i{c}^2}{e{B}_0}=\frac{c}{{\omega }_{p{e}^{-}}}{\left(\frac{{\sigma }_{i,c}}{q}\frac{m_i}{m_e}\right)}^{1/2}{\left(1+\frac{{\Theta }_e}{{\hat{\gamma }}_e-1}\right)}^{-1/2}.$$

$$⟨{\gamma }_e-1⟩=a(\sigma ,\Theta ){\kappa }^{-\chi (\sigma ,\Theta )}+b(\sigma ,\Theta ),$$

Full text

active galaxies – magnetic reconnection – plasmas – acceleration of particles

††journal: ApJ††software: TRISTAN-MP (Buneman, 1993; Spitkovsky, 2005), IDL version 8.6 (Exelis Visual Information Solutions, Boulder, Colorado).

1 Introduction

A fundamental question in the physics of astrophysical relativistic outflows is how their energy, which is initially carried in the form of Poynting flux, is first transferred to the plasma, and then radiated away to power the observed emission. Magnetic field dissipation via reconnection has been often invoked to explain the non-thermal signatures of pulsar wind nebulae (PWNe; e.g., Lyubarsky & Kirk, 2001; Pétri & Lyubarsky, 2007; Sironi & Spitkovsky, 2011a; Cerutti et al., 2012; Philippov & Spitkovsky, 2014, see Sironi & Cerutti 2017 for a recent review), gamma-ray bursts (GRBs; e.g., Thompson, 1994; Usov, 1994; Spruit et al., 2001; Drenkhahn & Spruit, 2002; Lyutikov & Blandford, 2003; Giannios, 2008; Beniamini & Giannios, 2017), and jets from active galactic nuclei (AGN; e.g., Romanova & Lovelace, 1992; Giannios et al., 2009, 2010; Giannios, 2013; Petropoulou et al., 2016; Nalewajko et al., 2018; Christie et al., 2019).

In most relativistic astrophysical outflows, reconnection proceeds in the so-called relativistic regime in which the Alfvén  velocity of the plasma approaches the speed of light (or equivalently the plasma magnetization, defined as the ratio of magnetic to particle enthalpy densities, is $\sigma\gtrsim 1$). The physics of reconnection can only be captured from first principles by means of fully-kinetic particle-in-cell (PIC) simulations. Extensive numerical work on relativistic reconnection of electron-positron (pair) plasmas has been performed in two dimensions (2D; e.g., Zenitani & Hoshino, 2001, 2007; Daughton & Karimabadi, 2007; Cerutti et al., 2012; Sironi & Spitkovsky, 2014; Guo et al., 2014, 2015; Liu et al., 2015; Nalewajko et al., 2015; Sironi et al., 2015, 2016; Werner et al., 2016; Kagan et al., 2018; Petropoulou & Sironi, 2018; Hakobyan et al., 2018) and in three dimensions (3D; e.g., Zenitani & Hoshino, 2005, 2008; Liu et al., 2011; Sironi & Spitkovsky, 2011a, 2012; Kagan et al., 2013; Cerutti et al., 2014; Sironi & Spitkovsky, 2014; Guo et al., 2015; Werner & Uzdensky, 2017), whereas the study of trans-relativistic and relativistic reconnection in 2D electron-proton plasmas became possible more recently (e.g., Melzani et al., 2014; Sironi et al., 2015; Guo et al., 2016; Rowan et al., 2017; Werner et al., 2018; Ball et al., 2018).

In contrast to other astrophysical outflows, such as PWNe, the plasma composition of astrophysical jets is largely unknown. On the one hand, there is no direct way of probing the plasma composition in jets and, on the other hand, there are large theoretical uncertainties about the jet baryon loading mechanisms (for recent kinetic simulations of black-hole jet launching, see Parfrey et al., 2019). As a result, any attempts to infer the jet plasma composition rely on the modeling of the emitted radiation (e.g., Ghisellini, 2012; Ghisellini et al., 2014), which, however, suffers from degeneracies that are inherent in the radiative models. For AGN jets, in particular, both pair and electron-proton compositions have been discussed in the literature. A pure pair composition in powerful AGN jets (e.g., in flat spectrum radio quasars, FSRQs) is disfavored, since bulk Comptonization of the ambient low-energy photons by the pairs would result in luminous spectral features in X-rays that are not observed (e.g., Sikora et al. (1997); Sikora & Madejski (2000); see however Kammoun et al. (2018)). This argument does not apply to less powerful jets (such as BL Lac type sources), since the ambient radiation fields are weak or even absent, and a pure pair plasma cannot be excluded in this case. If jets are devoid of pairs, namely they are composed of electron-proton plasmas, the inferred power (which is dominated by the kinetic power of protons) is large, usually exceeding the accretion power (e.g. Ghisellini et al., 2014; Madejski et al., 2016). A mixed composition with tens of pairs per proton may be more realistic, as it can reduce the inferred jet power by a factor equal to the lepton-to-proton number ratio, the so-called pair multiplicity (e.g. Ghisellini et al., 2010; Ghisellini, 2012; Madejski et al., 2016). The presence of pairs in the dissipation regions of jets is also expected to affect the average energy per lepton available for particle heating as well as the efficiency with which non-thermal particles are accelerated.

The goal of this work is to study the general properties of relativistic reconnection in the unexplored regime of plasmas with mixed composition. We focus on electron-positron-proton (or pair-proton) plasmas, as they bridge the gap between the pair plasma and electron-proton plasma cases that have been extensively studied in the past. We perform a suite of large-scale 2D PIC simulations using the realistic proton-to-electron mass ratio ($m_{i}/m_{e}=1836$) while varying three physical parameters, namely the plasma magnetization ($\sigma=1,3$ and 10), the plasma temperature ($\Theta_{e}\equiv kT_{e}/m_{e}c^{2}=0.1-100$ with equal electron and proton temperatures), and the number of pairs per proton ($\kappa=1-199$). In this study, even in cases where the pairs dominate by number, the plasma rest mass energy is governed by protons. We study, for the first time, the inflows and outflows of plasma in the reconnection region, the energy partition between pairs, protons, and magnetic fields, and the energy distributions of accelerated particles as a function of the pair multiplicity.

This paper is organized as follows. In Sect. 2 we describe the setup of our simulations. In Sect. 3 we present the structure of the reconnection layer for different pair multiplicities. In Sect. 4 we focus on the inflow and outflow motions of the plasma and in Sect. 5 we discuss the energy partition between magnetic fields and different particle species in the reconnection region. In Sect. 6 we focus on the evolution of the particle energy spectrum, illustrating how the lepton power-law slope depends on the pair multiplicity. In Sect. 7 we discuss the astrophysical implications of our findings and conclude in Sect. 8 with a summary of our results. Readers interested primarily in the application of our results to jetted AGN can move directly to Sect. 7.

2 Numerical setup

We use the 3D electromagnetic PIC code TRISTAN-MP (Buneman, 1993; Spitkovsky, 2005) to study magnetic reconnection in pair-proton plasmas. We explore anti-parallel reconnection, i.e., we set the guide field perpendicular to the alternating fields to be zero. The reconnection layer is initialized as a Harris sheet of length $L$ with the magnetic field $\bm{B}=-B_{0}\,\tanh\,(2\pi y/\Delta)\,\bm{\hat{x}}$ reversing at $y=0$ over a thickness $\Delta$. Here, we set $\Delta$=80 $c/\omega_{p}$, where $\omega_{p}$ is the all-species plasma frequency defined in Eq. (A8), and choose a spatial resolution of $\,c/\omega_{p}$=3 computational cells.

The field strength $B_{0}$ is defined through the (total) plasma magnetization $\sigma=B_{0}^{2}/4\pi h$, where $h$ is the enthalpy density of the unreconnected plasma including all species (see Eq. (A1)). The Alfvén speed is related to the magnetization as $v_{A}/c=\sqrt{\sigma/(\sigma+1)}$. We focus on the regime of relativistic reconnection (i.e., $v_{A}/c\sim 1$) and explore cases with $\sigma=1,3$ and 10 (see Table 1). The proton and pair plasmas outside the layer are initialized with the same temperature ($T_{i}=T_{e}$). We consider cases where the pairs are initially relativistically hot ($\Theta_{e}\equiv kT_{e}/m_{e}c^{2}=1$, 10, and 100), but for completeness we study also a few cases with initially colder pairs ($\Theta_{e}=0.1$). In all simulations, the protons are non-relativistic ($\Theta_{i}\equiv kT_{i}/m_{i}c^{2}=\Theta_{e}m_{e}/m_{i}\ll 1$).

Let $N_{ppc}$ denote the total number of computational particles per cell, which is equally partitioned between negatively and positively charged particles. If $q=2/(\kappa+1)$ denotes the physical number ratio of protons to electrons in a plasma with pair multiplicity $\kappa$, then the number of computational protons and positrons per cell is given, respectively, by $(q/2)N_{ppc}$ and $[(1-q)/2]N_{ppc}$111We fill the cells with particles by performing two cycles of injection. We first inject protons and electrons at equal numbers (i.e., $N_{1}=(q/2)N_{ppc}$ per cell) and then inject positrons and electrons, with a number of $N_{2}=[(1-q)/2]N_{ppc}$ per cell for each component. The injection is not done on a cell-by-cell basis, but in slabs partitioned along the $y$ direction with $N_{cells}$ each, which are handled by different computer cores. When either $N_{1}\times N_{cells}$ or $N_{2}\times N_{cells}$ is $<10$, the actual number of particles injected is randomly drawn from a Poisson distribution with mean value equal to $N_{1}\times N_{cells}$ or $N_{2}\times N_{cells}$, respectively.. We varied $N_{ppc}$ from 4 to 64 and checked the convergence of our results in regard to the reconnection rate, outflow four-velocity, and particle energy distributions. For pair-proton simulations with high pair multiplicity (e.g., $\kappa>10$), we need to use $N_{ppc}>16$ to achieve convergence (within a few percent in inflow rate and outflow four-velocity), whereas for electron-proton simulations we find that 4 particles per cell are sufficient. For cases with high $\kappa$ (low $q$) there is a low probability of proton “injection” in a given cell due to the small (physical) fraction of protons per electron. This introduces an appreciable level of shot noise in fluid quantities that are computed from (or governed by) the protons (e.g., outflow four-velocity), which can be mitigated by increasing the number of computational particles per cell.

The magnetic pressure outside the current sheet is balanced by the particle pressure in the sheet. This is achieved by adding a component of hot plasma with the same composition as in the upstream region and over-density $\eta=3$ relative to the all-species number density outside the layer. We exclude the hot particles initialized in the current sheet from the particle energy spectra and from all thermodynamical quantities (except the plasma number density), as their properties depend on our choice of the sheet initialization.

Our simulations are performed in a 2D domain, but all three components of the velocity and of the electromagnetic fields are tracked. We adopt periodic boundary conditions in the $x$ direction of the reconnection outflow and we employ an expanding simulation box in the $y$ direction (i.e., the direction of the reconnection inflow). We also use two moving injectors receding from $y=0$ along $\pm{\bm{\hat{y}}}$, which constantly introduce fresh magnetized plasma into the simulation domain (for details, see Sironi & Spitkovsky, 2011b, 2014; Sironi et al., 2016; Rowan et al., 2017; Ball et al., 2018). In all cases, the box size along the $y$ direction increases over time and by the end of the simulation it is comparable or larger than the $x$ extent.

We trigger reconnection at the center of the simulation domain by instantaneously removing the pressure of hot particles that were initialized in the sheet (Sironi et al., 2016; Ball et al., 2018). This causes a central collapse of the current sheet and the formation of two “reconnection fronts” that are are pulled along the layer towards the edges of the box by the magnetic tension force and reach the boundaries at $t\sim L/2v_{A}$. The main advantage of this simulation setup is that the results are independent of the initialization of the sheet (i.e., overdensity, temperature, and thickness222This is true if the sheet is thick enough so that it does not become spontaneously tearing unstable at locations that have not been swept up yet by the receding reconnection fronts.) in contrast to the untriggered cases where the influence of the initial conditions may affect the temporal evolution of the reconnection rate and the particle energy distributions at early times (see, e.g., Fig. 4 of Petropoulou & Sironi, 2018).

We choose as our typical unit of length the Larmor radius of electrons ($\rho_{Le}$) with Lorentz factor equal to the cold pair plasma magnetization ($\sigma_{e,c}$), namely $\rho_{Le}=\sigma_{e,c}m_{e}c^{2}/eB_{0}$, implicitly assuming that reconnection transfers all the magnetic energy to relativistic pairs (for definitions, see Eq. (A12) and Eq. (A4)). The proton Larmor radius is defined in a similar way, i.e., $\rho_{Li}=\sigma_{i,c}m_{i}c^{2}/eB_{0}$, where $\sigma_{i,c}$ is the cold proton plasma magnetization (see Eq. (A6)). The size of the computational domain along the reconnection layer $L$ ranges from hundreds to thousands of $\rho_{Le}$ and tens to hundreds of $\rho_{Li}$ (see Table 1). The fact that the Larmor radii change as a function of pair multiplicity is a direct result of our choice to fix the total $\sigma$ and electron thermal spread $\Theta_{e}$, as shown in Fig. 16 of Appendix A.

A key parameter in our study, as it will become clear in the following sections, is the hot pair plasma magnetization. This is defined as $\sigma_{e,h}\equiv B_{0}^{2}/4\pi h_{\pm}$, where $h_{\pm}$ is the enthalpy density of the upstream pair plasma, and it relates to the total $\sigma$ as:

[TABLE]

where $q=2/(\kappa+1)$ is the ratio of proton-to-electron number densities and $\hat{\gamma}_{i,e}$ are the adiabatic indices of protons and leptons. Eq. (1) can be simplified in the following asymptotic regimes:

• •

relativistically cold electrons ($\Theta_{e}\ll 1$). Here, $\sigma_{e,h}\approx\sigma[m_{i}/m_{e}+\kappa]/\kappa$. For electron-proton plasmas (or in general, if $\kappa\ll m_{i}/m_{e}$) this reduces to the well-known result $\sigma_{e,h}\approx\sigma m_{i}/m_{e}$, whereas for pair-dominated plasmas with $\kappa\gg m_{i}/m_{e}$, we find $\sigma_{e,h}\approx\sigma$. Although pairs are cold, if their number density is sufficiently high, like in the latter case, their pressure (which is $\propto\kappa\Theta_{e}$) can be more important than the proton rest-mass energy density.

• •

relativistically hot electrons ($1<\Theta_{e}<m_{i}/m_{e}$). Here, $\sigma_{e,h}\approx\sigma[m_{i}/m_{e}+4\Theta_{e}\kappa]/4\Theta_{e}\kappa$. This reduces to $\sigma_{e,h}\approx\sigma m_{i}/(4m_{e}\Theta_{e})$ for $\kappa=1$, while for $\kappa\gg(m_{i}/m_{e})/4\Theta_{e}$ we find $\sigma_{e,h}\approx\sigma$. In the latter case, the pressure of the hot pairs is large enough to dominate over the rest-mass energy density of protons. Note that the critical pair multiplicity here is lower by a factor of $\sim 4\Theta_{e}$ compared to the cold electron case (see first bullet point).

• •

relativistically hot protons ($\Theta_{e}\gg m_{i}/m_{e}$). In this ultra-relativistic regime, all fundamental plasma scales (e.g., the plasma frequencies and skin depths) become independent of the particle rest mass. They depend only on the average particle energy which, in this regime, is similar for protons and pairs. Here, $\sigma_{e,h}\approx\sigma$ independent of $\kappa$.

In this study, we focus on cases where the protons are non-relativistic and dominate the mass density. We refer the reader to Appendix A, for the full list of parameters and their definitions.

3 Structure of the reconnection layer

3.1 Temporal evolution

To illustrate the temporal evolution of the reconnection region we show in Fig. 1 snapshots of the 2D structure of the particle number density from one of our simulations in a $\sigma=1,\Theta_{e}=1$ pair-proton plasma with pair multiplicity $\kappa=19$ (A2 in Table 1). The localized (at the center) removal of pressure from the hot particle population initialized in the sheet (see Sect. 2) causes its collapse, thus leading to the formation of a central (or primary) X-point. Two reconnection fronts form on opposite sides of the primary X-point and move outwards due to the tension of the magnetic field lines. Plasmoid and secondary X-point formation takes place in the low-density region between the moving fronts, as shown in panels (a) and (b). The fronts reach the boundaries of the simulation domain at $t\approx L/2v_{A}\approx 0.7L/c$ and form the so-called boundary island, whose size eventually becomes a significant fraction of the layer length (here, $\sim 0.4\,L$ as shown in panels d and e). The formation of such a large plasmoid, which is the result of periodic boundary conditions, will eventually inhibit the inflow of fresh plasma into the layer, thus shutting off the reconnection process. We verified that the reconnection process remains active333We characterize the reconnection process as active, as long as the inflow rate of plasma into the reconnection region does not show a monotonically decreasing trend with time and remains $\gtrsim 0.01v_{A}$ at all times. for the entire duration of all simulations listed in Table 1 except A3, C6, and E3.

3.2 Dependence on pair multiplicity

The effect that the pair multiplicity $\kappa$ has on the appearance of the reconnection region is illustrated in Fig. 2, where 2D snapshots of the all-species particle density, including particles initially present in the sheet, are plotted for increasing values of $\kappa$ (top to bottom) in plasmas with $\sigma=1,\Theta_{e}=1$ (left panel) and $\sigma=1,\Theta_{e}=10$ (right panel). In cases with fixed $\sigma$ and $\Theta_{e}$ but increasing $\kappa$ we find that the plasma outflows along the layer become more uniform (i.e., fewer X-points and plasmoids form in the layer) and the typical size of the plasmoids decreases. For fixed $\sigma$ and $\kappa$, an increasing upstream plasma temperature also leads to smaller plasmoids and less fragmentation in the reconnection region (compare left and right panels in Fig. 2).

One might argue that the differences in the appearance of the layer as a function of $\kappa$ are merely a result of the different box sizes in terms of the proton skin depth or alternatively $\rho_{Li}$ (see Table 1). To check this possibility, we compare cases with different physical conditions, but similar box sizes in terms of $\rho_{Li}$. We find that the plasma conditions have a major effect on the appearance of the layer (for details, see Appendix B) and that the differences seen in Fig. 2 are not just a numerical artifact.

Empirically, we find that the most important parameter controlling the appearance of the layer turns out to be $\sigma_{e,h}$. We find that the layer structure is similar for different values of the pair multiplicity and temperature, as long as $\sigma_{e,h}$ is nearly the same. For example, compare panel (e) on the left side to panel (c) on the right side of Fig. 2. Typically, the density profile is smoother and the plasmoid sizes are smaller for lower $\sigma_{e,h}$ values444 The apparent correlation of the plasmoid size with $\sigma_{e,h}$ is likely related to the dependence of the electron Larmor radius on $\sigma_{e,h}$ (i.e., $\rho_{Le}\propto\sigma_{e,h}^{1/2}$). (e.g., compare panels (a) and (f) on the left side of Fig. 2).

Similar results have been presented by Ball et al. (2018) (see Fig. 4 therein) for trans-relativistic electron-proton reconnection and an increasing electron plasma $\beta_{e}$, defined as the ratio of upstream electron plasma pressure and magnetic pressure (see Eq. (A7)). The similarity of our findings is not unexpected and can be understood as follows. The increasing pair multiplicity corresponds to a decreasing hot pair plasma magnetization $\sigma_{e,h}$ (see Eq. (1) and Fig. 15), which in turn is inversely proportional to $\beta_{e}$ in the limit of $\kappa\gg 1$ (see Eq. (A7)).

Henceforth, we choose $\sigma_{e,h}$ over $\beta_{e}$ to perform our parameter study, since the relative contribution of the rest-mass and internal energy densities to the enthalpy density of the upstream plasma varies among our simulations. In Sections 5 and 6 we will also demonstrate that $\sigma_{e,h}$ is the main parameter that regulates the energy partition and the power-law slope of the lepton energy spectrum.

4 Inflows and outflows

To compute the reconnection rate in our simulations, we average at each time the inflow speed over a slab centered at $x=0.5\,L$ with width $0.2\,L$ across the layer (i.e., along the $y$ direction) and length $0.5\,L$ (along the $x$ direction). Our results are nearly insensitive to the choice of the slab dimensions as long as the region occupied by the boundary island, where the inflow rate is inhibited, is excluded from the averaging process. The spatially averaged inflow rate is then averaged over time for $t>L/2v_{A}$, i.e., excluding times when the reconnection fronts are still in the slab.

Our results for simulations with different $\sigma,\,\Theta_{e}$, and $\kappa$ (Table 1) are presented in Fig. 3, where the average inflow speed $v_{in}$ (normalized to $v_{A}$) is plotted as a function of the hot pair plasma magnetization $\sigma_{e,h}$. Results for pair-proton and electron-proton cases are indicated with filled and open symbols, respectively. The error bars, which indicate the standard deviation of the reconnection rate over the duration of the simulation, become typically larger with increasing $\sigma_{e,h}$ and fixed $\sigma,\Theta_{e}$. This suggests that the layer becomes accordingly more structured (see also Fig. 2), since the temporal variations of the reconnection rate about its average value relate to the motion and coalescence of plasmoids (see also Petropoulou & Sironi, 2018). We find a weak dependence of the average reconnection rate on $\sigma_{e,h}$, as this changes only by a factor of $\sim 3$ ($\sim 0.05-0.15$) over more than three orders of magnitude in $\sigma_{e,h}$. Despite this weak dependence, our results reveal a clear trend of lower reconnection rates at lower $\sigma_{e,h}$ (i.e., at higher $\beta_{e}$), in agreement with the findings of Ball et al. (2018).

The four-velocity of the plasma outflows in the reconnection region along the $x$ direction, $\Gamma v_{out}/c$, is computed using all particle species, although it is controlled by the protons that contribute most to the plasma inertia. The Lorentz factor $\Gamma$ takes into account the motion in all three directions, but the bulk motion along $x$ dominates. To estimate the maximum four-velocity we compute at each time the 95th percentile555We compute the absolute values of the four-velocity measured at different locations along the layer at $y=0$, sort them in descending order, and determine the value below which 95% of the measurements falls. of all values of $\Gamma v_{out}/c$ at $y=0$, and show in Fig. 4 its temporal evolution from simulations with $\sigma=1,\Theta_{e}=1$, and different pair multiplicities. The outflowing plasma accelerates soon after the onset of reconnection, its motion becomes relativistic, and its maximum four-velocity approaches the asymptotic value $\sqrt{\sigma}$ (Lyubarsky, 2005). We note that the 95th percentile of $\Gamma v_{out}/c$ values in the layer provides a more conservative estimate of the maximum outflow four-velocity than the one derived using, for example, the fifth (or tenth) largest value (see e.g. Sironi et al., 2016). We verified that with the latter method the peak four-velocity is even closer to $\sqrt{\sigma}$. We find no systematic dependence of the maximum outflow four-velocity on the pair multiplicity, apart from the fact that the bulk acceleration is more gradual in plasmas with $\kappa=199$ (see black line in Fig. 4); this is also true for other values of $\Theta_{e}$ and $\sigma=1-3$.

5 Energy partition in the reconnection region

The question of how the available energy is shared between particles and magnetic fields in the region where plasma has undergone reconnection (henceforth, the reconnection region) is of particular astrophysical importance, since it is related to the intensity and spectrum of the associated electromagnetic radiation. Here, we study the energy partition in pair-proton plasmas post reconnection, as a function of pair multiplicity, magnetization, and temperature of the unreconnected plasma.

To identify the reconnection region we use a mixing criterion, as proposed by Daughton et al. (2014). Particles are tagged with an identifier (0 or 1) based on their initial location (below or above) with respect to the current sheet. Particles from these two regions get mixed in the course of the reconnection process. We identify the reconnection region by the ensemble of computational cells with mixing fraction above a certain threshold $\epsilon$ and below $1-\epsilon$; here, we employed $\epsilon=0.01$666We verified that our results are insensitive to the exact value, except for very early times (i.e., $\lesssim 0.15L/c$) where the small size of the reconnection region makes the computation of quantities therein sensitive to the choice of $\epsilon$. (for more details, we refer the reader to Rowan et al., 2017; Ball et al., 2018). 2D snapshots of the mixing fraction from two indicative simulations (see A1 and A6 in Table 1) are presented in Fig. 5, where the reconnection region is identified by the mixed colors (green and red).

We compute the kinetic energy of each particle species by summing up the contributions from all computational cells that define the reconnection region, namely $u_{j}=m_{j}c^{2}\sum\limits_{cells}n_{j}(\gamma_{j}-1)$, where $\gamma_{j}$ is the average Lorentz factor of particles of species $j$ in a computational cell. We then normalize $u_{j}$ to the total energy $u_{tot}=u_{B}+\!\!\!\sum\limits_{j=i,e^{\pm}}\!\!\!u_{j}$, where $u_{B}=\sum\limits_{cells}B^{2}/8\pi$. In Fig. 6 we show the temporal evolution of $u_{j}/u_{tot}$ for the same cases as those shown in Fig. 5. At very early times, when the reconnection region is small (see grey-colored region in Fig. 6), the plasma properties therein depend on how exactly the reconnection region is identified. Yet, neither the time-averaged properties nor their late-time evolution are sensitive to the definition of the reconnection region. Given that there might be also other factors affecting the early time evolution (e.g., initial setup), we henceforth ignore this transitional early period. At later times, the ratio of post-reconnection magnetic energy to the total energy decreases gradually with time, whereas the pair energy density ratio reaches an almost constant value very soon after the onset of reconnection (i.e., already at $0.4L/c$). The proton energy ratio asymptotes to a constant value typically at later times compared to the pairs, but our simulations are long enough to capture the steady-state values of all energy ratios. We find similar temporal trends for other cases as well.

The time-averaged energy ratios of protons, pairs, and magnetic fields in the reconnection region are presented in Fig. 7. The leftmost point in each series with a given color corresponds to pair-proton plasmas with $\kappa=199$ and the rightmost point corresponds to the pure electron-proton case with $\kappa=1$ (open symbols). The fraction of energy that remains in the post-reconnection magnetic field is $\sim 1/3$ and is approximately constant for a wide range of $\sigma_{e,h}$ values, spanning almost three orders of magnitude (bottom panel). Only for $\sigma_{e,h}<3$, we find sub-equipartition values, i.e., $u_{B}/u_{tot}<1/3$. In this parameter regime, the pairs in the plasma carry most of the upstream total energy. Upon entering the reconnection region, the pair kinetic energy increases even further at the expense of magnetic energy due to field dissipation. As a result, the post-reconnection magnetic energy for $\sigma_{e,h}<3$ is only a small fraction of the total energy ($u_{B}/u_{tot}\sim 0.1-0.2$).

One can empirically define two regimes of interest for the particle energy ratios: a low-$\sigma_{e,h}$ regime ($\sigma_{e,h}\lesssim 30$), where $u_{i}/u_{tot}\propto\sigma_{e,h}$ and $u_{e^{\pm}}/u_{tot}\propto\sigma_{e,h}^{-1/2}$, and a high-$\sigma_{e,h}$ regime ($\sigma_{e,h}>30$), where both ratios are almost independent of the hot pair plasma magnetization. In both regimes, there is no dependence of the particle energy ratios on $\Theta_{e}$, but a weak dependence on the total plasma magnetization $\sigma$ is evident. This can be more clearly seen in the middle panel of Fig. 7, where points with the lowest $\sigma$ (black and cyan symbols) systematically lie below points with higher $\sigma$. Finally, energy equipartition between magnetic fields, protons, and pairs is asymptotically achieved for $\sigma\gg 1$ and $\sigma_{e,h}\gtrsim 30$, with each component carrying $\sim 1/3$ of the total energy.

The dependence of the particle energy densities on $\sigma_{e,h}$ could originate from either changes in the number density or in the mean particle Lorentz factor, or both. A proxy of the average post-reconnection particle Lorentz factor, $\langle\gamma_{j}-1\rangle=\sum\limits_{cells}u_{j}/\sum\limits_{cells}n_{j}m_{j}c^{2}$, is plotted against $\sigma_{e,h}$ in the left panel Fig. 8 for protons (top panel) and pairs (bottom panel). In all cases, we find that the post-reconnection mean proton Lorentz factor is almost independent of $\sigma_{e,h}$ and $\Theta_{e}$, but has a dependence on $\sigma$, with larger values leading to higher mean proton Lorentz factors. Indeed, when $\langle\gamma_{i}-1\rangle$ is normalized to $\alpha_{i}\sigma$ (with $\alpha_{i}=1/3$ for $\sigma=1,3$ and $\alpha_{i}=1/5$ for $\sigma=10$) all curves coincide, as shown in the right panel of Fig. 8. In contrast to the protons, the mean lepton Lorentz factor depends on $\Theta_{e}$, $\sigma$, and $\sigma_{e,h}$, as shown in the left panel of Fig. 8. We empirically find for $\sigma\geq 1$ that the mean lepton Lorentz factor can be approximated as (see also right panel in Fig. 8):

[TABLE]

The asymptotic value of the mean lepton Lorentz factor for $\sigma_{e,h}\ll 30$ implies that, in this regime, the pairs in the reconnection region still bear memory of their initial (pre-reconnection) conditions (and, in particular, of $\Theta_{e}$), in agreement with the discussion on Fig. 7. In the high-$\sigma_{e,h}$ regime, the mean lepton Lorentz factor scales almost linearly with $\sigma_{e,h}$. This asymptotic behavior of $\langle\gamma_{e}-1\rangle$ can be understood as follows. For fixed $\sigma$ and $\Theta_{e}$ (i.e., fixed amount of post-reconnection energy available for the particles), the energy per lepton increases as the number of leptons per proton decreases, or equivalently, as $\sigma_{e,h}$ increases (see also Eq. (1)). We refer the reader to Appendix D, for a quantitative discussion on the dependence of the mean lepton Lorentz factor on the physical parameters $\sigma,\Theta_{e}$, and $\kappa$ of the upstream plasma.

6 Particle energy distributions

After having discussed the general properties of reconnection in pair-proton plasmas (see sections 3-5), we continue our study by examining the particle energy distributions and their dependence on physical parameters, most notably on $\sigma_{e,h}$.

6.1 Temporal evolution of particle energy spectra

The energy distribution of each particle species is defined as $f_{j}(E)\equiv dN_{j}/dE$, where $E$ is the particle kinetic energy and $j=i,e^{-},e^{+}$. Henceforth, all particle energies are kinetic (i.e., excluding rest mass), unless stated otherwise.

As a representative example, we present in Fig. 9 the temporal evolution of the electron, positron, and proton energy distributions (from top to bottom) from a simulation with $\sigma=1,\Theta_{e}=1$, and $\kappa=19$ (see also Fig. 1, for a depiction of the layer structure). The energy distributions of each particle species are normalized to the total number of particles of that species in the reconnection region at the end of the simulation. The displayed spectra exclude the particle population that was initialized in the current sheet. For reference, the spectrum obtained at the time the reconnection fronts reach the boundaries (i.e., $t=L/2v_{A}$) is shown with a dashed black line.

Soon after the onset of reconnection, the electron and positron energy spectra in the reconnection region begin to deviate from their initial Maxwell-Jüttner distributions. They develop a non-thermal component even before the time the reconnection fronts reach the boundaries of the layer (i.e., at $ct/L\sim 0.7$). The non-thermal part of the spectrum of pairs can be described by a power law above a characteristic energy where the post-reconnection energy spectrum $Ef_{j}(E)$ obtains its peak value. For the adopted parameters, we find $E_{pk,e}/m_{e}c^{2}\sim 10$ in agreement with the value of the mean post-reconnection Lorentz factor that we derived in Sect. 5 (see third black symbol from the left in bottom panel of Fig. 8).

There is a clear difference between the temporal evolution of the lepton and proton energy distributions. More specifically, the non-thermal component of the proton spectrum begins to emerge only at $t>L/2v_{A}$, after the fronts have reached the boundaries. At earlier times, the proton energy spectrum shows a narrow peak that evolves with time. We interpret this early-time spectral feature as a result of heating and bulk motion of the proton plasma, whose outflow four-velocity evolves strongly for $t<0.5L/c$ (see blue curve in the top panel of Fig. 4). Similar results were obtained by Ball et al. (2018) for trans-relativistic reconnection in electron-proton plasmas (see Fig. 3 therein).

The late-time development of the power law in the proton distribution can be understood in terms of the interactions of particles with various structures in the layer. X-points are typically smaller than the proton Larmor radius, so direct proton acceleration by the non-ideal reconnection electric field is not very efficient. We find evidence of proton acceleration only when the boundary island, which is the biggest structure in the layer, begins to form. We argue that in a much larger simulation domain, where bigger secondary plasmoids could form, protons should show signs of acceleration even before the reconnection fronts interact with the boundaries.

6.2 Effects of pair multiplicity

To illustrate the dependence of the particle energy distributions on pair multiplicity, we show in Fig. 10 the energy spectra from a set of simulations with $\sigma=1$, $\Theta_{e}=1$ and different values of $\kappa$ marked on the plots. Thick solid and thin dashed lines show the spectra from the reconnection region and the whole simulation domain, respectively. The spectra are computed at the end of each simulation and are normalized to the total number of protons within the reconnection region. The vertical dash-dotted line in each panel marks the energy of particles with Larmor radius777The Larmor radius is computed using the upstream magnetic field strength. $0.1\,L$, i.e., comparable to the size of the largest plasmoids in the layer.

The peak energy of the pair energy distributions depends strongly on the pair multiplicity for $\kappa<6$, and becomes approximately constant (here, $E_{pk,e}/m_{e}c^{2}\sim 10$) for higher pair multiplicities. On the contrary, the peak proton energy is approximately constant for all $\kappa$ values we explored. The dependence of the peak particle energy on $\kappa$ is more clearly illustrated in Fig. 11, where the energy distributions of each particle species are plotted for different values of $\kappa$. These findings are in agreement with those presented in Fig. 8 for the mean particle Lorentz factor. The fact that the mean and the peak lepton energies are comparable is not surprising; most of the energy is expected to reside at the peak of the energy distribution, given that the power-law slopes of the lepton energy spectra are typically $\gtrsim 2$ (see below and Sect. 6.3).

Above the peak energy $E_{pk,e}$, the pair energy spectra can be approximated by a power law with slope $p$ (i.e., $f(E)\propto E^{-p}$) followed by a cutoff. The power-law segment used for the estimation of the slope (see Sect. 6.3) is overplotted (dash-dotted blue lines) for guiding the eye. Inspection of the figure (see also Fig. 11) shows that the power law of the pair distributions becomes steeper (i.e., larger $p$ values) as the pair multiplicity increases (for details, see Sect. 6.3). The power law of the pair distributions extends well beyond their peak energy for all the cases we explored, except for the cases with the highest $\Theta_{e}$, which are discussed in Appendix E. For protons a well-developed power law forms only for small pair multiplicities (here, for $\kappa<19$), while their energy distribution shows a steep drop above the peak energy $E_{pk,i}/m_{e}c^{2}\sim 10^{3}$ for $\kappa=66$ and 199. This should not be mistakenly interpreted as a limitation of reconnection in accelerating protons in plasmas with high pair multiplicities. It is merely a result of the limited size of the computational domain in terms of the proton Larmor radius: $L/\rho_{Li}$ drops by a factor of ten between the simulations with $\kappa=3$ and $\kappa=199$, as shown in Table 1 (the dependence of the particle energy distributions on the box size is discussed in Appendix C). For these reasons, we do not attempt to study the spectral properties of the proton energy distributions and, in what follows, we focus on the energy distributions of pairs.

6.3 Power-law slope of pair energy spectra

We compute the slope of the power-law segment of the pair energy distributions and explore its dependence on the physical parameters. Due to the similarity between the energy distributions of positrons and electrons (see Fig. 10) it is sufficient to use one of the two for computing the slope. Henceforth, we use for this purpose the electron energy spectrum obtained at the end of each simulation.

The electron energy distribution can be generally described by two components: a low-energy broad component that forms due to heating and a high-energy component, which can be described as a truncated power law at low energies with an exponential cutoff at higher energies (see e.g., panels in middle row of Fig. 10). A detailed fit to the simulation data is very challenging due to the degeneracy in the model parameters describing the two components. For example, the choice of the low-energy end of the power law affects the broadness and normalization of the low-energy component and vice versa. The slope inferred from the two-component fit to the data can vary at most by $0.2$ depending on the other model parameters. Given the inherit uncertainties in the fitting procedure, in what follows, we identify the power-law segment by eye and fit it with a single power law (see dash-dotted blue lines in Fig. 10).

The extent of the power law is, in most cases, sufficient to allow a reliable estimation of its slope. We assign a systematic error of $\pm 0.2$ to the derived slope, which dominates the statistical error from the fits, to account for the subjective choice of the fitting energy range. For simulations with duration much larger than all others (see electron-proton cases in Table 1), we computed the slope also at earlier times (i.e., comparable to the duration of all other cases) and found no difference in the inferred $p$ value within the systematic error. Although a hard power law can be safely distinguished from the thermal part of the energy distribution, for very steep power laws with $p\gtrsim 4$, we cannot exclude the possibility that what we are identifying as a power law is in fact the tail of a thermal-like distribution or a multi-temperature distribution (see e.g., bottom right panel in Fig. 10). Detailed modeling of the energy distributions which is important for determining the temporal evolution of the cutoff energy or the shape of the exponential cutoff (Werner et al., 2016; Kagan et al., 2018; Petropoulou & Sironi, 2018) lies beyond the scope of this paper.

\startlongtable

Our results are summarized in Fig. 12 and Table 2, where the slope of the electron energy distribution $p$ is plotted as a function of $\sigma_{e,h}$. We do not include the results from runs H1-H3 with the highest plasma temperature (see Table 1), since the energy spectra are qualitatively different from all other cases (for details, see Appendix E and Ball et al. (2018)). The inferred power-law slopes fall onto two branches (dashed grey lines) that track each other for $\sigma_{e,h}\sim 30-300$, but merge in the asymptotic regime of $\sigma_{e,h}\gtrsim 10^{3}$, where both protons and pairs start to behave as one particle species (i.e., their Larmor radii become similar). The upper branch (i.e., larger $p$ values) is composed of results from $\kappa=1$ simulations, whereas results for larger multiplicities ($\kappa\gtrsim 6$) fall onto the lower branch (i.e., smaller $p$ values). For a fixed pair of $\Theta_{e}$ and $\sigma$ values, a transition from the lower to the upper branch, which is accompanied by a steepening of the power law, occurs at $\kappa\sim 3-6$. No transition is found for $\sigma=10$. The power-law slopes derived for the majority of the simulations lie on the lower branch for a wide range of $\sigma_{e,h}$ values, spanning more than three orders in magnitude, despite the differences in the total plasma magnetization, temperature, and pair multiplicity. This suggests that $\sigma_{e,h}$ is a key physical parameter in regard to the pair energy distribution.

In general, higher $\sigma_{e,h}$ lead to the production of harder power laws (i.e., smaller $p$ values), which is similar to the trend reported by Ball et al. (2018) for a decreasing electron plasma $\beta_{e}$ in electron-proton reconnection (see Fig. 13 therein). By tracking a large number of particles, Ball et al. (2018) showed that at low $\beta_{e}$ particles primarily accelerate by the non-ideal electric field at X-points. Since their number was found to decrease with increasing $\beta_{e}$, the authors argued that lower acceleration efficiencies and steeper power laws are expected at high $\beta_{e}$. The dependence of our derived power-law slopes on $\sigma_{e,h}$ can be qualitatively understood in the same context, since at high $\sigma_{e,h}$ (or equivalently low $\beta_{e}$) more X-points and secondary plasmoids are formed (see Sect. 3.2). A quantitative description of our results requires a detailed study of the electron acceleration, which is beyond the scope of this paper.

7 Astrophysical implications

In this section, we discuss the findings of our simulations in the context of AGN jets. We focus on blazars, the most extreme subclass of AGN, with jets closely aligned to our line of sight. The blazar jet emission has a characteristic double-humped shape with a broad low-energy component extending from radio wavelengths up to UV or X-ray energies, and a high-energy component extending across the X-ray and $\gamma$-ray bands (Ulrich et al., 1997; Fossati et al., 1998; Costamante et al., 2001). The low-energy hump is believed to be produced by synchrotron emission of relativistic pairs with a power-law (or broken power-law) energy distribution (e.g. Celotti & Ghisellini, 2008), which is suggestive of non-thermal particle acceleration. The synchrotron-emitting pairs can also inverse Compton scatter low-energy photons to $\gamma$-ray energies, which can explain the high-energy component of the blazar spectrum888This is true in leptonic scenarios where the broadband jet emission is attributed to relativistic pairs. This is our working hypothesis and our results should be interpreted in this framework.. In blazars with TeV $\gamma$-ray emission, electrons should accelerate up to Lorentz factors $10^{5}-10^{6}$ to explain the highest photon energies (e.g. Aleksić, 2012; Ahnen et al., 2018).

7.1 Properties of radiating particles

A key parameter in blazar emission models is the shape of the non-thermal pair distribution (e.g., power law, broken power law, log-parabolic, and others). The assumed distribution in most cases is phenomenological, as it is not derived from a physical scenario. Upon adopting a specific model for the energy distribution of accelerated pairs, its properties (e.g., power-law slope, minimum, and maximum Lorentz factors) are inferred by modeling the broadband blazar photon spectrum (e.g., Celotti & Ghisellini, 2008; Ghisellini et al., 2014). However, not all the model parameters can be uniquely determined due to degeneracies that are inherent in the radiative models (e.g., Cerruti et al., 2013).

Bearing in mind the aforementioned caveats, we continue with a tentative comparison of our results (see Sect. 5-6) with those inferred by radiative leptonic models. As an indicative example, we use the results of Celotti & Ghisellini (2008). The accelerated lepton distribution that was used for the modeling was assumed to be a broken power-law:

[TABLE]

where $s_{1}=1$, and $s_{2}$, $\gamma_{inj}$ were determined by the fit to the data. There is some degeneracy in the low-energy index, since distributions with even flatter spectra than the one above (i.e., $s_{1}<1$) cannot be usually distinguished by the data (see also Ghisellini et al., 2014).

In our simulations, we find that the post-reconnection pair energy distributions exhibit a power law extending well beyond a broad thermal-like component that peaks at $E_{pk,e}$ (see e.g., Figs. 10 and 11). At $E<E_{pk,e}$, the pair spectra in the reconnection region generally follow the low-energy tail of a Maxwell-Jüttner distribution (see e.g., Fig. 23), which can be modeled by an inverted power law (i.e., $s_{1}<0$). For the purposes of making a general comparison to the modeling results, we can phenomenologically describe the lepton energy spectra from our simulations by Eq. (5), with $s_{1}<1$, $s_{2}=p$, and a peak Lorentz factor $\gamma_{inj}=1+E_{pk,e}/m_{e}c^{2}$, which depends on the total magnetization and temperature of the plasma (see e.g., Fig. 10 and Fig. 23). Using the fitting results of Celotti & Ghisellini (2008) (see Table A1 therein), we compute the mean Lorentz factor of the accelerated distribution (i.e., without radiative cooling) and compare it against the one determined by our simulations (see e.g., Fig. 8).

Our results are summarized in Fig. 13, where the power-law index $p$ above the peak Lorentz factor of the distribution is plotted against the mean Lorentz factor $\langle\gamma_{e}\rangle$ of the distribution (for a tabulated list of our results, see Table 2). Open and filled triangles indicate the values from the leptonic modeling of Celotti & Ghisellini (2008) for FSRQs and BL Lac objects, respectively. The predictions of reconnection are shown with colored symbols (for details, see figure caption). The degeneracy of the power-law index $p$ on the physical parameters, such as $\sigma$ and $\Theta_{e}$ shown in Fig. 12, is lifted when $p$ is plotted against the mean lepton Lorentz factor. This is illustrated in Fig. 13, where, for fixed $\sigma$, curves corresponding to higher $\Theta_{e}$ values are shifted towards larger $\langle\gamma_{e}\rangle$ and $p$ values (upper right corner of the plot). For fixed plasma temperature but increasing $\sigma$, the curves are shifted towards lower $p$ values (i.e., harder power laws) and larger mean particle energies, regardless of the pair multiplicity.

Interestingly, the values from our simulations fall in the same range with those inferred by leptonic radiation models. More specifically, the numerically obtained curves for $\Theta_{e}=1$ and 10 enclose most of the results for FSRQs (open triangles). One can envision different families of curves that pass through the data points for FSRQs, which can be obtained by simply changing the temperature of the upstream plasma from $\Theta_{e}=1$ to $10$. For example, some FSRQ results could be interpreted by reconnection in pair-proton plasmas with $\Theta_{e}=3$, $\sigma=1$, and $\kappa\sim 1-10$ (imagine the blue line with circles shifted to the right and upwards). The relevant range of multiplicities would be somewhere between $\sim 10-70$, for $\Theta_{e}=3$ and $\sigma=3$ (imagine the red line with circles shifted to the right and upwards). We find that reconnection in cold pair-proton plasmas ($\Theta_{e}\ll 1$) with $\sigma\leq 10$ typically results in slopes and mean lepton energies that are not compatible with the FSRQ results.

BL Lac sources with $\langle\gamma_{e}\rangle\sim 10^{2}-10^{3}$ are compatible with our simulation results for reconnection in pair-proton plasmas with $\sigma\sim 3-10$, $\Theta_{e}\sim 10$ and $\kappa\sim 1-10$. The majority of BL Lac sources, however, requires mean Lorentz factors $>10^{3}$. Reconnection in strongly magnetized plasmas ($\sigma>10$) can lead to high values of the mean Lorentz factor, but at the same time produces hard power laws ($p<2$) above the peak Lorentz factor $\gamma_{inj}$ (see e.g., orange curves) that do not agree with the fitting results for $\langle\gamma_{e}\rangle\gg 10^{3}$ (filled triangles). In this regime, however, we argue that $\gamma_{inj}$ could be interpreted as the maximum Lorentz factor of a hard power law with $p<2$, as found in our high-$\sigma$ models, with $p$ now corresponding to the index $s_{1}$ (see Eq. (5)). Because the determination of the maximum Lorentz factor from the simulation spectra is not trivial (see e.g. Werner et al., 2018), we refrain from drawing strong conclusions from the comparison of our results to the BL Lac sources in the sample of Celotti & Ghisellini (2008).

7.2 Equipartition conditions

One of the reasons that makes the principle of energy equipartition between particles and magnetic fields attractive is that it leads to minimum power solutions for blazar jets (e.g., Dermer et al., 2014; Petropoulou et al., 2016). The energy density ratio of radiating particles and magnetic fields in the blazar emitting region is usually a free parameter, which is determined by the fitting of photon spectra. Leptonic emission models typically find $0.03\lesssim u_{e^{\pm}}/u_{B}\lesssim 30$, although specific sources may require even higher values (e.g., Celotti & Ghisellini, 2008; Tavecchio et al., 2010; Ghisellini et al., 2014). Alternatively, one can impose the constraint of rough energy equipartition between pairs and magnetic fields while searching for the best fit model, as demonstrated successfully by Cerutti et al. (2014); Dermer et al. (2014, 2015).

The post-reconnection ratio $u_{e^{\pm}}/u_{B}$ obtained from our simulations is plotted in Fig. 14 as a function of the mean lepton Lorentz factor $\langle\gamma_{e}\rangle$. We find that $0.2\lesssim u_{e^{\pm}}/u_{B}\lesssim 10$, with higher values obtained for hotter upstream plasmas. Even larger ratios, as those inferred by modeling of TeV BL Lacs (e.g., Tavecchio et al., 2010), would require a pool of ultra-relativistically hot particles entering the reconnection region. The presence of a guide field (i.e., of a magnetic field component that does not reconnect) would make the reconnection region more magnetically dominated, thus leading to $u_{e^{\pm}}/u_{B}<0.2$. More specifically, for electron-proton reconnection it was demonstrated that the fraction of magnetic energy transferred to non-thermal electrons can decrease from $\sim 50\%$ (in the absence of guide field) to $\sim 10\%$ for a guide field with strength comparable to that of the reconnecting field component (Sironi et al., 2015; Werner & Uzdensky, 2017). Yet, dissipation efficiencies as low as a few percent are still compatible with the global energetic requirements for AGN emission (Ghisellini et al., 2014; Sironi et al., 2015). A systematic study of the effects of the guide-field in pair-proton reconnection will be the topic of a future study.

8 Summary

For the first time, we have investigated magnetic reconnection in electron-positron-proton plasmas with a suite of large-scale 2D PIC simulations, covering a wide range of pair multiplicities ($\kappa=1-199$) for different values of the all-species plasma magnetization ($\sigma=1,3$ and 10) and plasma temperature ($\Theta_{e}=0.1,1,10$, and 100). In all cases we explored, protons in the upstream plasma have non relativistic temperatures and dominate the total mass.

The inflow rate of plasma into the reconnection region (i.e., the reconnection rate) ranges between $\sim 0.05v_{A}$ and $0.15v_{A}$ for a wide range of values of the hot pair plasma magnetization $\sigma_{e,h}$, with a weak trend towards higher rates for larger $\sigma_{e,h}$ values. The motion of the plasma outflow in the reconnection region, which is governed by the proton inertia, is relativistic with a maximum four-velocity that approaches the expected asymptotic value of $\sqrt{\sigma}$. We found no significant dependence of the outflow four-velocity on the pair multiplicity or temperature.

We showed that $\sim 1/3$ of the total energy remains in the post-reconnection magnetic field for $\sigma_{e,h}\gtrsim 3$, with the remaining $2/3$ of the energy being shared between pairs and protons. Energy equipartition between protons and pairs is achieved for $\sigma\gg 1$ and $\sigma_{e,h}\gtrsim 30$. For $\sigma_{e,h}\lesssim 3$, most of the energy in the reconnection region is carried by the pairs, with protons and magnetic fields contributing $\sim(1-10)\%$ to the total energy.

The reconnection process produces non-thermal particle energy distributions. We found that the mean Lorentz factor of the proton distribution (or, more accurately $\langle\gamma_{i}-1\rangle$) is almost independent of the pair multiplicity and plasma temperature, but it is approximately equal to $\sigma/3$. The mean Lorentz factor of the pair distribution can be described by a simple analytical expression (see Eq. (2)) for different values of $\sigma_{e,h}$, $\sigma$, and $\Theta_{e}$.

The electron and positron energy distributions in the reconnection region are similar and can be modeled as a power law with slope $p$ above a peak Lorentz factor, which, in most cases, is comparable with the mean Lorentz factor given by Eq. (2). The energy distribution below the peak can be, in general, approximated by a flat power law (with index $<0$). We showed that $p$ is mainly controlled by $\sigma_{e,h}$ (with harder power laws obtained for higher magnetizations) for a wide range of $\sigma,\Theta_{e}$, and $\kappa$ values. There is, however, a dependence of $p$ on pair multiplicity, with power laws getting steeper as $\kappa$ decreases from a few to unity.

We discussed the implications of our results in the context of AGN jets. We showed that reconnection in pair-proton plasmas naturally produces power-law pair distributions with slopes and average Lorentz factors similar to those obtained by leptonic modeling of the broadband jet emission. In general, we find that the majority of the modeling results can be explained in the context of reconnection in pair plasmas with multiplicities $\kappa\sim 1-20$, magnetizations $\sigma\sim 1-10$, and temperatures $\Theta_{e}\sim 1-10$.

Acknowledgments. The authors thank the anonymous referee for their report that helped to improve the manuscript. MP acknowledges support from the Lyman Jr. Spitzer Postdoctoral Fellowship and the Fermi Guest Investigation grant 80NSSC18K1745. LS acknowledges support from DoE DE-SC0016542, NASA ATP NNX-17AG21G, NSF ACI-1657507, and NSF AST-1716567. AS is supported by NASA ATP 80NSSC18K1099 and the Simons Foundation (grant 267233). DG acknowledges support from the NASA ATP NNX17AG21G and the NSF AST-1816136 grants. The simulations were performed using computational resources at the TIGRESS high-performance computer center at Princeton University, the Texas Advanced Computing Center (TACC) at The University of Texas at Austin (TG-AST180038, TG-AST180001), and at the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.

Appendix A Parameter definitions

We summarize the basic physical parameters that are relevant for this study (see Table 3) and provide their definitions below. The total (all-species) plasma magnetization is defined as:

[TABLE]

where $B_{0}$ is the upstream magnetic field strength and $n_{i},n_{e^{\pm}}$ are the number densities of protons and pairs, respectively, in the upstream region. Particles are initialized with temperatures $T_{i}=T_{e}$. The adiabatic indices for pairs and protons are computed iteratively using the equation of state by Synge (1957). We find that $\hat{\gamma}_{e}\approx 4/3$ $\hat{\gamma}_{i}\approx 5/3$, except for $\Theta_{e}=0.1$ where $\hat{\gamma}_{e}\approx 1.5$. The cold plasma magnetization, which neglects the enthalpy terms is defined by:

[TABLE]

A key parameter in the study of the post-reconnection particle energy distributions (see Sect. 5 and Sect. 6) is the hot pair plasma magnetization, which relates to the total $\sigma$ as:

[TABLE]

The cold pair plasma magnetization is identical to $\sigma_{e,h}$ only for non-relativistically hot plasmas ($\Theta_{e}\ll 1$) and is defined as:

[TABLE]

Similar to the pair plasma, one can define the hot proton plasma magnetization:

[TABLE]

which is $\approx\sigma$ for all our cases. The cold proton plasma magnetization is written as:

[TABLE]

and it is the same as $\sigma_{i,h}$ as long as $\Theta_{e}\ll m_{i}/m_{e}$.

The ratio of the electron plasma pressure and the magnetic pressure (plasma $\beta_{e}$), which is a key parameter in studies of electron-proton reconnection, relates to $\sigma_{e,h}$ as:

[TABLE]

If all particle species are relativistically hot ($\Theta_{e}\gg m_{i}/m_{e}$), then $\sigma_{e,h}\approx 2\sigma/(2-q)$ and $\beta_{e}$ reaches its maximum value $\approx 1/4\sigma$.

Let $\omega_{p}$ denote the all-species plasma frequency:

[TABLE]

where the electron, positron, and proton plasma frequencies are given by:

[TABLE]

[TABLE]

and

[TABLE]

Finally, we define the Larmor radius of electrons and protons with Lorentz factors $\sigma_{e,c}$ and $\sigma_{i,c}$, respectively, assuming that all the magnetic energy is transferred to the particles:

[TABLE]

and

[TABLE]

Figures 15 and 16 show the various magnetizations and particle Larmor radii as a function of the proton fraction $q$, which is related to the pair multiplicity $\kappa$ as $q=2/(\kappa+1)$.

Appendix B Appearance of the reconnection layer

In Sect. 3.2 we explored the effects of the pair multiplicity on the appearance of the layer. More specifically, we showed that the layer becomes more structured (i.e., more secondary plasmoids) as the pair multiplicity decreases, for all other parameters kept the same. One could argue that these differences are merely a result of the different box sizes in terms of the proton skin depth. A straightforward way of checking this possibility is to compare cases with different physical conditions, but similar box sizes in terms of $\rho_{Li}$. Snapshots of the density structure from three such pairs of simulations are presented in Figs. 17-19. These comparative plots clearly show that the appearance of the layer is significantly affected by the plasma conditions.

Appendix C Effects of box size

We discuss the effect of the box size on the inflow and outflow rates as well as on the post-reconnection particle energy distributions.

We selected two simulations (see runs A3-A5, C4-C5 in Table 1) and varied the box size in the $x$-direction, as indicated in Fig. 20. Although the peak inflow rate is systematically higher for smaller box sizes, the difference is less than $\sim 3-5\%$. The temporal evolution of the reconnection rate is similar for all box sizes (top panel in Fig. 20), until the formation of the boundary island inhibits the inflow of plasma in the reconnection region, as shown in the bottom panel (blue line). The asymptotic outflow four-velocity is independent of the box size, even for layer lengths of only a few hundred $\rho_{Le}$.

Snapshots of the post-reconnection particle energy distributions from simulations with different box sizes are shown in Fig. 21. The power-law segment of the pair energy spectra is similar for the different cases, suggesting a saturation of the power-law slope already for boxes as small as $L\sim 300\rho_{Le}$ (see also Ball et al., 2018). Thus, we are confident that the power-law slopes we report in Sect. 6.3 (Fig. 12), which were obtained for the spectra plotted with blue lines in Fig. 21, are robust. The high-energy cutoff of the pair distribution, however, increases (almost linearly) with increasing box size, as shown more clearly in the right plot of Fig. 21. Even larger domains are needed for capturing the asymptotic temporal evolution of the cutoff energy. The proton distribution depends strongly on the box size, for both $\sigma$ values we considered. A well-developed power-law forms above the peak proton energy in the largest simulations, thus supporting the argument that reconnection results in extended non-thermal proton distributions (see also Sect. 6.3).

The effects of the box size on the quantities discussed above and in Sect. 5 are summarized in Fig. 22. The outflow four-velocities are not included in this plot, because they are almost the same for the box sizes we considered.

Appendix D Dependence of the mean lepton Lorentz factor on physical parameters

The mean energy of the relativistic pair distribution is of astrophysical importance, as it can be imprinted on the radiated non-thermal photon spectra (for details, see Sect. 7). We therefore attempted to quantify the dependence of mean lepton Lorentz factor on the physical parameters ($\sigma,\Theta_{e}$, and $\kappa$) using a proxy of $\langle\gamma_{e}-1\rangle$, as defined in Sect. 5. We caution the reader that the latter does not necessarily refer to a pure power-law energy distribution. In fact, the definition of $\langle\gamma_{e}-1\rangle$ is agnostic to the shape of the lepton energy distribution.

In general, we find that $\langle\gamma_{e}-1\rangle$ can be described by:

[TABLE]

where $a,b$, and $\chi$ are obtained from a $\chi^{2}$ fit to the data. The best-fit values and the associated $1\sigma$ statistical errors are summarized in Table 4. We note that cases with $\sigma=1,\Theta_{e}=0.1$; $\sigma=10,\Theta_{e}=1$; and $\sigma=10,\Theta_{e}=10$ are excluded from the fit, since the number of $\kappa$ values is the same or less than the free parameters of Eq. (D1). Nevertheless, we still find that $\langle\gamma_{e}-1\rangle\propto\kappa^{-1}$.

Appendix E Effects of plasma temperature on pair energy spectra

The post-reconnection particle energy distributions obtained for the highest temperature simulations (H1-H3 in Table 1) show a high-energy component that forms at late times, as illustrated in Fig. 23 (left panel). This can be described by a power law with slope $p\sim 3.2-3.6$ (see Table 2), which is harder than the power laws obtained for lower temperatures but similar $\sigma_{e,h}$ values (see Fig. 12).

Snapshots of the pair energy distributions from simulations with the same magnetization and multiplicity, but different plasma temperatures, are shown in the middle and right panels of Fig. 23. Although for $\Theta_{e}=100$ there is a prominent high-energy component in the distributions that is independent of $\kappa$, we see a hint of this component at lower temperatures ($\Theta_{e}=10$) only at $\kappa=199$ (right panel). These results imply that the high-energy component of the spectrum is not just related to the plasma temperature. The common denominator in all the cases that show the high-energy component is the high $\beta_{e}$ (i.e., $\beta_{e}>0.1$; see Table 2).

Similar results have been reported by Ball et al. (2018) for trans-relativistic reconnection in electron-proton plasmas with high $\beta_{e}$ approaching the maximum value $1/4\sigma$ (when both electrons and protons start as relativistically hot). The formation of the high-energy component was attributed to a Fermi-like acceleration of particles with initial energy $\sim kT_{e}$ bouncing between the reconnection outflow and the stationary boundary island (see Sect. 6.3 in Ball et al. (2018)). The fact that it takes some time for the boundary island to grow, it is in agreement with the late-time formation of the high-energy component in the spectrum.