TY - JOUR AU1 - Ortuño-Macías,, José AU2 - Nalewajko,, Krzysztof AB - ABSTRACT We present the results of two-dimensional particle-in-cell (PIC) simulations of relativistic magnetic reconnection (RMR) in electron–positron plasma, including the dynamical influence of the synchrotron radiation process, and integrating the observable emission signatures. The simulations are initiated with a single Harris current layer with a central gap that triggers the RMR process. We achieve a steady-state reconnection with unrestricted outflows by means of open boundary conditions. The radiative cooling efficiency is regulated by the choice of initial plasma temperature Θ. We explore different values of Θ and of the background magnetization σ0. Throughout the simulations, plasmoids are generated in the central region of the layer, and they evolve at different rates, achieving a wide range of sizes. The gaps between plasmoids are filled by smooth relativistic outflows called minijets, whose contribution to the observed radiation is very limited due to their low-particle densities. Small-sized plasmoids are rapidly accelerated; however, they have lower contributions to the observed emission, despite stronger relativistic beaming. Large-sized plasmoids are slow but produce most of the observed synchrotron emission, with major part of their radiation produced within the central cores, the density of which is enhanced by radiative cooling. Synchrotron light curves show rapid bright flares that can be identified as originating from mergers between small/fast plasmoids and large/slow targets moving in the same direction. In the high-magnetization case, the accelerated particles form a broken power-law energy distribution with a soft tail produced by particles accelerated in the minijets. acceleration of particles, magnetic reconnection, plasmas, relativistic processes, methods: numerical 1 INTRODUCTION Signatures of non-thermal particle energy distributions are commonly observed in the high-energy astrophysical phenomena such as gamma-ray bursts (e.g. Piran 2004), pulsar wind nebulae (e.g. Buehler et al. 2012), and blazars (e.g. Madejski & Sikora 2016). What these diverse sources have in common is very broad photon spectra indicating broad non-thermal energy distributions of the radiating particles. Blazars are a subclass of active galactic nuclei (AGN), in which a relativistic jet emerging from the supermassive black hole (SMBH) is aligned with our line of sight (Urry & Padovani 1995). They are extremely luminous sources of radiation that spans the entire electromagnetic spectrum from radio up to very high-energy γ-rays (∼TeV). The extreme broadness and power-law appearance of the blazar spectra are signatures of an efficient non-thermal particle acceleration (NTPA) process. Two major spectral components are typically observed with regular characteristics (Fossati et al. 1998; Ghisellini et al. 2017), with the low-energy one (radio to optical/UV/X-rays) interpreted as synchrotron emission, and the high-energy one (mainly γ-rays) due to either leptonic (inverse Compton) or hadronic processes (e.g. Sikora et al. 2009; Böttcher et al. 2013). The emission of blazars is also characterized by strong variability on time-scales ranging from decades (e.g. Ahnen et al. 2016; Goyal et al. 2017) to minutes (e.g. Aharonian et al. 2007,2009; Albert et al. 2007; Meyer, Scargle & Blandford 2019). The shortest variability time-scales are much shorter than the light-crossing time (∼ hours) of the gravitational radius of the SMBH contained in these type of AGN (MBH ∼ 109M⊙). Such short variability time-scales and high luminosities require very compact emitting regions within the jet (Levinson 2007; Begelman, Fabian & Rees 2008; Katarzyński et al. 2008). Compact regions containing relativistic particles and magnetic fields have been often invoked to explain the multiavelength emission of blazars (e.g. Kirk, Rieger & Mastichiadis 1998; Chiaberge & Ghisellini 1999). Theoretical models for the formation of collimated jets with relativistic velocities require relativistically magnetized environments1 (Blandford & Znajek 1977; Li, Chiueh & Begelman 1992; Komissarov et al. 2007). In such environments, shock waves are weak structures that do not provide enough energy dissipation to power the observed emission (Sironi, Petropoulou & Giannios 2015). Magnetic reconnection is a dissipation mechanism that is most efficient in strongly magnetized plasmas. During the reconnection, free magnetic energy is transferred to the thermal heating and bulk acceleration of plasma, as well as to the NTPA process (Zweibel & Yamada 2009). Numerous arguments point to the relativistic magnetic reconnection (RMR) as the process responsible for the multiwavelength and multi-time-scale variability emission of blazars (Sironi et al. 2015). A basic requirement for reconnection to occur in relativistic jets is to have local inversions of the magnetic field lines, which may be caused by internal jet instabilities, in particular, the current-driven kink modes (Begelman 1998; Giannios & Spruit 2006; Alves, Zrake & Fiuza 2018; Bromberg et al. 2019), or by global magnetic polarity reversals (Lovelace, Newman & Romanova 1997; Giannios & Uzdensky 2019). Traditional models of magnetic reconnection (Parker 1957; Sweet 1958; Petschek 1964) have been adapted to the relativistic regime by Lyutikov & Uzdensky (2003) and Lyubarsky (2005). The latter work became a basis for the minijets model of relativistic bulk outflows driven by the RMR (Giannios, Uzdensky & Begelman 2009; Nalewajko et al. 2011). Other works emphasized the fundamental role of plasmoids (or magnetic flux tubes) that arise spontaneously in sufficiently long and thin current layers due to the tearing instability (Loureiro, Schekochihin & Cowley 2007; Giannios 2013). Plasmoids trap energized particles and evacuate them together with the reconnected magnetic field from the active magnetic X-points, enhancing the overall reconnection rate2 (Uzdensky, Loureiro & Schekochihin 2010). The NTPA mechanism in the case of RMR has been studied extensively by means of kinetic particle-in-cell (PIC) numerical simulations. In the most basic case of electron–positron pair plasma, it has been established that NTPA produces power-law energy distributions of particles N(E) ∝ E−p with localized power-law indices as hard as p ≃ 1 in the limit of σ ≫ 1 (Zenitani & Hoshino 2001; Guo et al. 2014; Sironi & Spitkovsky 2014; Werner et al. 2016), converging, however, to p ≃ 2 in the long term (Petropoulou & Sironi 2018). Most PIC simulations of the RMR are initiated from Harris-type current layers contained in periodic domains. In the relatively small domains, this leads to the cancellation of reconnection outflows momenta and to the formation of artificially large single plasmoids. An alternative approach is to use open boundaries that allow outgoing particles to escape freely and absorb the outgoing Poynting flux (Daughton, Scudder & Karimabadi 2006). This approach allows to investigate magnetic reconnection as a sustained steady-state process with unimpeded outflows approaching the terminal Alfv´en velocity and continuous generation and evolution of plasmoids (Daughton & Karimabadi 2007). The open-boundary approach has been first applied to the case of RMR by Sironi, Giannios & Petropoulou (2016), who investigated the statistical properties of plasmoids, determined their size distribution as a power law, and demonstrated an anticorrelation between plasmoid growth and bulk acceleration. The corresponding predictions for the variability of non-thermal radiation were applied to explain the observational characteristics of blazar flares (Petropoulou, Giannios & Sironi 2016). Based on the plasmoid scaling laws, a stochastic model for the evolution of plasmoid chains was developed by Petropoulou et al. (2018). The results of these PIC simulations were also post-processed to calculate light curves of synchrotron and inverse Compton emission as would be observed if the reconnection region was located in a relativistic jet (Christie et al. 2019a, 2020). A semi-analytical model of broad-band emission from reconnection plasmoids in the context of blazar flares was also developed by Morris, Potter & Cotter (2019). Other numerical studies of the RMR were performed in the context of γ-ray flares from the Crab Nebula, taking into account the effect of radiation reaction on individual particles (e.g. Cerutti et al. 2013, 2014). Recent PIC simulations that included the effects of synchrotron and inverse Compton cooling have shown that the particle acceleration can be decreased in comparison with previous non-radiative results (Nalewajko 2018; Hakobyan, Philippov & Spitkovsky 2019; Schoeffler et al. 2019; Sironi & Beloborodov 2019). We may therefore expect that evolution of plasmoids will be affected by radiative cooling. In this work, we present the results of PIC simulations of the steady-state RMR process allowed by the open boundaries under strong radiative energy losses due to the synchrotron radiation reaction. We investigate in detail the evolution of individual plasmoids under different levels of radiative cooling and calculate accurate light curves of synchrotron radiation. We demonstrate a connection between rapid radiation flares and tail-on mergers of small/fast and large/slow plasmoids. Plan of this work: Section 2: simulations setup; Section 3: analysis methods; Section 4.1: initial sequence; Section 4.2: evolved reconnection layer; Section 4.3: spacetime diagrams of the current layer; Section 4.4: particle density and mean energy distributions; Section 4.5: analysis of individual plasmoids; Section 4.6: acceleration and cooling of selected individual particles; Section 4.7: synchrotron light curves; Section 4.8: global energy conservation; Section 4.9: energy distributions of particles and photons; Section 5: discussion; and Section 6: conclusions. 2 SIMULATION SETUP We make use of a custom version of the PIC code Zeltron (Cerutti et al. 2013). We perform two-dimensional simulations of relativistically magnetized pair plasma in a fixed tall domain of dimensions Lx × Ly with Ly = 4Lx, within the coordinate ranges 0 < x < Lx and −Ly/2 < y < Ly/2. Our simulations are initiated from an equilibrium configuration |$\boldsymbol {B}(t=0) = \boldsymbol {B}_{\rm ini}$| and |$\boldsymbol {E}(t=0) = \boldsymbol {E}_{\rm ini} = 0$| that involves a single Harris-type current layer placed in the middle of the computational domain: $$\begin{eqnarray*} B_{\rm ini,x} = -B_0\tanh {(y/\delta)}\, , \end{eqnarray*}$$(1) where B0 is the characteristic value of magnetic field strength, and δ is the Harris layer half-thickness. We include no initial guide field; hence, Bini, z = 0. The magnetic field gradient across the Harris layer is supported by the electric current and pressure provided by the drifting population of particles that is characterized by the Maxwell–Jüttner energy distribution with dimensionless temperature Θd = kBTd/(mec2), Lorentz-boosted with the drift velocity βd = vd/c = 0.3, and number density $$\begin{eqnarray*} n_{\rm d}(y) = n_{\rm d,0}\cosh ^{-2}(y/\delta)\, , \end{eqnarray*}$$(2) where |$n_{\rm d,0} = \gamma _{\rm d} B_0^2 / (8\pi \Theta _{\rm d}m_{\rm e}c^2)$|⁠, and |$\gamma _{\rm d} = (1-\beta _{\rm d}^2)^{-1/2}$| is the drift Lorentz factor. The Harris layer thickness is related to the nominal plasma gyroradius ρ0 = Θdmec2/(eB0) as δ = 2ρ0/(γdβd) (Kirk & Skjæraasen 2003). The current layer is immersed in a background plasma with the Maxwell–Jüttner energy distribution with dimensionless temperature Θb = kBTb/(mec2) = Θd ≡ Θ, and number density nb that is determined by the assumed value of the magnetization |$\sigma _0 = B_0^2/(4\pi n_{\rm b} \Theta _{\rm b} m_{\rm e}c^2)$|⁠. The corresponding initial density contrast can be expressed as nd, 0/nb = γdσ0/2. Our choice of drift velocity means that the current layer is relatively thick and hence stable to the tearing modes. We experimented with higher values βd ≃ 0.5, which results in spontaneous formation of slowly evolving plasmoids. In order to speed up the evolution of the current layer, we follow Sironi et al. (2016) in triggering fast magnetic reconnection by placing a narrow gap (Δx = 8δ) in the middle of the current layer. Within the trigger gap, the distribution of drifting particles is initialized with temperature reduced by a factor of 10, and with no drift velocity. The dimensions of our computational domain are up to Nx × Ny = 4608 × 18432 cells with equal spatial resolution for both dimensions dx = dy = ρ0/3, which results in the physical scales of Lx × Ly = (1536 × 6144)ρ0. The temporal resolution is |${\rm d}t = 0.9\, {\rm d}t_{\rm CFL}$|⁠, where dtCFL = [(dx)−2 + (dy)−2]−1/2/c is the Courant–Friedrichs–Lewy (CFL) time-step. The initial distributions of background and drifting electrons/positrons are represented on average by 64 macroparticles per cell per species. Our computational domain is open at left/right boundaries, where outgoing particles are removed, and fresh particles are injected at every time-step at a fixed rate calculated to maintain the initial number density of both background and drifting particles. At every time-step, we inject a fixed number of particles along each boundary dNinj/d(ct) = (〈β〉/4)dNini/dx, where dNini is the number of initial particles in a single grid column, and |$\left\langle \beta \right\rangle = \left\langle |\boldsymbol \beta |\right\rangle = \left\langle v\right\rangle /c \simeq 1$| is the mean particle dimensionless velocity module.3 The injected particles are located at fixed distance just within each boundary (x0 = Δx/100 in case of left boundary; x0 = Lx − Δx/100 in case of right boundary). Their distribution along the y coordinate is proportional to the distribution of initial particles. The injected particles are likewise divided into background and drifting populations. A Lorentz boost is applied to the drifting particles in order to match the current density |$j_{z,\rm ini}$| of the initial Harris layer (otherwise, the initial current layer would evolve strongly from the boundaries inwards). Within the left/right boundaries, we place the field-absorbing layers of thickness |$\Delta _{\rm abs} = 30\, {\rm d}x$|⁠. In the absorbing layers, in addition to the standard time advance of magnetic and electric fields, at every time-step we perform the following operation: $$\begin{eqnarray*} \boldsymbol {B}(x) &\rightarrow & \boldsymbol {B}(x) + \lambda (x)\left[\boldsymbol {B}_{\rm ini}(x) - \boldsymbol {B}(x)\right]\, , \end{eqnarray*}$$(3) $$\begin{eqnarray*} \boldsymbol {E}(x) &\rightarrow & \boldsymbol {E}(x) + \lambda (x)\left[\boldsymbol {E}_{\rm ini}(x) - \boldsymbol {E}(x)\right]\, , \end{eqnarray*}$$(4) where λ(x) = 0.5(|x − xabs|/Δabs)3, and xabs is the position of the absorbing layer inner edge. For the top/bottom boundaries, we apply periodic conditions for the particles and fields, while only the z-component of the electric field is absorbed using the above rule. We apply the synchrotron radiation reaction to every particle at every time-step,4 following the prescription of Cerutti et al. (2013): $$\begin{eqnarray*} \frac{\partial \boldsymbol {u}}{\partial t} = -\frac{P_{\rm syn} \boldsymbol {u} }{\gamma m_{\rm e}c^2}, \,\, \,\,\,\,\,\, P_{\rm syn} = \frac{\sigma _{\rm T}c}{4\pi }\left[(\gamma \boldsymbol {E} + \boldsymbol {u}\times \boldsymbol {B})^2 - (\boldsymbol {E}\cdot \boldsymbol {u})^2\right]\, , \end{eqnarray*}$$(5) where |$\sigma _{\rm T} = (8\pi /3)r_{\rm e}^2$| is the Thomson cross-section with re = e2/mec2 the classical electron radius, and |$\boldsymbol {u} = \gamma \boldsymbol \beta$| is the dimensionless momentum of a particle with dimensionless velocity |$\boldsymbol \beta = \boldsymbol {v}/c$| and Lorentz factor γ = (1 − β2)−1/2 = (1 + u2)1/2. Noting that for γ ≫ 1 we have u ≃ γ, the radiative cooling rate can be estimated as: $$\begin{eqnarray*} \frac{\partial \gamma }{\partial t} = \frac{\boldsymbol {u}}{\gamma }\cdot \frac{\partial \boldsymbol {u}}{\partial t} \simeq -\frac{P_{\rm syn}}{m_{\rm e} c^2}\, . \end{eqnarray*}$$(6) In the limit of zero electric field and isotropic particle distribution, equation (5) reduces to Psyn = (4/3)σTcu2UB, 0, where |$U_{B,0} = B_0^2/8\pi$| is the initial background magnetic energy density. Taking into account that for the initial Maxwell–Jüttner distribution we have 〈γ〉 ≃ 3Θ and 〈γ2〉 ≃ 12Θ2, we define the nominal cooling length as (Nalewajko, Yuan & Chruślińska 2018): $$\begin{eqnarray*} l_{\rm cool} = c\tau _{\rm cool} = \frac{\langle \gamma \rangle }{\langle \vert {\rm d}\gamma /c{\rm d}t\vert \rangle } \simeq \frac{\langle \gamma \rangle }{\langle \gamma ^2\rangle } \frac{3m_{\rm e}c^2}{4\sigma _{\rm T}U_{\rm B,0}} \simeq \frac{(3\pi /2)e}{\sigma _{\rm T} \Theta ^2 B_0} \rho _0\, . \end{eqnarray*}$$(7) This can be compared with the radiation reaction limit on electron Lorentz factor γrad = (3ρ0/2Θre)1/2 that can be achieved locally under crossed electric and magnetic fields of equal strengths E = B (e.g. Cerutti et al. 2014). We find a relation between lcool and γrad by eliminating B0: $$\begin{eqnarray*} \frac{\left\langle \gamma \right\rangle }{\gamma _{\rm rad}} &=& \left(\frac{27}{8}\frac{\rho _0}{l_{\rm cool}}\right)^{1/2}\, . \end{eqnarray*}$$(8) The total synchrotron power emitted in all directions per volume element of the initial background plasma is given by: $$\begin{eqnarray*} \mathcal {P}_{\rm syn,b,0} \equiv \frac{{\rm d}E_{\rm syn}}{{\rm d}x\, {\rm d}y\, {\rm d}t} = 16c\sigma _{\rm T}n_{\rm b}\Theta ^2U_{\rm B,0}\, . \end{eqnarray*}$$(9) We also collect the synchrotron radiation spectra that would be measured by observers placed at the left and right sides of the simulation domain. The total spectrum is calculated as the sum of contributions from all macroparticles present in the domain (see Nalewajko et al. 2018, and references therein). $$\begin{eqnarray*} L_{\rm syn}(\nu) = \frac{\sqrt{3}e^2}{c} \sum _{e^+ e^-} N_{\rm e,1}F(\xi)\Omega _1\, , \end{eqnarray*}$$(10) where |$F(\xi) = \xi \int _{\xi }^{\infty }K_{5/3}(x){\rm d}x$|⁠, ξ = 4πν/3γ2Ω1, |$\Omega _1 = (e/m_{\rm e}c)\left|(\boldsymbol {E} + \boldsymbol {n} \times \boldsymbol {B}) \times \boldsymbol {n}\right|$|⁠, |$\boldsymbol {n}=\boldsymbol {v}/|\boldsymbol {v}|$|⁠, Ne, 1 is the number of charged particles (either electrons or positrons) represented by a single macroparticle, and K5/3 is the modified Bessel function of the second kind. A characteristic synchrotron frequency for monoenergetic electrons of Lorentz factor γ can be derived from the isotropic synchrotron kernel5: $$\begin{eqnarray*} \nu _{\rm syn}(\gamma) \simeq (3/2\pi)\tilde{\xi }_0\gamma ^2\Theta \omega _0 \equiv \left(\frac{\gamma }{\Theta }\right)^2\nu _{\rm syn0}\, . \end{eqnarray*}$$(11) We performed several large simulations for different values of background magnetization σ0 and initial particle distribution temperature Θ. In Table 1, we report these parameters, as well as the corresponding ratios lcool/Lx and 〈γ〉/γrad. In our strong-cooling case s10Th, lcool is comparable to Lx, while the initial mean particle energy 〈γ〉 is factor ≃ 25 below γrad. Table 1. Designations and key parameter values for our largest simulations performed in the domain of physical width Lx = 1536ρ0. Common parameters include: the initial background magnetic field strength B0 = 1 G, and the initial dimensionless drift velocity βd = 0.3. The nominal cooling length lcool is calculated from equation (7), and the radiation reaction limit ratio 〈γ〉/γrad is calculated from equation (8). Name . σ0 . Θ . lcool/Lx . 〈γ〉/γrad . s10Tl 10 2 · 105 55 0.006 s10Tm 10 5 · 105 8.9 0.016 s10Th 10 1.25 · 106 1.4 0.039 s50Tm 50 5 · 105 8.9 0.016 Name . σ0 . Θ . lcool/Lx . 〈γ〉/γrad . s10Tl 10 2 · 105 55 0.006 s10Tm 10 5 · 105 8.9 0.016 s10Th 10 1.25 · 106 1.4 0.039 s50Tm 50 5 · 105 8.9 0.016 Open in new tab Table 1. Designations and key parameter values for our largest simulations performed in the domain of physical width Lx = 1536ρ0. Common parameters include: the initial background magnetic field strength B0 = 1 G, and the initial dimensionless drift velocity βd = 0.3. The nominal cooling length lcool is calculated from equation (7), and the radiation reaction limit ratio 〈γ〉/γrad is calculated from equation (8). Name . σ0 . Θ . lcool/Lx . 〈γ〉/γrad . s10Tl 10 2 · 105 55 0.006 s10Tm 10 5 · 105 8.9 0.016 s10Th 10 1.25 · 106 1.4 0.039 s50Tm 50 5 · 105 8.9 0.016 Name . σ0 . Θ . lcool/Lx . 〈γ〉/γrad . s10Tl 10 2 · 105 55 0.006 s10Tm 10 5 · 105 8.9 0.016 s10Th 10 1.25 · 106 1.4 0.039 s50Tm 50 5 · 105 8.9 0.016 Open in new tab 3 ANALYSIS METHODS 3.1 Plasmoids identification We describe here an automated procedure that we use for identifying individual plasmoids along the reconnection layer. During the simulations, we frequently save one-dimensional x-profiles of various parameters along the reconnection layer; more precisely, we extract these data from narrow rectangular stripes extending over 0 < x < Lx and |y| < δ/2. For each grid cell along the x-axis, we average the data over all grid cells along the y-axis that fit within the stripe. The x-profiles were then Gaussian-smoothed with the standard deviation of |$3\, {\rm d}x$|⁠. An example of such profiles is shown in Fig. 1. Figure 1. Open in new tabDownload slide Profiles of particle density (blue line), magnetic field component By (green line), and magnetic vector potential component Az (magenta line) measured along the reconnection layer for a selected moment of simulation s10Tm (cf. Fig. 3). The green vertical stripes mark the horizontal limits of plasmoid cores, and the light magenta areas mark the plasmoid layers. Figure 1. Open in new tabDownload slide Profiles of particle density (blue line), magnetic field component By (green line), and magnetic vector potential component Az (magenta line) measured along the reconnection layer for a selected moment of simulation s10Tm (cf. Fig. 3). The green vertical stripes mark the horizontal limits of plasmoid cores, and the light magenta areas mark the plasmoid layers. We applied the scipy.signal.find_peaks algorithm to identify the maxima of particle number density n (with the minimum prominence parameter of 1.5nb), the minima and maxima of magnetic field component By (with the minimum prominence of 0.02B0), and the minima and maxima of the magnetic vector potential component Az (with the minimum prominence of 0.22B0ρ0). Next, we searched for all ordered pairs consisting of a By minimum followed (to the right) by a By maximum. For each such pair, we checked whether between those two points there is exactly one maximum of n and one maximum of Az. We also require that the distance between the n maximum and the Az maximum is no more than |$4\, {\rm d}x$|⁠. Only the structures that satisfied all of the above conditions were classified as proper plasmoids. We distinguish two structures within each plasmoid, a core delimited by the minimum and maximum peaks of By (green stripes in Fig. 1), and a layer where Az exceeds a threshold level equal to the higher of the closest Az minima located on both sides of the core and excluding the core itself (light-magenta areas in Fig. 1; this is consistent with the definition of plasmoid limits adopted by Sironi et al. 2016). Qualitatively, the plasmoid cores are compact structures of high particle density, while the plasmoid layers are extended structures of relatively low particle density and closed magnetic field lines. 3.2 Local reference frames We convert certain parameters to their values in the local reference frames. As discussed in Werner et al. (2018), in relativistically hot fluid there are two alternative ways to define a ‘co-moving’ reference frame: (1) zero-current (Eckart) frame |$\mathcal {O}^{\prime }$| of velocity |$\boldsymbol \beta _1 = \left\langle \boldsymbol \beta \right\rangle$|⁠; and (2) zero-momentum (Landau) frame |$\mathcal {O}^{\prime \prime }$| of velocity |$\boldsymbol \beta _2 = \left\langle \gamma \boldsymbol \beta \right\rangle / \left\langle \gamma \right\rangle$|⁠. These averages are calculated locally for all particles (electrons and positrons) in a given grid cell as 〈a〉 = (∑iaini)/(∑ini), where ni are the multiplicities of individual macroparticles such that ∑ini ≡ n is the local particle number density. For example, we calculate the electric field component in the Eckart frame as: $$\begin{eqnarray*} E_z^{\prime } = \Gamma _1 \left(E_z + \left\langle \beta _x\right\rangle B_y - \left\langle \beta _y\right\rangle B_x\right) \end{eqnarray*}$$(12) where Γ1 = (1 − 〈βx〉2 − 〈βy〉2)−1/2 is the Eckart bulk Lorentz factor. On the other hand, we calculate the mean particle energy in the Landau frame as 〈γ″〉 = 〈γ〉/Γ2, where |$\Gamma _2 = (1 - \beta _{2,x}^2 - \beta _{2,y}^2)^{-1/2}$| is the Landau bulk Lorentz factor. 4 RESULTS 4.1 Initial sequence Fig. 2 shows the initial evolution of the reconnection layer in our reference simulation s10Tm. Each panel shows a particle number density map and magnetic field lines for a central region of our computational domain (only a 1/64 fraction of its total vertical extent), in which the current layer is contained. A centrally located gap in the current layer expands rapidly sideways towards the left/right boundaries and the initial drifting plasma component is pulled towards the boundaries by the tension of closed magnetic field lines. Magnetic reconnection is triggered in the central low-density region, sustained by a thin current layer that is unstable to tearing modes, which lead to the emergence of plasmoids. Figure 2. Open in new tabDownload slide Initial sequence of simulation s10Tm, presenting the logarithm of particle number density n/nb (see the top panel of Fig. 3 for the colour scale) and the magnetic field lines (solid white lines). The dashed magenta lines indicate the limits of the field-absorbing boundary layers. Figure 2. Open in new tabDownload slide Initial sequence of simulation s10Tm, presenting the logarithm of particle number density n/nb (see the top panel of Fig. 3 for the colour scale) and the magnetic field lines (solid white lines). The dashed magenta lines indicate the limits of the field-absorbing boundary layers. At ct ≃ Lx, the swept-up fronts of the initial drifting plasma leave the left/right boundaries. From that moment on, we can describe the simulated reconnection process as steady-state, with newer plasmoids continuously generated around the centre of the layer, and older plasmoids escaping through the left/right boundaries. The upper panels of Fig. 2 demonstrate how our implementation of open boundaries works when large dense structures move across them. We also find that the initial trigger gap induces a quasi-circular wave that propagates in all directions. The horizontal wave fronts are effectively absorbed by the left/right boundaries, while the vertical wave fronts reach the top/bottom boundaries by ct ≃ 2Lx and are partially absorbed there; some very weak reflections return to the reconnection layer by about ct ≃ 4Lx without a significant effect. 4.2 Evolved reconnection layer Fig. 3 presents a selected evolved state of our reference simulation s10Tm illustrated with multiple parameters. In this snapshot centred at the reconnection midplane, we observe diverse substructures. The main current layer is observed at 0.34 < x/Lx < 0.49, characterized by high values of the specific electric current |jz|/jmax (electric current density jz close to its maximum value jmax = cen) despite low particle number density n, low bulk velocity vx, and strong electric field |$E_z^{\prime }$| (except for the two minor plasmoids located at x ≃ 0.39Lx and x ≃ 0.425Lx). Figure 3. Open in new tabDownload slide Selected evolved state of simulation s10Tm (the same as shown in Fig. 1). From the top, the panels show: (1) particle number density n/nb in log scale, (2) out-of-plane electric current per particle number density jz/(cen), (3) Landau velocity component along the current layer vx/c, (4) Landau velocity component across the current layer vy/c, (5) mean particle Lorentz factor measured in the simulation frame 〈γ〉/Θ, (6) mean particle Lorentz factor measured in the Landau co-moving frame 〈γ″〉/Θ, (7) out-of-plane electric field component measured in the Eckart co-moving frame |$E_z^{\prime }$|⁠, and (8) total synchrotron power Esyn in log scale. The solid white lines show the magnetic field lines, and the dashed magenta lines indicate the limits of the field-absorbing boundary layers. Figure 3. Open in new tabDownload slide Selected evolved state of simulation s10Tm (the same as shown in Fig. 1). From the top, the panels show: (1) particle number density n/nb in log scale, (2) out-of-plane electric current per particle number density jz/(cen), (3) Landau velocity component along the current layer vx/c, (4) Landau velocity component across the current layer vy/c, (5) mean particle Lorentz factor measured in the simulation frame 〈γ〉/Θ, (6) mean particle Lorentz factor measured in the Landau co-moving frame 〈γ″〉/Θ, (7) out-of-plane electric field component measured in the Eckart co-moving frame |$E_z^{\prime }$|⁠, and (8) total synchrotron power Esyn in log scale. The solid white lines show the magnetic field lines, and the dashed magenta lines indicate the limits of the field-absorbing boundary layers. To the left of the main current layer, behind a medium-sized plasmoid at x ≃ 0.3Lx, we find a relativistically fast reconnection outflow (strongly negative vx for x < 0.25Lx; 〈γ〉 ≫ 〈γ″〉 due to Lorentz transformation), also characterized by moderate particle density and intrinsic temperature 〈γ″〉/Θ, specific electric current |jz|/jmax decreasing systematically with distance, very weak intrinsic electric field |$E_z^{\prime }$|⁠, and weak synchrotron emission. This is a structure that has all characteristics of minijets – regular reconnection outflows of relativistic bulk velocity (Giannios et al. 2009; Nalewajko et al. 2011). We can easily recognize the conical geometry of the outflow region with parallel outflow velocity field (very low values of vy), and oblique magnetic field lines crossing the outflow boundaries, as has been described by an analytical model of relativistic Petschek-type reconnection by Lyubarsky (2005). There is one qualitative difference from that model – the magnetic field lines in the outflow region are not vertical and the magnetic field gradient ∂Bx/∂y does not vanish in that region and it is supported by the non-zero electric current density jz. We note that there is a roughly uniform vertical inflow of background plasma into the minijet region with velocity vy (reconnection rate) of the same order as that of the inflow into the main current layer. To the right of the main current layer, we find a group of several plasmoids of various sizes, all propagating to the right at different velocities vx > 0. The largest plasmoid can be seen centred at x ≃ 0.725Lx; it is clearly slower than its smaller neighbour centred at x ≃ 0.6Lx. The smaller plasmoid is in the process of merging with the large one, even though they both propagate in the same direction. This is a natural consequence of the inverse relation between the growth and bulk acceleration of plasmoids that has been first noticed by Sironi et al. (2016). 4.3 Spacetime diagrams Most of the information on evolution of current layers, and especially on the plasmoids, is contained along the reconnection midplane; this information can be presented very efficiently in the form of spacetime diagrams. Following the practice of Nalewajko et al. (2015), the one-dimensional parameter x-profiles described in Section 3.1 are combined into spacetime diagrams of high time resolution. Fig. 4 shows spacetime diagrams for several parameters for our reference simulation s10Tm. Figure 4. Open in new tabDownload slide Spacetime diagrams extracted along the reconnection midplane for the simulation s10Tm. From left to right we plot: (1) the magnetic field component By/B0, (2) the electric field component Ez/B0, (3) the Landau-type bulk velocity component βx, (4) logarithm of the cold magnetization parameter log (σ/σ0), (5) mean particle Lorentz factor 〈γ″〉/Θ measured in the Landau frame, and (6) logarithm of the particle number density log (n/nb). In the first five panels, we show the particle number density contours n = 7nb with green solid lines. We indicate the simulation time ct = 3.36Lx corresponding to the simulation state shown in Figs 1 and 3 with the magenta–dashed horizontal lines. The black–dotted horizontal lines in the last panel indicate the initial simulation states shown in Fig. 2. Figure 4. Open in new tabDownload slide Spacetime diagrams extracted along the reconnection midplane for the simulation s10Tm. From left to right we plot: (1) the magnetic field component By/B0, (2) the electric field component Ez/B0, (3) the Landau-type bulk velocity component βx, (4) logarithm of the cold magnetization parameter log (σ/σ0), (5) mean particle Lorentz factor 〈γ″〉/Θ measured in the Landau frame, and (6) logarithm of the particle number density log (n/nb). In the first five panels, we show the particle number density contours n = 7nb with green solid lines. We indicate the simulation time ct = 3.36Lx corresponding to the simulation state shown in Figs 1 and 3 with the magenta–dashed horizontal lines. The black–dotted horizontal lines in the last panel indicate the initial simulation states shown in Fig. 2. The spacetime diagrams reveal a sustained bifurcating outflow along the x coordinate (regions of negative and positive velocity component βx) and a variety of plasmoids (indicated by enhanced plasma density and sharp positive gradients of magnetic field component dBy/dx > 0), most of which are generated along the βx ≃ 0 line. There are a few large plasmoids that evolve very slowly, their bulk velocities are non-relativistic, acceleration time-scale is long, and hence they spend more than the light-crossing time Lx/c before leaving the simulation domain. On the other hand, there are many small plasmoids that accelerate quickly to relativistic bulk velocities, and they spend less than Lx/c in the simulation domain. Small plasmoids may either escape through the boundaries or merge with a larger plasmoid. As the velocity field is generally divergent, plasmoid mergers are relatively rare; typically, they involve plasmoids of different sizes moving in the same direction (tail-on mergers). A large plasmoid can attract nearby small plasmoids and can even reverse their motion, e.g. Fig. 3 shows that at ct/Lx = 3.36, a small plasmoid centred at x ≃ 0.85Lx has βx > 0; however, the spacetime diagrams reveal that it will merge with the large plasmoid located on its left side by ct/Lx ≃ 3.6. The regions between plasmoids are characterized by roughly uniform electric field component Ez ≃ 0.1B0 and by two-value magnetic field component By ≃ ±0.15B0 (positive where βx < 0 and negative where βx > 0). The uniformity of electric field measured in the simulation frame is consistent with the uniform reconnection rate indicated in the βy panel of Fig. 3. It is remarkable, given the non-uniformity of electric field measured in local Eckart frames |$E_z^{\prime }$|⁠, as shown in another panel of Fig. 3. This indicates a smooth connection between the minijet outflows (where |$E_z^{\prime }$| is close to zero) and the proper magnetic diffusion areas (where |$E_z^{\prime }$| is strong). High-density regions of the plasmoids – essentially the plasmoid cores – are initially characterized by enhanced temperatures (mean particle energies measured in the local Landau frames reaching values of 〈γ″〉 ≃ 10Θ) due to the heating of particles in the reconnection process. Our spacetime diagram of 〈γ″〉/Θ demonstrates that the cores of large plasmoids cool down gradually. In Fig. 5, we compare the spacetime diagrams of 〈γ″〉/Θ for three simulations with different nominal cooling lengths lcool. In the simulation s10Tl, in which lcool > Lx, we find that plasmoid temperatures are the highest, showing no signs of cooling over the light-crossing time-scale Lx/c. On the other hand, in the simulation s10Th, the plasmoid temperatures are the lowest and they eventually become even lower than for the plasma between the plasmoids. Figure 5. Open in new tabDownload slide Spacetime diagrams of average particle Lorentz factor 〈γ″〉/Θ measured in the Landau frame compared for three simulations characterized by different radiative cooling efficiencies. Figure 5. Open in new tabDownload slide Spacetime diagrams of average particle Lorentz factor 〈γ″〉/Θ measured in the Landau frame compared for three simulations characterized by different radiative cooling efficiencies. 4.4 Plasma density and temperature distributions The left-hand panel of Fig. 6 compares the distributions of particle number density n/nb based on the spacetime diagrams of three simulations with σ0 = 10 and with different efficiencies of radiative cooling. The distributions form power-law tails with the slope of ∼−2.3, independent of the cooling efficiency. However, the highest values of the particle density increase with stronger cooling from nmax ≃ 150nb for the simulation s10Tl to nmax ≃ 1500nb for the simulation s10Th. The highest particle densities are realized in the cores of large plasmoids. Strong radiative cooling removes the gas pressure support, allowing the plasmoid cores to contract further. Figure 6. Open in new tabDownload slide Distributions of particle density n/nb (left-hand panel) and mean particle energy 〈γ″〉/Θ calculated in the Landau co-moving frame (right-hand panel) extracted from the spacetime diagrams probing narrow regions along the reconnection layers for the three simulations with σ0 = 10 that differ by initial gas temperature Θ and hence by the radiative cooling efficiency (higher Θ corresponds to more efficient cooling). Figure 6. Open in new tabDownload slide Distributions of particle density n/nb (left-hand panel) and mean particle energy 〈γ″〉/Θ calculated in the Landau co-moving frame (right-hand panel) extracted from the spacetime diagrams probing narrow regions along the reconnection layers for the three simulations with σ0 = 10 that differ by initial gas temperature Θ and hence by the radiative cooling efficiency (higher Θ corresponds to more efficient cooling). The right-hand panel of Fig. 6 compares the distributions of mean particle energy 〈γ″〉/Θ measured in the local Landau frames, also extracted from the spacetime diagrams of simulations with σ0 = 10. This shows that the distributions are strongly peaked at values that decrease slightly with increasing cooling efficiency from 〈γ″〉peak ≃ 5.4Θ for the simulation s10Tl to 〈γ″〉peak ≃ 4.8Θ for the simulation s10Th. 4.5 Individual plasmoids Figs 7–9 present the histories of individual plasmoids compared for the four simulations listed in Table 1. The plasmoid histories are compiled from the parameters determined independently for each time-step in automated plasmoid identification algorithm described in Section 3.1. As such, they are subject to some irregularities, in particular, due to shifting positions of the minima of magnetic potential Az (that define the plasmoid boundaries), especially for large plasmoids merging with the small ones. In order to reduce these irregularities, we chose not to record plasmoid histories in the following circumstances: (1) when there is no gap between two separate plasmoids (i.e. they share the same local Az minimum) and (2) when the plasmoid layer boundary is less than 100 dx from the left/right domain boundary. Figure 7. Open in new tabDownload slide The histories of individual plasmoids compared for the four large simulations. From the top, the rows present: (1) plasmoid core width wc/ρ0; (2) total plasmoid width w/ρ0; (3) mean core density nc/nb; and (4) peak magnetic field strength |By|/B0. For each plasmoid, time is measured relative to its first appearance. The line colour indicates the peak total plasmoid width. Figure 7. Open in new tabDownload slide The histories of individual plasmoids compared for the four large simulations. From the top, the rows present: (1) plasmoid core width wc/ρ0; (2) total plasmoid width w/ρ0; (3) mean core density nc/nb; and (4) peak magnetic field strength |By|/B0. For each plasmoid, time is measured relative to its first appearance. The line colour indicates the peak total plasmoid width. Figure 8. Open in new tabDownload slide Further histories of individual plasmoids. From the top, the rows present: (1) mean core velocity |βc, x|; (2) mean particle energy 〈γ〉c/Θ measured in the simulation frame, averaged over the plasmoid core; (3) mean Landau-frame particle energy 〈γ″〉c/Θ of the plasmoid core; and (4) mean Landau-frame particle energy 〈γ″〉l/Θ of the plasmoid layer. See Fig. 7 for more description. Figure 8. Open in new tabDownload slide Further histories of individual plasmoids. From the top, the rows present: (1) mean core velocity |βc, x|; (2) mean particle energy 〈γ〉c/Θ measured in the simulation frame, averaged over the plasmoid core; (3) mean Landau-frame particle energy 〈γ″〉c/Θ of the plasmoid core; and (4) mean Landau-frame particle energy 〈γ″〉l/Θ of the plasmoid layer. See Fig. 7 for more description. Figure 9. Open in new tabDownload slide Further histories of individual plasmoids. From the top, the rows present: (1) mean synchrotron emissivity |$\left\langle \mathcal {P}_{\rm syn,c}\right\rangle / \mathcal {P}_{\rm syn,b,0}$| of the plasmoid core and (2) mean synchrotron emissivity |$\left\langle \mathcal {P}_{\rm syn,l}\right\rangle / \mathcal {P}_{\rm syn,b,0}$| of the plasmoid layer. See Fig. 7 for more description, and equation (9) for the definition of |$\mathcal {P}_{\rm syn,b,0}$|⁠. Figure 9. Open in new tabDownload slide Further histories of individual plasmoids. From the top, the rows present: (1) mean synchrotron emissivity |$\left\langle \mathcal {P}_{\rm syn,c}\right\rangle / \mathcal {P}_{\rm syn,b,0}$| of the plasmoid core and (2) mean synchrotron emissivity |$\left\langle \mathcal {P}_{\rm syn,l}\right\rangle / \mathcal {P}_{\rm syn,b,0}$| of the plasmoid layer. See Fig. 7 for more description, and equation (9) for the definition of |$\mathcal {P}_{\rm syn,b,0}$|⁠. We present separately the total plasmoid widths w (dominated by the widths of plasmoid layers) and the widths of their cores wc, noting that there is no correlation between them. In all simulations, the plasmoid core widths, after a brief formation phase, show a decreasing trend. On the other hand, the total plasmoid widths grow systematically in time for every simulation. We find that the densities of plasmoid cores nc grow systematically in time for all simulations. The larger the plasmoid, the denser its core becomes. For simulations with σ0 = 10, higher core densities are reached for higher radiative efficiencies. In the case of σ0 = 50, the core densities are even higher, roughly in proportion to σ0. The histories of plasmoid core velocities βc, x confirm the picture discussed before of small plasmoids being accelerated very rapidly to relativistic velocities and of large plasmoids being accelerated slowly only to mildly relativistic velocities. The mean particle energies of plasmoid cores, measured in the simulation frame, are somewhat higher for small plasmoids. However, when measured in the Landau frame of the core, they are instead higher for large plasmoids until the radiative cooling effects become significant. The difference is mainly due to relativistic bulk motions of the small plasmoids. In particular, in the simulations s10Tl and s10Tm characterized by weak/moderate cooling, small plasmoids have 〈γ″〉c ≃ 6Θ, roughly constant in time, while large plasmoids can reach 〈γ″〉c ≃ 11Θ in the case s10Tl. In the simulation s10Th characterized by strong cooling, the intrinsic mean particle energies of the core are reduced down to 〈γ″〉c ≃ 2Θ by Lx/c. In the simulation s50Tm, the cores of large plasmoids are heated up to 〈γ″〉c ≃ 23Θ before cooling down to 〈γ″〉c ≲ 10Θ. In contrast, the plasmoid layers are characterized by similar and stable intrinsic temperatures 〈γ″〉c ≃ (4 − 6)Θ for all simulations with σ0 = 10 (and 〈γ″〉c ∼ 15Θ for σ0 = 50), showing virtually no signs of radiative cooling. Most of the plasmoids show their synchrotron emissivity building up in time, both for the core and for the layer. Some of the smaller plasmoids show a very slow decline. The cores of large plasmoids produce significantly more (2–3 orders of magnitude) synchrotron emission per volume element than the cores of small plasmoids and also about 2 orders of magnitude higher emission density than the large plasmoid layers. Our basic finding is that synchrotron emission of the plasmoids is sustained over a long term, irrespective of the cooling efficiency. For the cores of large plasmoids that undergo the most efficient radiative cooling in our high-temperature simulations, it appears that systematically increasing core density, as well as systematically increasing magnetic field strength, is able to offset the reduction in particle mean energy. In the case of σ0 = 50, we find a systematic increase of the emission from the cores of large plasmoids to the levels up to |$10^7 \mathcal {P}_{\rm syn,b,0}$|⁠. 4.6 Individual particles Here we describe the behaviour of four selected energetic positrons (denoted as particles #1–#4) for which a detailed history has been recorded. The spacetime tracks of these particles are indicated in Fig. 10. From this, one can see that particles #1 and #2 become trapped in different plasmoids, while particles #3 and #4 become accelerated in the low-density regions between plasmoids. Figure 10. Open in new tabDownload slide Spacetime diagram of the tracks of selected energetic particles, the acceleration of which is characterized in detail in Fig. 11. The line colour indicates the instantaneous particle energy measured in the simulation frame. Particle density contours n = 7nb are indicated with grey lines. Figure 10. Open in new tabDownload slide Spacetime diagram of the tracks of selected energetic particles, the acceleration of which is characterized in detail in Fig. 11. The line colour indicates the instantaneous particle energy measured in the simulation frame. Particle density contours n = 7nb are indicated with grey lines. Fig. 11 presents the detailed acceleration histories of these four particles, all in the same time window. Particle #1 becomes accelerated to Lorentz factor γ ≃ 22Θ in a single episode. The beginning of acceleration episode coincides with the particle becoming trapped in the reconnection midplane (|y| < 5ρ0). We can also see that the acceleration episode coincides with a formation of a new plasmoid that traps the particle also in the x coordinate. The particle experiences mainly the electric field component Ez ∼ 0.1B0, and its momentum gain is at first mainly in the z direction, later also in the x direction. After the acceleration episode, the particle oscillates around the plasmoid centre, which results in oscillations of its Lorentz factor γ measured in the simulation frame. However, its Lorentz factor γ′ measured in the plasmoid frame shows a slow gradual decline in time, which we attribute to the radiative cooling. Figure 11. Open in new tabDownload slide Acceleration histories for selected energetic particles, the x-positions of which are shown in Fig. 10. For each particle, we present detailed information as functions of simulation time on five panels, from the left: (1) the y-position measured from the reconnection midplane, (2) particle Lorentz factor γ (also indicated with a colour scale on panels 1 and 2) normalized to Θ, (3) three components of the particle 4-velocity |$\boldsymbol {u}$| (x – red, y – green, and z – blue) normalized to Θ, (4) three components of the local electric field |$\boldsymbol {E}$| in units of B0, and (5) three components of the local magnetic field |$\boldsymbol {B}$| in units of B0. For particles #1 and #2, the black lines indicate parameter values (Lorentz factor γ′, electric field |$E_\parallel ^{\prime }$|⁠, and magnetic field |$B_\perp ^{\prime }$|⁠) measured in the co-moving frame of a plasmoid to which the particle is attached. Figure 11. Open in new tabDownload slide Acceleration histories for selected energetic particles, the x-positions of which are shown in Fig. 10. For each particle, we present detailed information as functions of simulation time on five panels, from the left: (1) the y-position measured from the reconnection midplane, (2) particle Lorentz factor γ (also indicated with a colour scale on panels 1 and 2) normalized to Θ, (3) three components of the particle 4-velocity |$\boldsymbol {u}$| (x – red, y – green, and z – blue) normalized to Θ, (4) three components of the local electric field |$\boldsymbol {E}$| in units of B0, and (5) three components of the local magnetic field |$\boldsymbol {B}$| in units of B0. For particles #1 and #2, the black lines indicate parameter values (Lorentz factor γ′, electric field |$E_\parallel ^{\prime }$|⁠, and magnetic field |$B_\perp ^{\prime }$|⁠) measured in the co-moving frame of a plasmoid to which the particle is attached. Particle #2 shows a similar behaviour to particle #1; the acceleration episode also coincides with the formation of a new plasmoid. In this case, the acceleration episode is shorter and the energy gain is also lower, with acceleration proceeding in similar electric fields. The radiative cooling is less efficient, and it should be noted that the co-moving perpendicular magnetic field is about twice weaker. At ct ≃ 2.1Lx, the plasmoid in which particle #2 is trapped merges into a larger plasmoid on its right side. We see that this results in particle #2 losing most of its energy measured in the simulation frame. It appears that this particle did not experience direct deceleration by strong electric field, instead it just happened to be at its lower energy level in the oscillation cycle at the moment of plasmoids merger. It has also been kicked out from the reconnection midplane, settling briefly at y ∼ 15ρ0, trailing behind the merged plasmoid before it exits the right boundary of the domain. Particle #3 was accelerated in a long acceleration episode under constant electric field Ez ∼ 0.12B0 after becoming trapped in the y coordinate to the reconnection midplane. The y and Bx data reveal a typical Speiser orbit with gradually increasing period. The acceleration episode is interrupted at ct ≃ 2.1Lx, when the particle begins to interact with a large plasmoid, which forces the particle away from the reconnection midplane. The particle bounces twice off the trailing edge of the plasmoid before it exits the left boundary of the domain. Particle #4 shows another example of Speiser-orbit acceleration in the low-density region of reconnection midplane that is enabled by trapping the particle in the y coordinate. 4.7 Synchrotron emissivity and light curves Fig. 12 presents the spacetime distribution of the total synchrotron power emitted from the reconnection midplane in all directions for simulation s10Tm, and light curves that would be received by two distant observers placed at the opposite sides of the simulation domain, one on the left (−x-axis) and one on the right (+x-axis). While the spacetime diagram of synchrotron power is based on the x-profiles integrated over narrow stripes of |y| < δ/2, the light curves are calculated from a much wider region of |y| < Lx/2 to ensure contribution from the whole plasmoids and their broader surroundings. Figure 12. Open in new tabDownload slide The left-hand and right-hand panels show synchrotron light curves measured by observers placed at left and right sides of the reconnecting layer, respectively. The cyan and orange lines correspond to the two selected frequency bands indicated in Fig. 15. The middle panel shows the spacetime distribution of total synchrotron power normalized to the value |$\mathcal {P}_{\rm syn,b,0}$| defined in equation (9). The particle number density contours n = 7nb are shown with green solid lines. The black dashed lines represent the light-cones corresponding to selected features in either light curve. The magenta dashed line indicates the simulation time ct = 3.36Lx corresponding to the simulation state shown in Figs 1 and 3. Figure 12. Open in new tabDownload slide The left-hand and right-hand panels show synchrotron light curves measured by observers placed at left and right sides of the reconnecting layer, respectively. The cyan and orange lines correspond to the two selected frequency bands indicated in Fig. 15. The middle panel shows the spacetime distribution of total synchrotron power normalized to the value |$\mathcal {P}_{\rm syn,b,0}$| defined in equation (9). The particle number density contours n = 7nb are shown with green solid lines. The black dashed lines represent the light-cones corresponding to selected features in either light curve. The magenta dashed line indicates the simulation time ct = 3.36Lx corresponding to the simulation state shown in Figs 1 and 3. The synchrotron emission is strongly concentrated along the plasmoid trajectories. As we have noted in Section 4.5, large plasmoids produce significantly more synchrotron emission. The emissivity of the large plasmoid cores exceeds the emissivity of the background plasma by four orders of magnitude. Light curves received by either observer are completely different, indicating a high level of anisotropy. The most conspicuous features of the light curves are very sharp and bright flares. We can identify the origin of these flares by drawing the corresponding light-cones on the spacetime diagram to the events of small plasmoids approaching the observer with relativistic velocity and merging with a large target plasmoid that is also approaching the observer. Zooming up on these flares, we find their characteristic time-scales of τ ∼ ρ0/c, so they are basically unresolved.6 Because of their very short duration, the contribution of these flares to the overall radiation fluence is rather insignificant. The light curves also contain smooth structures characterized by relatively long rise and short decline. These can be attributed to large plasmoids propagating with mildly relativistic velocities, and the sharp declines observed in the light curves coincide with these plasmoids exiting the simulation domain. The contribution of these plasmoids to the light curves is not complete, since their emission is not contained within the simulation boundaries, as we have already noted in Section 4.5. Comparing the light curves recorded in two frequency bands,7 we find that in general the contribution from large plasmoids is higher in the lower frequency band (cyan lines), while the flares produced by small plasmoids are stronger in the higher frequency band (orange lines). Light curves obtained from different simulations are qualitatively very similar, in particular, the level of radiative cooling efficiency does not clearly affect the light-curve characteristics. 4.8 Energy conservation Because of the use of open boundaries with steady injection of fresh particles, our simulations do not conserve energy globally. In order to evaluate the efficiency of energy transformations, we defined a fixed region |$\mathcal {R}$| centred around the reconnection midplane between the left/right absorbing boundary layers, defined by 2Δabs < x < (Lx − 2Δabs) and −Lx/4 < y < Lx/4. In addition to the instantaneous energy contained in |$\mathcal {R}$| in the form of magnetic and electric fields, as well as in the particles, we also calculate the cumulative energy emitted by all particles in the synchrotron process and the fluxes of particles and electromagnetic fields (i.e. the Poynting flux) inflowing/outflowing across the |$\mathcal {R}$| boundaries. In Fig. 13, we present the time evolution of different forms of energy contained in the region |$\mathcal {R}$| for the simulation s10Tm. At the beginning of the simulation, the region |$\mathcal {R}$| is dominated by magnetic energy (⁠|$\mathcal {E}_{\rm B,0} \simeq 0.6\, \mathcal {E}_{\rm tot}$|⁠). The initial (ct/Lx < 0.6) energization of particles at the cost of magnetic energy is due to the trigger mechanism. This is followed by the somewhat erratic variation of the particle energy, which reflects systematic heating by magnetic reconnection and episodic escapes of large plasmoids. Over the course of the simulation (ct/Lx ≃ 4.5), the magnetic energy of the region |$\mathcal {R}$| decreases by |$\simeq 40{{\ \rm per\ cent}}$|⁠, while the particle energy decreases only by about |$\simeq\! 15{{\ \rm per\ cent}}$|⁠. At the same time, we measure a large net influx of electromagnetic energy (accumulating to |$\simeq 1.5\, \mathcal {E}_{\rm B,0}$|⁠), mainly through the top/bottom boundaries of the region |$\mathcal {R}$|⁠, and even larger net outflow of particle energy, mainly through the left/right boundaries. The net energy outflow (particle minus electromagnetic) through the region boundaries amounts to |$\simeq 0.25\, \mathcal {E}_{B,0}$| of the initial magnetic energy, which is slightly less than the particle energy lost to the synchrotron radiation (⁠|$\simeq 0.3\, \mathcal {E}_{B,0}$|⁠). Accounting for all these energy components and flows, the total energy in |$\mathcal {R}$| is conserved at the |$\sim 0.1{{\ \rm per\ cent}}$| level. Figure 13. Open in new tabDownload slide Top panel: mean energy densities, normalized to initial magnetic energy density UB, 0, calculated for the analysis region |$\mathcal {R}$|⁠, defined by 2Δabs < x < Lx − 2Δabs and −Lx/4 < y < Lx/4, for the simulation s10Tm. We plot the contributions from magnetic fields (green), electric fields (magenta), all particles (e+ and e−; red), and total conserved energy including inflows and outflows (black). The black dotted line shows the total instantaneous energy contained in |$\mathcal {R}$|⁠, including contributions from the magnetic and electric fields, and all particles. Middle panel: energy flux inflows and outflows: the inflowing Poynting flux (green), outflowing particle energy flux (red), and total synchrotron emission (blue). Bottom panel: conservation of the total energy of region |$\mathcal {R}$|⁠, including inflows and outflows. Figure 13. Open in new tabDownload slide Top panel: mean energy densities, normalized to initial magnetic energy density UB, 0, calculated for the analysis region |$\mathcal {R}$|⁠, defined by 2Δabs < x < Lx − 2Δabs and −Lx/4 < y < Lx/4, for the simulation s10Tm. We plot the contributions from magnetic fields (green), electric fields (magenta), all particles (e+ and e−; red), and total conserved energy including inflows and outflows (black). The black dotted line shows the total instantaneous energy contained in |$\mathcal {R}$|⁠, including contributions from the magnetic and electric fields, and all particles. Middle panel: energy flux inflows and outflows: the inflowing Poynting flux (green), outflowing particle energy flux (red), and total synchrotron emission (blue). Bottom panel: conservation of the total energy of region |$\mathcal {R}$|⁠, including inflows and outflows. 4.9 Energy distributions of particles and photons Fig. 14 shows the energy distributions of all particles: electrons and positrons. For each simulation, it is averaged over a period of time that excludes only the initial stage (ct/L ≲ 0.85). In all studied cases, the particle energy distributions established after the initial period show no significant evolution in time. As energetic particles escape across the open left/right boundaries, other particles are energized across the current layer and are subject to radiative energy losses within the plasmoids. The balance between these processes is maintained regardless of the efficiency of radiative cooling. In all studied cases, a small fraction of particles reach energies of γcutoff = 4σ0Θ established as a cutoff energy in a previous study of non-radiative Harris-layer reconnection within periodic boundaries (Werner et al. 2016). In the case of σ0 = 10, we find only a minor effect of radiative cooling in limiting the high-energy excess for Θ ≳ 106. In the case of σ0 = 50, the high-energy component can be described as a broken power law with a hard slope of p1 ≃ 1.5 extending up to γ ≃ 25Θ and a soft tail of p2 ≃ 3.6 extending up to γ ≃ 150Θ. In that case, we also have γcutoff ≃ γrad. Figure 14. Open in new tabDownload slide Energy distributions of all particles contained in the simulation domain, compared for the four main simulations (higher Θ corresponds to more efficient cooling). For each simulation, the distribution is averaged over simulation time, excluding the initial stage (ct/L ≲ 0.85), in the space of flux logarithm. The thin dashed lines indicate the corresponding standard deviation values. The thick grey dashed line represents the initial Maxwell–Jüttner distribution for Θ ≫ 1. The distributions are presented in arbitrary units and they are normalized to match the low-energy sections. The vertical dotted lines indicate the characteristic values of γ/Θ: 4 (the γ2N(γ) peak for the initial background particles; grey), 4 × 10 (red),and 4 × 50 (black). The oblique black dotted lines indicate two power-law slopes p (N(γ) ∝ γ−p) along the σ0 = 50 distribution. The four stars indicate the values of γrad/Θ for each simulation. Figure 14. Open in new tabDownload slide Energy distributions of all particles contained in the simulation domain, compared for the four main simulations (higher Θ corresponds to more efficient cooling). For each simulation, the distribution is averaged over simulation time, excluding the initial stage (ct/L ≲ 0.85), in the space of flux logarithm. The thin dashed lines indicate the corresponding standard deviation values. The thick grey dashed line represents the initial Maxwell–Jüttner distribution for Θ ≫ 1. The distributions are presented in arbitrary units and they are normalized to match the low-energy sections. The vertical dotted lines indicate the characteristic values of γ/Θ: 4 (the γ2N(γ) peak for the initial background particles; grey), 4 × 10 (red),and 4 × 50 (black). The oblique black dotted lines indicate two power-law slopes p (N(γ) ∝ γ−p) along the σ0 = 50 distribution. The four stars indicate the values of γrad/Θ for each simulation. Fig. 15 shows the spectral energy distributions (SED) νFν of the synchrotron emission produced by all particles in all directions, averaged over the same periods of time as the particle energy distributions presented in Fig. 14. In the case of σ0 = 10, the SED are dominated by the contribution from low-energy particles peaking around ν ≃ 19νsyn0, with a high-frequency excess extending beyond a characteristic value of |$\nu _{\rm cutoff} \simeq 19\sigma _0^2\nu _{\rm syn0}$|⁠. The level of this high-frequency excess increases with decreasing gas temperature Θ, which means that radiative cooling suppresses the high-frequency radiation component more clearly than it affects the high-energy particle tail. In the case of σ0 = 50, the SED is strongly dominated by the contribution from energetic particles with the maximum photon energies consistent with a cutoff at |$\nu _{\rm cutoff} \simeq 19\sigma _0^2\nu _{\rm syn0}$|⁠, which coincides with the radiation reaction limit νsyn, max = (γrad/Θ)2νsyn0. We note that the SED shape around its peak is not described by a broken power law corresponding directly to that indicated in the electron distribution (with slopes νFν ∝ ν−s; s = (p − 3)/2; and characteristic frequencies νi/νsyn0 ≃ (γi/Θ)2). This is because the extent of the electron energy distribution that can be described as a broken power law is too short to result in a broken power-law photon spectrum when folded with the synchrotron kernel. Figure 15. Open in new tabDownload slide Isotropic spectra of the synchrotron radiation emitted across the simulation domain, compared for the four main simulations. For each simulation, the distribution is averaged over simulation time, excluding the initial stage (ct/L ≲ 0.85), in the space of flux logarithm. The thin dashed lines indicate the corresponding standard deviation values. The thick grey dashed line represents the synchrotron spectrum of the initial Maxwell–Jüttner distribution. The frequencies are normalized to the characteristic synchrotron frequency νsyn0 defined in equation (11). The distributions are presented in arbitrary units and they are normalized to match the low-frequency sections. The vertical dotted lines indicate the characteristic values of ν/νsyn0: 19 (the νF(ν) peak for the initial background particles; grey), 19 × 102 (red), and 19 × 502 (black). The oblique black dotted lines indicate two power-law slopes s = (3 − p)/2 (νFν ∝ ν−s) that would be expected for the corresponding power laws p1, p2 marked in Fig. 14. The four stars indicate the values of MHD synchrotron frequency limit νsyn, max/νsyn0 = (γrad/Θ)2 for each simulation. The cyan and orange stripes indicate the frequency bands from which the light curves shown in Fig. 12 were extracted. Figure 15. Open in new tabDownload slide Isotropic spectra of the synchrotron radiation emitted across the simulation domain, compared for the four main simulations. For each simulation, the distribution is averaged over simulation time, excluding the initial stage (ct/L ≲ 0.85), in the space of flux logarithm. The thin dashed lines indicate the corresponding standard deviation values. The thick grey dashed line represents the synchrotron spectrum of the initial Maxwell–Jüttner distribution. The frequencies are normalized to the characteristic synchrotron frequency νsyn0 defined in equation (11). The distributions are presented in arbitrary units and they are normalized to match the low-frequency sections. The vertical dotted lines indicate the characteristic values of ν/νsyn0: 19 (the νF(ν) peak for the initial background particles; grey), 19 × 102 (red), and 19 × 502 (black). The oblique black dotted lines indicate two power-law slopes s = (3 − p)/2 (νFν ∝ ν−s) that would be expected for the corresponding power laws p1, p2 marked in Fig. 14. The four stars indicate the values of MHD synchrotron frequency limit νsyn, max/νsyn0 = (γrad/Θ)2 for each simulation. The cyan and orange stripes indicate the frequency bands from which the light curves shown in Fig. 12 were extracted. In order to clarify the connection between the electron energy distribution and synchrotron SED in the case of σ0 = 50, we analysed a sample of individually tracked particles. We selected particles over two energy ranges: (1) a medium-energy range 11 < γ/Θ < 22, corresponding to the hard power-law section of index p1 ≃ 1.5; and (2) a high-energy range 50 < γ/Θ < 150, corresponding to the soft power-law section of index p2 ≃ 3.6. In Fig. 16, we show the distributions of these particles along coordinate x and over local magnetic field strength B, as well as their contributions to the synchrotron SED, taking into account accurate local electromagnetic fields felt by each particle. These distributions are averaged over multiple simulation time-steps for ct/Lx > 0.85. For the medium-energy particles, we find that their distribution along x is fairly uniform, and their distribution over B is broad, with some particles found in strongly amplified magnetic field B > 10B0 characteristic for the plasmoid cores. For the high-energy particles, we find that they are clearly concentrated towards the left/right boundaries, and that they are found almost exclusively in magnetic fields of moderate strength B < 5B0. In addition, we observe that for individual simulation time-steps, the medium-energy particles are clearly concentrated within the plasmoids, while the high-energy particles are diffused over x. The medium-energy particles dominate the synchrotron SED for most frequencies, except the highest values of ν > 5 × 104νsyn0. At low frequencies, the relative contribution of the medium-energy particles is ≃ 60 times higher than that of the high-energy ones. Figure 16. Open in new tabDownload slide Analysis of individual tracked particles for the simulation s50Tm. Particles are selected over two energy ranges – 11 < γ/Θ < 22 (red) and 50 < γ/Θ < 150 (blue) – corresponding to the two power-law sections indicated in Fig. 14. The top panel compares their normalized distributions along coordinate x; the middle panel compares their normalized distributions over magnetic field strength B; and the bottom panel compares their contributions to the isotropic synchrotron SED (arbitrary units). The distributions are averaged over multiple simulation time-steps for ct/Lx > 0.85. The vertical dotted lines in the bottom panel correspond to those in Fig. 15. Figure 16. Open in new tabDownload slide Analysis of individual tracked particles for the simulation s50Tm. Particles are selected over two energy ranges – 11 < γ/Θ < 22 (red) and 50 < γ/Θ < 150 (blue) – corresponding to the two power-law sections indicated in Fig. 14. The top panel compares their normalized distributions along coordinate x; the middle panel compares their normalized distributions over magnetic field strength B; and the bottom panel compares their contributions to the isotropic synchrotron SED (arbitrary units). The distributions are averaged over multiple simulation time-steps for ct/Lx > 0.85. The vertical dotted lines in the bottom panel correspond to those in Fig. 15. 5 DISCUSSION Our results are consistent with the basic picture of steady-state relativistic plasmoid reconnection that has been established since the work of Sironi et al. (2016). The open boundaries allow to evacuate the reconnected plasma and develop unimpeded reconnection outflows that reach relativistic bulk velocities of the order of the background Alfv´en velocity. In the centrally positioned active magnetic X-point, plasmoids are generated spontaneously over a wide range of sizes. The smallest discernible plasmoids, with the width of order ∼20ρ0 (deep blue lines in Figs 7–9), are rapidly accelerated to relativistic speeds and their lifetimes are about ∼0.2Lx/c. The largest plasmoids, with the width of order ∼(150−200)ρ0 (brown lines in Figs 7–9) accelerate slowly, reaching mildly relativistic speeds ∼c/2 only after ∼1.5Lx/c. This is a confirmation of the anticorrelation between growth and acceleration of plasmoids identified by Sironi et al. (2016). An important consequence of the relation between plasmoid sizes and their acceleration length scale is that it is not possible to contain the bulk acceleration of large plasmoids within a domain with open boundaries. We have found that as we increase the size of the simulation domain, given proportionally more simulation time, we obtain ever larger plasmoids that require ever more acceleration time. In other words, we are unable to isolate the bulk acceleration of plasmoids from the domain boundaries and obtain a coasting phase of uniformly Alfv´enic reconnection outflows. A similar problem applies to the production of synchrotron light curves. In our simulations, even for the highest plasma temperatures investigated, we were not able to contain synchrotron emission from the largest plasmoids within the domain boundaries. As we show in Figs 8 and 9, even though in the case of short cooling length lcool ≪ Lx, the plasmoid cores indeed undergo efficient radiative cooling within the Lx/c time-scale, the total synchrotron emission of plasmoid cores and layers is not decreasing significantly. The nominal cooling time-scale cannot be reduced indefinitely by a further increase of the particle distribution temperature. We found that, with unrestricted radiation reaction, already for Θ = 1.25 × 106, the particle densities in the cores of the largest plasmoids increase to the level at which the Debye length eventually becomes unresolved λD = (Θmec2/4πne2)1/2 < dx, which leads to the development of numerical electrostatic instabilities. Our ad hoc restriction of radiation reaction [see a footnote preceding equation (5)] helps to avoid developing these instabilities. We note that a similar restriction of radiation physics (pair creation) in the plasmoid cores has been applied by Hakobyan et al. (2019). Christie et al. (2019a) calculated light curves of synchrotron and inverse Compton emission by post-processing the results of non-radiative PIC simulations of relativistic steady-state reconnection of Sironi et al. (2016). Their radiation transfer model is based on several assumptions that can be verified by the results of our radiative PIC simulations. One of their key assumptions is that the intrinsic structure of plasmoids is not important and can be approximated by using their average parameters that in addition are constant in time. Our results suggest that the synchrotron emissivity is strongly concentrated in the central parts of the plasmoids (see the bottom panel of Fig. 3), and that in the radiatively efficient regime, the plasmoid cores undergo significant time evolution with systematic increase of plasmoid core density and peak magnetic field strength (Fig. 7). We suggest that small plasmoids and the cores of large plasmoids are important for understanding the production of rapid radiation flares. Investigation of these structures is also the most challenging from the numerical perspective. Our study suggests that properly resolving the cores of large plasmoids will be critical for understanding the radiative signatures of plasmoid reconnection. Recent non-radiative PIC simulations of relativistic reconnection demonstrated an important role of large plasmoids in extending the high-energy tail of the particle energy distribution along a power law of slope ≃ 2 (Petropoulou & Sironi 2018). However, taking into account radiative cooling, which is expected to be particularly strong in the plasmoid cores, the maximum energy achievable in the plasmoids may be significantly limited. Plasmoid mergers have been suggested previously to be important for particle acceleration in relativistic reconnection, based on the results of relatively modest PIC simulations of reconnection in periodic boundaries (Nalewajko et al. 2015). Subsequent numerical studies emphasized a more fundamental role of magnetic X-points as a crucial first step for particles that eventually achieve the highest energies (Ball, Sironi & Özel 2019; Guo et al. 2019). Here, we would like to point out that even if plasmoid mergers may not dominate particle acceleration in non-radiative reconnection, they are important for the production of transient radiation signals. Even if major head-on collisions of large plasmoids are rare events, tail-on collisions of unequal plasmoids, arguably more frequent events due to the growth–acceleration anticorrelation, can be responsible for most of the sharpest features observed in the resulting light curves. Recent results of Christie et al. (2020) show that large plasmoids can attract many small plasmoids originating on both sides of their trajectory, enhancing the rates of both tail-on and head-on mergers. Our results also demonstrate the co-existence of plasmoids and minijets in the same reconnection layer. We find minijets persisting in the gaps forming between plasmoids, plasmoids are able to slide along a minijet without causing much disturbance, and a minijet reforms behind a passing plasmoid. The structure of minijets found in our simulations is qualitatively very similar to the analytical model of Lyubarsky (2005). Although the minijets contain some highly energetic particles, their contribution to the observed radiative signatures appears to be very weak. There may be two reasons behind this: (1) the minijets are characterized by much lower particle density than the plasmoids and (2) energetic particles propagating along the minijets show only weak radiative energy losses due to weak perpendicular magnetic field component.8 In the high-magnetization case of σ0 = 50, we found that the particle energy distribution is maintained in the form of a broken power law with a hard power-law slope p1 ≃ 1.5 breaking around γ ≃ 25Θ into a soft power-law slope p2 ≃ 3.6. The hard power-law slope is consistent with the results of PIC simulations of non-radiative highly relativistic Harris-layer reconnection in periodic boundaries (Guo et al. 2014; Sironi & Spitkovsky 2014; Werner et al. 2016). The soft power-law slope is reminiscent of that seen in the recent PIC simulations of relativistic reconnection with strong inverse-Compton cooling, both in periodic boundaries (Werner, Philippov & Uzdensky 2019) and in the open ones (Sironi & Beloborodov 2019). In the latter work, the electron energy distribution has been decomposed into contributions from particles accelerated at different sites, which suggests that the soft power law arises from a combined action of primary X-points of the main reconnection layer, secondary X-points formed by merging plasmoids, and unstructured outflows (the minijets). This is consistent with our analysis of the distribution of individual tracked particles belonging to the energy range of 50 < γ/Θ < 150. We found that these particles are strongly concentrated towards the boundaries and are found exclusively outside the plasmoid cores at normal magnetic field strengths B < 5B0. Typical examples of such particles are Particles #3 and #4 presented in Figs 10 and 11. Low relative numbers and diffuse spatial distribution imply that these particles are not important to the overall synchrotron SED (except for the very highest frequencies), nor to the production of rapid radiation flares. 6 CONCLUSIONS We presented the results of the first kinetic simulations of RMR within open boundaries that enable steady-state plasmoid reconnection and including the synchrotron radiation reaction. We confirm the general picture of steady-state relativistic plasmoid reconnection established by Sironi et al. (2016) and subsequent works. We find that synchrotron emission of plasmoids cannot be contained within open boundaries. The cores of large plasmoids are the main sites of synchrotron emission and their particle densities are significantly enhanced due to radiative pressure losses. Rapid flares of synchrotron radiation can be produced by tail-on mergers between small/fast plasmoids with large/slow targets. The plasmoids are also found to co-exist with the minijets that do not produce a lot of radiation due to their low particle densities. In the high-magnetization case (σ0 = 50), the energy distribution of accelerated particles can be described as a broken power law with a hard medium-energy section produced mainly by particles accelerated in the plasmoids and a soft high-energy tail produced by diffuse particles accelerated in the minijets. ACKNOWLEDGEMENTS We thank the reviewer for helpful comments. We acknowledge discussions with Benoît Cerutti, Dimitrios Giannios, and Maria Petropoulou. The original version of the Zeltron code was created by Benoît Cerutti and co-developed by Gregory Werner at the University of Colorado Boulder (http://benoit.cerutti.free.fr/Zeltron/). These results are based on numerical simulations performed at the supercomputer Prometheus located at Cyfronet AGH, Poland (PLGrid grant recjose18; PI: K. Nalewajko). This work was supported by the Polish National Science Centre grant 2015/18/E/ST9/00580. DATA AVAILABILITY The data underlying this article will be shared on reasonable request to the corresponding author. Footnotes 1 In general, magnetization parameter is defined as σ = B2/(4πw), where B is the magnetic field strength and w is the relativistic enthalpy density of the matter. 2 Inflow velocity of the background plasma. 3 The factor 〈β〉/4 is derived from the rate at which particles of uniform density and isotropic momentum distribution |$N(\boldsymbol {u}) = N f(u) g(\mu)$| cross the boundary. Here, f(u) is the Maxwell–Jüttner distribution (normalized to unity) of momentum module |$u = |\boldsymbol {u}| = \gamma \beta$| with γ = (1 − β2)−1/2 = (1 + u2)1/2 the Lorentz factor, and g(μ) is the uniform distribution of the cosine parameter μ = ux/u = βx/β. For isotropic target particle distribution, g(μ) = 1/2 for μ ∈ [ − 1: 1]; hence, |$\int _{-1}^{1}{\rm d}\mu \, g(\mu) = 1$|⁠. For injected particles, this distribution is truncated to particles inflowing to the domain – μ ∈ (0: 1] in case of left boundary, and μ ∈ [ − 1: 0) in case of right boundary. For particles of momentum |$\boldsymbol {u}$|⁠, we find |${\rm d}N(\boldsymbol {u})/{\rm d}(ct) = \beta _x\, {\rm d}N(\boldsymbol {u})/{\rm d}x$|⁠, integrating this over |$\boldsymbol {u}$|⁠, we find |${\rm d}N/{\rm d}(ct) = {\rm d}N/{\rm d}x\, \int _{0}^{\infty }{\rm d}u\, \beta f(u)\, \int _{0}^{1}{\rm d}\mu \, \mu g(\mu) = {\rm d}N/{\rm d}x\, \left\langle \beta \right\rangle \, (1/4)$|⁠. 4 We apply a restriction such that particles cannot lose more than half of their energy or reverse their momentum over a single time-step. 5 Defined as |$\mathcal {R}(\tilde{\xi }) = \tilde{\xi }^2\lbrace K_{1/3}(\tilde{\xi })K_{4/3}(\tilde{\xi }) - (3/5)\tilde{\xi }[K_{4/3}^2(\tilde{\xi }) - K_{1/3}^2(\tilde{\xi })]\rbrace$|⁠, where |$\tilde{\xi }= 2\pi \nu /3\gamma ^2\Theta \omega _0$| (Crusius & Schlickeiser 1986). Note that |$\tilde{\xi }\, \mathcal {R}(\tilde{\xi }) \propto \nu \, F(\nu)$| peaks at |$\tilde{\xi }_0 \simeq 0.575$|⁠. 6 Contributions to every light curve are recorded at every simulation time-step at the temporal resolution of dt. 7 The light curves are presented in the same arbitrary units equivalent to the νFν flux density. 8 Radiative cooling in the minijets might be stronger if the guide field component Bz were included. REFERENCES Aharonian F. et al. , 2007 , ApJ , 664 , L71 10.1086/520635 Crossref Search ADS Crossref Aharonian F. et al. , 2009 , ApJ , 696 , L150 10.1088/0004-637X/696/2/L150 Crossref Search ADS Crossref Ahnen M. L. et al. , 2016 , A&A , 593 , A91 10.1051/0004-6361/201628447 Crossref Search ADS Crossref Albert J. et al. , 2007 , ApJ , 669 , 862 10.1086/521382 Crossref Search ADS Crossref Alves E. P. , Zrake J., Fiuza F., 2018 , Phys. Rev. Lett. , 121 , 245101 Crossref Search ADS PubMed Ball D. , Sironi L., Özel F., 2019 , ApJ , 884 , 57 10.3847/1538-4357/ab3f2e Crossref Search ADS Crossref Begelman M. C. , 1998 , ApJ , 493 , 291 10.1086/305119 Crossref Search ADS Crossref Begelman M. C. , Fabian A. C., Rees M. J., 2008 , MNRAS , 384 , L19 10.1111/j.1745-3933.2007.00413.x Crossref Search ADS Crossref Blandford R. D. , Znajek R. L., 1977 , MNRAS , 179 , 433 10.1093/mnras/179.3.433 Crossref Search ADS Crossref Böttcher M. , Reimer A., Sweeney K., Prakash A., 2013 , ApJ , 768 , 54 10.1088/0004-637X/768/1/54 Crossref Search ADS Crossref Bromberg O. , Singh C. B., Davelaar J., Philippov A. A., 2019 , ApJ , 884 , 39 10.3847/1538-4357/ab3fa5 Crossref Search ADS Crossref Buehler R. et al. , 2012 , ApJ , 749 , 26 10.1088/0004-637X/749/1/26 Crossref Search ADS Crossref Cerutti B. , Werner G. R., Uzdensky D. A., Begelman M. C., 2013 , ApJ , 770 , 147 10.1088/0004-637X/770/2/147 Crossref Search ADS Crossref Cerutti B. , Werner G. R., Uzdensky D. A., Begelman M. C., 2014 , ApJ , 782 , 104 10.1088/0004-637X/782/2/104 Crossref Search ADS Crossref Chiaberge M. , Ghisellini G., 1999 , MNRAS , 306 , 551 10.1046/j.1365-8711.1999.02538.x Crossref Search ADS Crossref Christie I. M. , Petropoulou M., Sironi L., Giannios D., 2019 , MNRAS , 482 , 65 10.1093/mnras/sty2636 Crossref Search ADS Crossref Christie I. M. , Petropoulou M., Sironi L., Giannios D., 2020 , MNRAS , 492 , 549 10.1093/mnras/stz3265 Crossref Search ADS Crossref Crusius A. , Schlickeiser R., 1986 , A&A , 164 , L16 Daughton W. , Karimabadi H., 2007 , Phys. Plasmas , 14 , 072303 10.1063/1.2749494 Crossref Search ADS Crossref Daughton W. , Scudder J., Karimabadi H., 2006 , Phys. Plasmas , 13 , 072101 10.1063/1.2218817 Crossref Search ADS Crossref Fossati G. , Maraschi L., Celotti A., Comastri A., Ghisellini G., 1998 , MNRAS , 299 , 433 10.1046/j.1365-8711.1998.01828.x Crossref Search ADS Crossref Ghisellini G. , Righi C., Costamante L., Tavecchio F., 2017 , MNRAS , 469 , 255 10.1093/mnras/stx806 Crossref Search ADS Crossref Giannios D. , 2013 , MNRAS , 431 , 355 10.1093/mnras/stt167 Crossref Search ADS Crossref Giannios D. , Spruit H. C., 2006 , A&A , 450 , 887 10.1051/0004-6361:20054107 Crossref Search ADS Crossref Giannios D. , Uzdensky D. A., 2019 , MNRAS , 484 , 1378 10.1093/mnras/stz082 Crossref Giannios D. , Uzdensky D. A., Begelman M. C., 2009 , MNRAS , 395 , L29 10.1111/j.1745-3933.2009.00635.x Crossref Search ADS Crossref Goyal A. et al. , 2017 , ApJ , 837 , 127 10.3847/1538-4357/aa6000 Crossref Search ADS Crossref Guo F. , Li H., Daughton W., Liu Y.-H., 2014 , PhRvL , 113 , 155005 10.1103/PhysRevLett.113.155005 Crossref Guo F. , Li X., Daughton W., Kilian P., Li H., Liu Y.-H., Yan W., Ma D., 2019 , ApJ , 879 , L23 10.3847/2041-8213/ab2a15 Crossref Search ADS Crossref Hakobyan H. , Philippov A., Spitkovsky A., 2019 , ApJ , 877 , 53 10.3847/1538-4357/ab191b Crossref Search ADS Crossref Katarzyński K. , Lenain J.-P., Zech A., Boisson C., Sol H., 2008 , MNRAS , 390 , 371 10.1111/j.1365-2966.2008.13753.x Crossref Search ADS Crossref Kirk J. G. , Skjæraasen O., 2003 , ApJ , 591 , 366 10.1086/375215 Crossref Search ADS Crossref Kirk J. G. , Rieger F. M., Mastichiadis A., 1998 , A&A , 333 , 452 Komissarov S. S. , Barkov M. V., Vlahakis N., Königl A., 2007 , MNRAS , 380 , 51 10.1111/j.1365-2966.2007.12050.x Crossref Search ADS Crossref Levinson A. , 2007 , ApJ , 671 , L29 10.1086/524708 Crossref Search ADS Crossref Li Z.-Y. , Chiueh T., Begelman M. C., 1992 , ApJ , 394 , 459 10.1086/171597 Crossref Search ADS Crossref Loureiro N. F. , Schekochihin A. A., Cowley S. C., 2007 , Phys. Plasmas , 14 , 100703 10.1063/1.2783986 Crossref Search ADS Crossref Lovelace R. V. E. , Newman W. I., Romanova M. M., 1997 , ApJ , 484 , 628 10.1086/304358 Crossref Search ADS Crossref Lyubarsky Y. E. , 2005 , MNRAS , 358 , 113 10.1111/j.1365-2966.2005.08767.x Crossref Search ADS Crossref Lyutikov M. , Uzdensky D., 2003 , ApJ , 589 , 893 10.1086/374808 Crossref Search ADS Crossref Madejski G. , Sikora M., 2016 , ARA&A , 54 , 725 10.1146/annurev-astro-081913-040044 Crossref Search ADS Crossref Meyer M. , Scargle J. D., Blandford R. D., 2019 , ApJ , 877 , 39 10.3847/1538-4357/ab1651 Crossref Search ADS Crossref Morris P. J. , Potter W. J., Cotter G., 2019 , MNRAS , 486 , 1548 10.1093/mnras/stz920 Crossref Search ADS Crossref Nalewajko K. , 2018 , MNRAS , 481 , 4342 10.1093/mnras/sty2549 Crossref Search ADS Crossref Nalewajko K. , Giannios D., Begelman M. C., Uzdensky D. A., Sikora M., 2011 , MNRAS , 413 , 333 10.1111/j.1365-2966.2010.18140.x Crossref Search ADS Crossref Nalewajko K. , Uzdensky D. A., Cerutti B., Werner G. R., Begelman M. C., 2015 , ApJ , 815 , 101 10.1088/0004-637X/815/2/101 Crossref Search ADS Crossref Nalewajko K. , Yuan Y., Chruślińska M., 2018 , JPlPh , 84 , 755840301 Parker E. N. , 1957 , J. Geophys. Res. , 62 , 509 Crossref Search ADS Petropoulou M. , Sironi L., 2018 , MNRAS , 481 , 5687 10.1093/mnras/sty2702 Crossref Search ADS Crossref Petropoulou M. , Giannios D., Sironi L., 2016 , MNRAS , 462 , 3325 10.1093/mnras/stw1832 Crossref Search ADS Crossref Petropoulou M. , Christie I. M., Sironi L., Giannios D., 2018 , MNRAS , 475 , 3797 10.1093/mnras/sty033 Crossref Search ADS Crossref Petschek H. E. , 1964 , NASA Spec. Publ. , 50 , 425 Piran T. , 2004 , Rev. Mod. Phys. , 76 , 1143 10.1103/RevModPhys.76.1143 Crossref Search ADS Crossref Schoeffler K. M. , Grismayer T., Uzdensky D., Fonseca R. A., Silva L. O., 2019 , ApJ , 870 , 49 10.3847/1538-4357/aaf1b9 Crossref Search ADS Crossref Sikora M. , Stawarz Ł., Moderski R., Nalewajko K., Madejski G. M., 2009 , ApJ , 704 , 38 10.1088/0004-637X/704/1/38 Crossref Search ADS Crossref Sironi L. , Beloborodov A. M., 2019 , preprint (arXiv:1908.08138) Sironi L. , Spitkovsky A., 2014 , ApJ , 783 , L21 10.1088/2041-8205/783/1/L21 Crossref Search ADS Crossref Sironi L. , Petropoulou M., Giannios D., 2015 , MNRAS , 450 , 183 10.1093/mnras/stv641 Crossref Search ADS Crossref Sironi L. , Giannios D., Petropoulou M., 2016 , MNRAS , 462 , 48 10.1093/mnras/stw1620 Crossref Search ADS Crossref Sweet P. A. , 1958 , Electromagn. Phenom. Cosmical Phys. , p. 123 Urry C. M. , Padovani P., 1995 , PASP , 107 , 803 10.1086/133630 Crossref Search ADS Crossref Uzdensky D. A. , Loureiro N. F., Schekochihin A. A., 2010 , Phys. Rev. Lett. , 105 , 235002 Crossref Search ADS PubMed Werner G. R. , Uzdensky D. A., Cerutti B., Nalewajko K., Begelman M. C., 2016 , ApJ , 816 , L8 10.3847/2041-8205/816/1/L8 Crossref Search ADS Crossref Werner G. R. , Uzdensky D. A., Begelman M. C., Cerutti B., Nalewajko K., 2018 , MNRAS , 473 , 4840 10.1093/mnras/stx2530 Crossref Search ADS Crossref Werner G. R. , Philippov A. A., Uzdensky D. A., 2019 , MNRAS , 482 , L60 10.1093/mnrasl/sly157 Crossref Search ADS Crossref Zenitani S. , Hoshino M., 2001 , ApJ , 562 , L63 10.1086/337972 Crossref Search ADS Crossref Zweibel E. G. , Yamada M., 2009 , ARA&A , 47 , 291 10.1146/annurev-astro-082708-101726 Crossref Search ADS Crossref © 2020 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) TI - Radiative kinetic simulations of steady-state relativistic plasmoid magnetic reconnection JF - Monthly Notices of the Royal Astronomical Society DO - 10.1093/mnras/staa1899 DA - 2020-09-11 UR - https://www.deepdyve.com/lp/oxford-university-press/radiative-kinetic-simulations-of-steady-state-relativistic-plasmoid-3KHmeeB0wX SP - 1365 EP - 1381 VL - 497 IS - 2 DP - DeepDyve ER -