Kinetic Simulations of Radiative Magnetic Reconnection in the Coronae of Accreting Black Holes

and

Published 2020 August 12 © 2020. The American Astronomical Society. All rights reserved.
, , Citation Lorenzo Sironi and Andrei M. Beloborodov 2020 ApJ 899 52 DOI 10.3847/1538-4357/aba622

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

0004-637X/899/1/52

Abstract

We perform 2D and 3D particle-in-cell simulations of reconnection in magnetically dominated e± plasmas subject to strong Compton cooling. Magnetic reconnection under such conditions can operate in accretion disk coronae around black holes, which produce hard X-rays through Comptonization. Our simulations show that most of the plasma in the reconnection layer is kept cold by Compton losses and locked in magnetically dominated plasmoids with a small thermal pressure. Compton drag clears cavities inside plasmoids and also affects their bulk motions. These effects, however, weakly change the reconnection rate and the plasmoid size distribution from those in nonradiative reconnection. This demonstrates that the reconnection dynamics is governed by similar magnetic stresses in both cases and weakly affected by thermal pressure. We examine the energy distribution of particles energized by radiative reconnection and observe two distinct components: (1) A mildly relativistic peak, which results from bulk motions of cooled plasmoids. This component receives most of the dissipated reconnection power and dominates the output X-ray emission. The peak has a quasi-Maxwellian shape with an effective temperature of ∼100 keV. Thus, it mimics thermal Comptonization used previously to fit hard-state spectra of accreting black holes. (2) A high-energy tail, which receives ∼20% of the dissipated reconnection power. It is populated by particles accelerated impulsively at X-points or "picked up" by fast outflows from X-points. The high-energy particles immediately cool, and their inverse Compton emission explains the MeV spectral tail detected in the hard state of Cyg X-1. Our first-principle simulations support magnetic reconnection as a mechanism powering hard X-ray emission from magnetically dominated regions of accreting black holes.

Export citation and abstract BibTeX RIS

1. Introduction

Luminous accretion disks around stellar-mass black holes are routinely observed in "soft" and "hard" X-ray states (e.g., Zdziarski & Gierliński 2004). The soft state is dominated by quasi-thermal emission from an optically thick accretion disk, while the hard state is usually explained by a phenomenological model of Comptonization of soft X-rays in a hot plasma ("corona") with electron temperature of ∼100 keV. The electron energization mechanism required to balance the fast inverse Compton (IC) losses of the corona is unknown. Recently, Beloborodov (2017, hereafter B17) proposed that this mechanism is the bulk acceleration of multiscale plasmoids generated by magnetic reconnection.

The corona above the accretion disk has a low mass density ρ and is likely dominated by the disk magnetic field B, with the magnetization parameter σ ≡ B2/4πρc2 > 1. Then, the Alfvén speed in the corona is close to the speed of light and magnetic reconnection proceeds in the relativistic regime. This regime was explored analytically (e.g., Lyutikov & Uzdensky 2003; Lyubarsky 2005) and studied with fully kinetic particle-in-cell (PIC) simulations. The simulations provided a detailed view of relativistic reconnection, including reliable measurements of the reconnection rate, and revealed the mechanisms of particle acceleration (see Kagan et al. 2015 for a review). With unprecedentedly large-scale simulations, Sironi et al. (2016) have investigated the properties of the chain of plasmoids (magnetic islands) that are copiously generated in the reconnection layer, as a self-consistent by-product of the system evolution (Uzdensky et al. 2010).

For conditions appropriate to coronae of accreting black holes, B17 argued that the plasmoid chain experiences fast IC losses accompanied by electron–positron pair creation, estimated the resulting distribution of plasmoid bulk speeds, and performed Monte Carlo simulations of photon Comptonization by the plasmoid motions. The self-regulated plasmoid chain was found to radiate its energy in hard X-rays with a spectrum similar to the hard-state spectra observed in black hole X-ray binaries and active galactic nuclei.

In the present work, we employ kinetic PIC simulations to study relativistic magnetic reconnection with strong IC losses. Numerous previous studies of relativistic electron–positron reconnection with the PIC technique focused on the regime with negligible cooling (e.g., Zenitani & Hoshino 2001; Lyubarsky & Liverts 2008; Kagan et al. 2013; Guo et al. 2014, 2015, 2019; Sironi & Spitkovsky 2014; Sironi et al. 2015, 2016; Nalewajko et al. 2015; Werner et al. 2016; Werner & Uzdensky 2017; Petropoulou & Sironi 2018). One exception is the models of reconnection in pulsars and pulsar wind nebulae, where synchrotron emission from extremely energetic electrons was included in the PIC simulations (e.g., Cerutti et al. 2013, 2014; Kagan et al. 2016, 2018; Hakobyan et al. 2019). Recently, Werner et al. (2019) incorporated IC effects in their PIC simulations (see also Nalewajko et al. 2018) and showed that radiative cooling considerably reduces the efficiency of relativistic reconnection in producing nonthermal particles.

In this work, we present a systematic investigation of the role of strong IC losses on the plasmoid chain dynamics and nonthermal particle acceleration. We perform large-scale 2D and 3D PIC simulations of relativistic reconnection in e± plasma immersed in a dense radiation field, which exerts Compton drag on the particles. For simplicity, we assume that the radiation field is fixed during the simulation and composed of low-energy photons, so that Compton scattering occurs in the Thomson regime. The efficiency of Compton drag is determined by the radiation density, which is the main parameter of the problem, besides the magnetization parameter σ. We employ outflow boundary conditions along the reconnection exhausts; then, the plasmoid chain dynamics can be reliably studied in a quasi-steady state. A small guide field is included in the setup of the reconnection layer.

This paper is organized as follows. In Section 2 we describe the setup of our simulations and our choice for parameterization of the radiation density. In Section 3 we justify the physical parameters of our simulations as appropriate for the coronae of accreting black holes. In Section 4 we present our 2D results, focusing on the dynamics of the plasmoid chain and the spectrum of post-reconnection particles. We compare a representative run with strong IC losses with a no-cooling simulation (which otherwise has identical parameters) and then show how various properties of the reconnection system depend on the degree of IC cooling by varying it from zero to extreme values. In Section 5 we show that our main findings also hold in 3D simulations. Our conclusions and observational implications are summarized in Section 6.

2. Simulation Setup

We use the 3D electromagnetic PIC code TRISTAN-MP (Buneman 1993; Spitkovsky 2005) to study relativistic reconnection in electron–positron (pair) plasmas. Most of our simulations employ a 2D spatial domain in the x-y plane, except the full 3D simulations presented in Section 5. In all simulations we track all three components of the particle velocities and of the electromagnetic fields. Our 2D simulation setup closely parallels that in Sironi et al. (2016). It is summarized below for completeness. Details of the 3D setup are deferred to Section 5.

The reconnection layer is set up in Harris equilibrium, with the initial magnetic field ${{\boldsymbol{B}}}_{\mathrm{in}}=-{B}_{0}\,\hat{{\boldsymbol{x}}}\tanh (2\pi y/{\rm{\Delta }})$ reversing at y = 0 over a thickness Δ that will be specified below. The field strength is parameterized by the magnetization

Equation (1)

where ωc = eB0/mec is the Larmor frequency and ${\omega }_{{\rm{p}}}=\sqrt{4\pi {n}_{0}{e}^{2}/{m}_{e}}$ is the plasma frequency for the cold electron–positron plasma outside the layer (n0 is the number density, including both species), which is initialized with a small thermal spread of ${kT}/{m}_{e}{c}^{2}={10}^{-4}$. The Alfvén speed is related to the magnetization as ${v}_{A}/c=\sqrt{\sigma /(\sigma +1)}$. We focus on the regime σ ≫ 1 (i.e., vA/c ∼ 1) of relativistic reconnection and take σ = 10 as our fiducial case. In addition to the reversing field, we initialize a uniform component of "guide field" along z (i.e., aligned with the electric current in the reconnection layer) with strength Bg = 0.1 B0. The guide field provides pressure support to the cores of strongly cooled (and somewhat compressed) plasmoids. An investigation of the dependence of our results on σ and Bg/B0 is left for a future work.

Magnetic pressure outside the current sheet is initially balanced by particle pressure in the sheet, by adding a component of hot plasma that is overdense by a factor of 3 relative to the number density n0 of cold particles outside the layer. All the quantities presented in Sections 4 and 5 (e.g., thermodynamic properties of the layer and particle energy spectra) do not include this particle population set up in the current sheet, since its properties depend on arbitrary choices at initialization. The hot component is used only to set up the initial reconnection. After approximately one light-crossing time of the box, the artificially hot particles get ejected from the system. Then, reconnection proceeds in a quasi-steady state, losing memory of the initial hydrostatic state with the artificial hot component.

We trigger reconnection near the center of the computational domain, by removing the pressure of the hot particles initialized in the current sheet. This triggers a local collapse of the layer, which generates an X-point at the center. The initial perturbation results in the formation of two "reconnection fronts" that propagate away from the center along $\pm \hat{{\boldsymbol{x}}}$ (i.e., along the current layer), at roughly the Alfvén speed. We choose the thickness of the current sheet Δ to be large enough such that reconnection does not get spontaneously triggered anywhere else in the current layer, outside of the region in between the two reconnection fronts.4 In units of $\,{r}_{0,\mathrm{hot}}=\sqrt{\sigma }\,c/{\omega }_{{\rm{p}}}$ (which corresponds to the Larmor radius of particles with energy $\sigma {m}_{e}{c}^{2}$ in field B0)5 the thickness Δ is chosen to be Δ ∼ 30 ${r}_{0,\mathrm{hot}}$ ∼ 100 $c/{\omega }_{{\rm{p}}}$.

After one Alfvénic crossing time (=L/vA, where L is the half-length of our box along the x-direction), the two reconnection fronts reach the x boundaries of the computational domain. Here, we employ absorbing boundary conditions in the x-direction of the reconnection outflow, to mimic an open boundary where no information is able to propagate back inward (Daughton et al. 2006; Belyaev 2015; Cerutti et al. 2015). We refer to Sironi et al. (2016) for a discussion of our implementation of outflow boundaries. Along the y-direction of the reconnection inflow, we employ two "moving injectors" (receding from y = 0 along $\pm \hat{{\boldsymbol{y}}}$) and an expanding simulation box, a technique that we have extensively employed in our studies of relativistic shocks (Sironi & Spitkovsky 2009, 2011; Sironi et al. 2013) and relativistic reconnection (e.g., Sironi & Spitkovsky 2014; Sironi et al. 2016; Petropoulou & Sironi 2018). The two injectors constantly introduce fresh magnetized plasma into the simulation domain. This permits us to evolve the system as far as the computational resources allow, retaining all the regions that are in causal contact with the initial setup.

In the 2D simulations, we typically employ four particles per cell (including both species). We have checked that the results are weakly changed when using 16 particles per cell (extensive tests in the uncooled case have also been performed in Sironi et al. 2016). In order to reduce noise in the simulations, we filter the electric current deposited by the particles onto the grid, effectively mimicking the role of a larger number of particles per cell (Spitkovsky 2005; Belyaev 2015). The plasma skin depth $c/{\omega }_{{\rm{p}}}$ is resolved with 5 cells, so that the Larmor gyration period $2\pi /{\omega }_{{\rm{c}}}=2\pi /\sqrt{\sigma }\,{\omega }_{{\rm{p}}}$ is resolved with at least a few time steps (the numerical speed of light is 0.45 cells per time step).

In order to properly extrapolate our results to astrophysical scales, large-scale computational domains are essential. In the following, we shall call L the half-length of the computational domain along the x-direction, i.e., along the reconnection layer. In units of the Larmor radius of hot particles $\,{r}_{0,\mathrm{hot}}=\sqrt{\sigma }\,c/{\omega }_{{\rm{p}}}$, the half-length L of our fiducial 2D runs is L ≃ 1062 ${r}_{0,\mathrm{hot}}$. This corresponds to L/( $c/{\omega }_{{\rm{p}}}$) ≃ 3360, or equivalently 33,600 cells for the overall box size along the x-direction (so, for the full length 2 L).6 We evolve our simulations up to ∼5 L/c, corresponding to ∼185,000 time steps. This provides sufficient statistics to study the steady state of the system, which is established after ∼L/c, when the reconnection fronts reach the x boundaries. For the two fiducial cases discussed in Section 4, one without IC cooling losses and the other with strong cooling, we have also performed larger simulations, with L ≃ 2125 ${r}_{0,\mathrm{hot}}$ ≃ 6720 $c/{\omega }_{{\rm{p}}}$ (as well as a number of smaller simulations, which we do not report in this paper).

2.1. Compton Drag

In addition to the standard Lorentz force, the particles in our simulations are subject to "Compton drag" force,

Equation (2)

where ${{\boldsymbol{\beta }}}_{e}$ is the particle velocity in units of the speed of light c, ${\gamma }_{e}={(1-{\beta }_{e}^{2})}^{-1/2}$ is its Lorentz factor, σT is the Thomson cross section, and Urad is the radiation energy density. Equation (2) assumes that the radiation is isotropic in the simulation frame (the frame of the reconnection layer, in which the net Poynting flux vanishes). It also assumes that Urad is made of a large number of low-energy photons that scatter in the Thompson regime. This implies that the particle loses a small fraction of its energy in each scattering event, and the dynamic effect of scattering may be approximated as a continuous drag or friction opposing the particle motion.

One can define a characteristic Lorentz factor ${\gamma }_{\mathrm{cr}}$ ≫ 1 at which the maximum electric field in the reconnection layer, Erec = ηrecB0 (ηrec ∼ 0.1), is just sufficient to balance the Compton drag,

Equation (3)

Particle acceleration beyond ${\gamma }_{\mathrm{cr}}$ is suppressed by prohibitive IC losses. Instead of Urad, it is convenient to use ${\gamma }_{\mathrm{cr}}$ as the parameter controlling Compton drag; ${\gamma }_{\mathrm{cr}}=\infty $ corresponds to Urad = 0, while a low ${\gamma }_{\mathrm{cr}}$ corresponds to a large Urad and a strong drag. In this work, we run models with ${\gamma }_{\mathrm{cr}}$ as low as 11.3. This lowest ${\gamma }_{\mathrm{cr}}$ is comparable to σ = 10 chosen in our simulations, which represents a characteristic Lorentz factor for particle acceleration at X-points. Our fiducial model has ${\gamma }_{\mathrm{cr}}$ = 16. For comparison, we also run a benchmark model with ${\gamma }_{\mathrm{cr}}=\infty $ (no Compton drag).

A particle with a Lorentz factor γe is decelerated by Compton drag on the timescale ${t}_{\mathrm{IC}}={\gamma }_{e}{\beta }_{e}{m}_{e}c/{F}_{\mathrm{IC}}$, which may be rewritten as

Equation (4)

Our temporal resolution ${\rm{\Delta }}t=0.09\,{\omega }_{{\rm{p}}}^{-1}=0.09\,{\sigma }^{1/2}{\omega }_{c}^{-1}$ safely resolves this timescale even at the highest γe ∼ ${\gamma }_{\mathrm{cr}}$, since

Equation (5)

This constraint is easily satisfied as long as ${\gamma }_{\mathrm{cr}}\gtrsim \sigma \gg 1$, i.e., in the parameter regime we investigate here.

3. Plasma Conditions in Black Hole Coronae

In this section, we discuss the plasma conditions expected in coronae of black hole accretion disks (we refer to B17 for further details). Magnetic loops above the disk are sheared by differential rotation, generating current sheets. We focus on a current sheet of characteristic length L ∼ rg, where ${r}_{g}\,=2{{GM}}_{\mathrm{BH}}/{c}^{2}$ is the Schwarzschild radius of the black hole. The typical thickness of the current sheet is h ∼ ηrecL ∼ 0.1L, which corresponds to the width of the largest plasmoids in the reconnection layer (Sironi et al. 2016).

The reconnection process is initiated in the current sheet and dissipates magnetic energy, similar to solar flares but with a much higher power. This energy release is accompanied by creation of e± pairs inside and around the reconnection layer (see Figure 1 in B17). The pairs can strongly dominate the plasma density, depending on the compactness parameter of the magnetic flare. The pair plasma is kept at the Compton temperature of the radiation field until it flows into the reconnection layer. For simplicity in this work we neglect the ion component and assume a pure electron–positron composition with a given (pre-reconnection) density n0. We also idealize the problem by assuming a zero Compton temperature of the radiation field, so the plasma is cold before being heated by reconnection.

In the radiative regime, most of the dissipated magnetic energy is quickly converted to radiation, which escapes at speed ∼c, as long as the Thomson optical depth of the reconnection layer ${\tau }_{{\rm{T}}}={n}_{0}{\sigma }_{{\rm{T}}}h\sim {n}_{0}{\sigma }_{{\rm{T}}}{\eta }_{\mathrm{rec}}L$ is not much greater than unity (τT ∼ 1 is typically inferred for black hole coronae, and such optical depths may be sustained by pair creation in magnetic flares; see B17). Thus, energy conservation implies that the radiation density is Urad ∼ ηrecUB ∼ 0.1UB, where UB = B02/8π is the energy density of the reconnecting magnetic field (i.e., excluding the guide field, which suffers no dissipation).

The magnetization parameter σ can be expressed as

Equation (6)

where B is the magnetic compactness, defined as

Equation (7)

in analogy with the radiation compactness ${{\ell }}_{\mathrm{rad}}={U}_{\mathrm{rad}}{\sigma }_{{\rm{T}}}L/{m}_{e}{c}^{2}$.

Note that in the presence of ions, σ can be significantly reduced, as Equation (6) would include the ion mass density. Accretion disk coronae can have ion-dominated regions with σ < 1 and pair-dominated regions with σ ≫ 1. GRMHD simulations (e.g., Chatterjee et al. 2019; Jiang et al. 2019) demonstrate the existence of coronal regions (where magnetic pressure dominates the gas pressure) with σ < 1. Regions with σ ≫ 1 likely exist near the black hole spin axis, although they are harder to reproduce in MHD simulations, which employ density floors for numerical stability. Magnetic reconnection can generate heat in both regions, σ < 1 and σ > 1, with a rate scaling as (1 + 1/σ)−1/2. It is not known which region dominates the observed hard X-ray emission. In this paper, we focus on the pair-dominated plasma with σ ≫ 1; radiative reconnection in low-magnetization plasmas with ions will be studied elsewhere.

A characteristic compactness for black holes accreting at a significant fraction of the Eddington limit is B ∼ mp/me (B17). Then, using ηrec ∼ 0.1 and τT ∼ 1, Equation (6) gives σ ∼ 400. Then, reconnection occurs in the ultrarelativistic regime, and the fastest plasmoids are pushed to Lorentz factors $\sim \sqrt{\sigma }\gg 1$, i.e., reconnection can give ultrarelativistic bulk motions. The characteristic energy for particle acceleration at X-points is comparable to the mean energy per particle released in reconnection, γX ∼ σ/4 (the magnetic energy per particle is (σ/2)mec2, and another factor of 1/2 takes into account that only about half of the magnetic energy is dissipated; see Sironi et al. 2015).

Given that Urad ∼ ηrec UB, we can quantify the expected value of ${\gamma }_{\mathrm{cr}}$. Using the relation ${\sigma }_{{\rm{T}}}=8\pi {r}_{e}^{2}/3$, where ${r}_{e}={e}^{2}/{m}_{e}{c}^{2}$ is the classical electron radius, we find from Equation (3) that

Equation (8)

where we have substituted B ∼ mp/me. The high value of ${\gamma }_{\mathrm{cr}}$ ≫ γX ∼ σ/4 implies that particle acceleration at X-points is not impeded by Compton drag. The same condition can be expressed by comparing the IC timescale (Equation (4)) with the timescale tX for particle acceleration at X-points,

Equation (9)

On the other hand, the IC cooling timescale is much shorter than the timescale for plasma advection along the layer tadv ∼ L/vA ∼ L/c. Their ratio is

Equation (10)

This condition must be satisfied in our simulations for any ${\gamma }_{e}$, including ${\gamma }_{e}$ ∼ 1, which requires a large size L of the computational box,

Equation (11)

In summary, in order to mirror the conditions expected in black hole coronae, our simulations need to have (i) σ ≫ 1, and our choice of σ = 10 is then suitable; (ii) ${\gamma }_{\mathrm{cr}}$ ≫ γX ∼ σ/4, which is satisfied by all our simulations; and (iii) $L/(\,c/{\omega }_{{\rm{p}}})\gg {\gamma }_{\mathrm{cr}}^{2}/({\eta }_{\mathrm{rec}}\sqrt{\sigma })$. This is achieved, e.g., in our fiducial 2D model by choosing ${\gamma }_{\mathrm{cr}}$ = 16 and L/( $c/{\omega }_{{\rm{p}}}$) ≃ 3360. Thus, although our runs do not have the huge ${\gamma }_{\mathrm{cr}}$ ∼ 104 expected in black hole coronae, they satisfy the required hierarchy of relevant scales and energies.

4. Results

In this section, we present the results of our 2D simulations. In Sections 4.14.3, we compare two simulations: one with no IC losses (${\gamma }_{\mathrm{cr}}=\infty $) and the other with strong losses (${\gamma }_{\mathrm{cr}}$ = 16). We first present the overall flow structure in Section 4.1, then discuss the effect of IC drag on the plasmoid chain dynamics in Section 4.2, and finally explore the particle energy spectrum in Section 4.3. Section 4.4 presents the results of a suite of 2D simulations aimed at demonstrating how basic parameters of the reconnection layer depend on ${\gamma }_{\mathrm{cr}}$ for a wide range between ${\gamma }_{\mathrm{cr}}$ = 11.3 ∼ σ and ${\gamma }_{\mathrm{cr}}=\infty $.

4.1. Flow Structure

The overall structure of the reconnection flow with strong IC losses (${\gamma }_{\mathrm{cr}}$ = 16) is presented in Figure 1. It should be compared with Figure 2, which shows the simulation without IC losses, with otherwise identical parameters. One evident difference is the low internal energy in the radiative case. This is expected, as heat is lost to IC radiation.

Figure 1.

Figure 1. 2D structure of the reconnection layer at time t = 3.4(L/c) in our fiducial run with strong IC losses (γcr = 16). We show the region $| y| /L\lt 0.125$ where reconnection occurs (the actual extent of the computational box along y grows with time as described in Section 2). (a) Magnetic energy density B2/8π normalized by the initial plasma rest-mass energy, ${\epsilon }_{{\rm{B}}}={B}^{2}/8\pi {n}_{0}{m}_{e}{c}^{2}$, with overplotted magnetic field lines. The initial epsilonB far from the reconnection layer is epsilonB ≈ σ/2 = 5. (b) Particle number density n, in units of the number density n0 far from the reconnection layer. (c) Local average $\langle {\gamma }_{e}^{2}{\beta }_{e}^{2}\rangle =\langle {\gamma }_{e}^{2}\rangle -1$, obtained for each cell by averaging over the particles in the neighboring 5 × 5 cells. This quantity is proportional to the local total IC power radiated per particle. (d) γ2 − 1, where γ is the plasma bulk Lorentz factor, as defined in the text. This quantity is proportional to the local bulk IC power per particle. (e) The difference $\langle {\gamma }_{e}^{2}\rangle -{\gamma }^{2}$ represents the contribution from internal particle motions. One can see that the IC power radiated by isolated plasmoids is dominated by their bulk motion, whereas internal particle motions dominate in the thin regions near the midplane X-points and in the transient reconnection layers formed between merging plasmoids.

Standard image High-resolution image

The magnetic field structure displays a striking similarity in the two runs. The opening angle of the magnetic field lines in the inflow region is closely related to the inflow speed (Lyubarsky 2005), and so the similarity of these angles in Figures 1 and 2 indicates that the reconnection rate weakly depends on ${\gamma }_{\mathrm{cr}}$. We find that the time-averaged reconnection rate is ηrec = vrec/c ≈ 0.122 in the nonradiative simulation and ηrec/c ≈ 0.135 in the radiative simulation. The weak dependence of ηrec on the degree of radiative losses is further demonstrated in Section 4.4.

Figure 2.

Figure 2. Same as Figure 1, but for a simulation with no IC losses (${\gamma }_{\mathrm{cr}}=\infty $). Apart from ${\gamma }_{\mathrm{cr}}$, all parameters of the simulation are identical to that in Figure 1, and the snapshot is taken at the same time t = 3.4(L/c). Internal particle motions now dominate over bulk motions, because particles keep the heat received in reconnection rather than radiate it. However, the overall shape of the plasmoid chain and the reconnection rate are remarkably similar to the run with strong IC losses.

Standard image High-resolution image

In both radiative and nonradiative simulations, the reconnection layer fragments into a chain of plasmoids of various sizes, which appear as strongly magnetized structures. At the centers of plasmoids, the parameter ${\epsilon }_{{\rm{B}}}\equiv {B}^{2}/8\pi {n}_{0}{m}_{e}{c}^{2}$ (the ratio of magnetic energy density to upstream rest-mass energy density) reaches ∼100, regardless of the level of radiative losses. This should be compared with its initial value in the inflow region, epsilonB ≈ σ/2 = 5. The magnetic field B in the plasmoid core is compressed by a factor of B/B0 ∼ 4. When averaged over the plasmoid volume, the mean magnetic energy density in the plasmoid rest frame is roughly twice that in the upstream (as shown by Sironi et al. 2016 for cases with zero guide field Bg = 0). The similarity of the magnetic field structure in Figures 1 and 2 demonstrates that the plasmoid chain formation is controlled by magnetic stresses, and plasma pressure is not crucial even in the hot (nonradiative) regime.

Plasma density structure, on the other hand, exhibits significant differences between the two runs (compare panels (b) of Figures 1 and 2). In the radiative simulation, the plasmoids are strongly asymmetric: the plasma trapped in magnetic islands tends to pile up in the back of the moving island, leaving a cavity at the front (see, e.g., the plasmoids at x/L ∼ −0.8 and at x/L ∼ 0.9 in Figure 1(b)). This effect is caused by Compton drag on the plasma particles while most of the island inertia is in the magnetic field. Since plasma can slide along the magnetic field lines, Compton drag partially succeeds in resisting the plasma motion and keeps it in the back of the island. This is analogous to a fast-moving car with open top and without a wind shield, exposing passengers to air resistance. The air drag (Compton drag) pushes the passengers (plasma particles) toward the back of the car (magnetic island).

The biggest change in the radiative case is that most of the power released by reconnection is radiated away rather than stored in the plasma. In particular, almost all internal energy received by the plasma is lost. As a result, the IC emission is generated predominantly by the bulk motions of cooled plasmoids (especially the large ones). The approximate equality of the total and bulk IC losses is demonstrated by the nearly identical colors of large plasmoids in panels (c) and (d) of Figure 1 (see, e.g., the plasmoid at x/L ∼ −0.8). This confirms the picture described in B17: most of the reconnection energy goes into work by magnetic tension forces to move cold plasmoids against Compton drag. As a result, most of the released energy converts directly to radiation through bulk Comptonization.

A smaller fraction of the dissipated reconnection power goes into impulsive acceleration of individual particles, followed by their quick cooling.7 Individual particles are energized outside the plasmoids, in particular at X-points or in fast outflows emanating from the X-points. This is clearly observed in the thin regions of the reconnection layer, where the freshly energized plasma has not cooled yet (see, e.g., the thin red layer at −0.2 ≲ x/L ≲ 0.1 in Figure 1(e)). The same happens at the interfaces between merging plasmoids where transient "secondary" current sheets are formed perpendicular to the main current sheet (see the vertical red layer at x/L ∼ −0.6 or x/L ∼ 0.3). Here particles are actively heated and accelerated just like in the primary current sheet.

4.2. Bulk Motions in the Reconnection Layer

We now describe in more detail the plasma bulk motions in the radiative reconnection layer. We measure the bulk motion (for each species, electron and positron) as follows. For each cell, we compute the average particle velocity including all the particles in the neighboring 5 × 5 cells: ${\boldsymbol{\beta }}=\langle {{\boldsymbol{\beta }}}_{e}\rangle $. In the frame moving with velocity ${\boldsymbol{\beta }}$, the plasma stress energy tensor has a vanishing energy flux (e.g., Rowan et al. 2019), and we define this frame as the comoving frame of the plasma (see, e.g., Zhdankin et al. 2018 for alternative definitions). The corresponding bulk Lorentz factor is given by γ = (1 − βb2)−1/2.

We have measured the statistics of bulk motions in the reconnection layer using the positron component as a proxy. We include only particles inside the reconnection layer. This layer is defined based on the degree of mixing between the two opposite vertical inflows toward the midplane y = 0 (see, e.g., Rowan et al. 2017; Ball et al. 2018). Particles coming from above (y > 0) and below (y < 0) are tagged as two different populations, and the reconnection layer is defined as the region where both populations contribute at least 1% to the local plasma density. The results shown below are nearly insensitive to the exact value of the mixing threshold, as long as it is significantly above zero.

There is some y-component of bulk motions, which is produced by vertical secondary reconnection layers between merging plasmoids.8 However, the dominant component of bulk motions is along the x-axis (Figure 3), and below we focus on the x-component. It strongly varies in space and time. The scatter of these variations may be described by an effective temperature, which influences the hard X-ray spectrum produced by the reconnection layer through Comptonization (B17).

Figure 3.

Figure 3. Density-weighted distribution of positron bulk momenta, from the same simulations as in Figure 1 (${\gamma }_{\mathrm{cr}}$ = 16; left panel) and Figure 2 (${\gamma }_{\mathrm{cr}}=\infty ;$ right panel). A similar distribution is found for the electron component (not shown in the figure). The statistics of bulk 4-velocities ${\bf{u}}=\gamma {\boldsymbol{\beta }}$ was accumulated by averaging snapshots over the time period 1.5 ≲ ct/L ≲ 5. Dashed white lines indicate the dimensionless 4-velocity $| {u}_{x,y}| =\sqrt{\sigma }$ expected for the fastest outflows in relativistic reconnection.

Standard image High-resolution image

Figure 4 presents the distribution of ux = γ${\beta }_{x}$ measured at each vertical slice of the reconnection layer during an extended time period 1.5 ≲ ct/L ≲ 5. The result depends on the slice position x. We observe that the horizontal motions have a regular component (represented by the median or the mean expectation of ux(x)) and a strong stochastic component. The regular component represents the plasma outflow away from x = 0 (the center of the layer). By symmetry, the outflow speed is zero at x = 0 and increases toward the edges x = ±L, approaching $| {u}_{x}| \approx 2$ in the nonradiative model and $| {u}_{x}| \approx 1$ in the radiative model. This regular component is below $| {u}_{x}| \approx \sqrt{\sigma }\approx 3$ predicted by analytical models (Lyubarsky 2005). The strong reduction of the mean outflow speed in radiative reconnection is the result of Compton drag.

Figure 4.

Figure 4. Scatter of positron bulk motions in the x-direction (the direction of the reconnection outflow) from the same simulations as in Figure 1 (${\gamma }_{\mathrm{cr}}$ = 16; top) and Figure 2 (${\gamma }_{\mathrm{cr}}=\infty ;$ bottom). The scatter results from spatial and temporal variations of the plasma flow in the reconnection region. The plot was obtained by averaging 70 snapshots over the extended time period 1.5 ≲ ct/L ≲ 5 to convey a representative quasi-steady-state picture. In each snapshot, the positron bulk motion was measured using the x-component of the positron bulk 4-velocity ux. The measurement was made in each grid cell of the reconnection region (as defined in the text) and assigned weight proportional to the local positron density. At a fixed x, the result varies in time and also depends on the y position of the cell. This variability generates the distribution of ux shown in the figure. The solid black curve shows the median value of ux as a function x (the mean expectation of ux gives nearly the same curve), and the dashed black curves indicate the 10th and 90th percentiles of the ux distribution. The horizontal dashed white lines indicate $\pm \sqrt{\sigma }$, the value of 4-velocity expected for the fastest reconnection outflows.

Standard image High-resolution image

Superimposed on the average outflow are substantial stochastic motions. They are particularly strong for small plasmoids, which are capable of approaching the maximum $| {u}_{x}| \sim \sqrt{\sigma }$ even for the radiative case. About 10% of the time, $| {u}_{x}| $ exceeds 2. A more typical variation is comparable to the mean expectation value. The stochastic bulk motions are caused by continuing tearing instabilities forming new plasmoids and their subsequent mergers. Large old plasmoids tend to attract and accrete small young plasmoids, and so some young plasmoids are pulled back by a large plasmoid behind them. As a result, a small fraction of plasma motions are directed toward the center x = 0, opposite to the mean expectation.

The observed distributions of 4-velocity ux do not, however, provide a useful physical picture of the Compton drag effect on the reconnection chain. They are ignorant of the fact that most of reconnection energy resides in distinct, isolated plasmoids of various sizes w, as clearly seen in panel (a) of Figure 1. Therefore, we extend our analysis of the reconnection layer by examining the motions of individual plasmoids. Similar to Sironi et al. (2016), we identify the individual magnetic islands using the magnetic vector potential. The islands are typically found to have a quasi-elliptical shape, with a transverse size w and a comoving length of ∼3w/2. The plasmoid 4-velocity $u$ is defined as the bulk ux of the center of the magnetic island (the peak of the vector potential).

The motion of a plasmoid is governed by two competing forces: it is accelerated by magnetic tension and decelerated by Compton drag. The outcome of this competition depends on the plasmoid size. B17 argued that the magnetic tension forces in radiative reconnection should be approximately the same as in nonradiative reconnection and estimated the chain properties based on this conjecture. The basic idea of this method is as follows. In the absence of Compton drag, the measured plasmoid motions directly reflect the magnetic forces fB accelerating them. Thus, the nonradiative simulation gives ${f}_{B}(w,{t}_{\mathrm{age}}^{{\prime} })$ acting on plasmoids of various sizes w and proper ages ${t}_{\mathrm{age}}^{{\prime} }$.9 If ${f}_{B}(w,{t}_{\mathrm{age}}^{{\prime} })$ indeed remains the same for radiative reconnection, one can simply add Compton drag fdrag and evaluate the dynamics of plasmoids subjected to two known forces, fB and fdrag. Below we test this method against the results of our numerical simulations.

Previous analysis of nonradiative reconnection simulations demonstrated a close relation between ${t}_{\mathrm{age}}^{{\prime} }$, w, and the plasmoid 4-velocity $u$ (Sironi et al. 2016). This relation is described by the empirical formula,

Equation (12)

It shows that plasmoids with ${\eta }_{\mathrm{rec}}{{ct}}_{\mathrm{age}}^{{\prime} }/w\gt 1$ approach the maximum, saturated 4-velocity $u\approx \sqrt{\sigma }$. It also shows that before the saturation occurs, i.e., when $u$ is well below $\sqrt{\sigma }$, plasmoids grow proportionally to their proper ages, $w\approx {\eta }_{\mathrm{rec}}{{ct}}_{\mathrm{age}}^{{\prime} }$, and the growth slows down after the saturation of $u$.

A large part of the plasmoid momentum is carried by the magnetic field (even in a hot, nonradiative chain). After averaging over the plasmoid volume, its momentum density is given by

Equation (13)

where we used the relation ${U}_{B}\sim {\gamma }^{2}{B}_{0}^{2}/4\pi $ measured by Sironi et al. (2016). Plasmoids acquire this momentum during time ${t}_{\mathrm{age}}=\gamma {t}_{\mathrm{age}}^{{\prime} }$. This implies an average force per unit volume fB ∼ ppl/tage and gives

Equation (14)

Equations (12) and (14) allow one to exclude ${t}_{\mathrm{age}}^{{\prime} }$ and find fB acting on a plasmoid with given w and $u$.

We can now compare fB with the Compton drag force fdrag (B17),

Equation (15)

Here npl is the mean density of the plasmoid, and ${\gamma }_{e}^{{\prime} }$ is the characteristic Lorentz factor of particles measured in the plasmoid rest frame; ${\gamma }_{e}^{{\prime} }-1$ is a measure of heat. For most plasmoids (except the smallest ones), the drag timescale tdrag = ppl/fdrag satisfies the following conditions (see B17):

Equation (16)

This implies that the plasmoids are efficiently cooled to ${\gamma }_{e}^{{\prime} }\approx 1$ (we check the accuracy of this condition for our simulations below). Then, one finds

Equation (17)

where ${n}_{\mathrm{pl}}^{{\prime} }\approx {n}_{\mathrm{pl}}/\gamma $ is the mean proper density of the plasmoid.

The main predicted effect of Compton drag is that the plasmoid acceleration should cease when fdrag = fB, and so plasmoids should be prevented from entering the regime where fdrag > fB. The boundary of this "avoidance zone" on the w-u plane can be found by setting ${f}_{\mathrm{drag}}/{f}_{B}=1$ in Equation (17), then expressing ${t}_{\mathrm{age}}^{{\prime} }$ from this equation, and substituting the result into Equation (12). This gives

Equation (18)

Using the compression factor ${n}_{\mathrm{pl}}^{{\prime} }/{n}_{0}\approx 6$ measured in the simulation with ${\gamma }_{\mathrm{cr}}$ = 16 and L = 3360 $c/{\omega }_{{\rm{p}}}$, one can solve Equation (18) for $u$ for plasmoids of a given size w/L. The solution represents the expected boundary of the avoidance zone on the w-u plane. Plasmoids are expected to congregate near this boundary, and we can quantitatively test this expectation by measuring the plasmoid speeds in the simulation.

Figure 5 shows the 4-velocity $u$ of plasmoids at the end of their lives, before they exit the computational box or merge with another plasmoid. More exactly, the measurement of $u$ was taken when the plasmoid reached its maximum transverse size wmax, and the results are presented as a function of wmax. Figure 5 shows the results for nonradiative (${\gamma }_{\mathrm{cr}}=\infty $) and radiative (${\gamma }_{\mathrm{cr}}$ = 16) reconnection.

Figure 5.

Figure 5. Outcome of plasmoid acceleration as a function of its size. Plasmoid 4-velocities u were measured when the plasmoid reached its maximum size wmax, i.e., typically at the end of its life in the simulation—either at a merger or when the plasmoid exits the computational box. Blue and red points show the results in the radiative (${\gamma }_{\mathrm{cr}}$ = 16) and nonradiative (${\gamma }_{\mathrm{cr}}=\infty $) simulation, respectively. The red curve shows the analytical formula in Equation (19). The avoidance zone expected in the radiative regime and described by Equation (18) is shaded in blue.

Standard image High-resolution image

In the nonradiative regime, the measurements are consistent with Equation (12) when we substitute10 $\gamma \beta {{ct}}_{\mathrm{age}}^{{\prime} }\sim L$, which gives

Equation (19)

In agreement with previous simulations (Sironi et al. 2016), we observe the saturation of $u\approx \sqrt{\sigma }$ for small plasmoids, with ${w}_{\max }/L\lesssim ({\eta }_{\mathrm{rec}}/\sqrt{\sigma })\approx 0.03$.

A similar saturation should hold in the radiative regime; however, now there is an additional limiting factor—Compton drag. For radiative reconnection, we compare the measured $u({w}_{\max })$ with Equation (18) and observe a good agreement in the range $0.003\lt {w}_{\max }/L\lt 0.1$. This quantitative test confirms the conjecture that magnetic tension forces are approximately the same in radiative and nonradiative regimes. Deviations are observed for the smallest plasmoids and caused by ${\gamma }_{e}^{{\prime} }\gt 1$.

The deviation of ${\gamma }_{e}^{{\prime} }$ from unity is small for most plasmoids in radiative reconnection, because radiative losses dramatically reduce the plasma internal energy. The cooling timescale measured in the plasmoid rest frame is given by (B17)

Equation (20)

where ${U}_{\mathrm{rad}}^{{\prime} }\sim {\gamma }^{2}{U}_{\mathrm{rad}}$ from the transformation of Urad to the plasmoid frame. Cooling occurs faster than Compton drag. The drag effect during the cooling timescale ${t}_{\mathrm{IC}}=\gamma {t}_{\mathrm{IC}}^{{\prime} }$ may be expressed as

Equation (21)

Nonradiative reconnection with σ ≫ 1 would create plasmoids with σpl ∼ 1. By contrast, radiative relativistic reconnection gives σpl ≫ 1 for almost all plasmoids, except the young/small ones that have not cooled yet but have already reached $\gamma \approx \sqrt{\sigma }$. The small plasmoids must cool to ${\gamma }_{e}^{{\prime} }\approx 1$ when their proper ages reach ${t}_{\mathrm{IC}}^{{\prime} }$ evaluated with ${\gamma }_{e}^{{\prime} }\sim 1$ and $\gamma \sim \sqrt{\sigma }$. This condition defines the "cooling age." The corresponding characteristic "cooling size" of plasmoids is given by

Equation (22)

In particular, in our fiducial model with ${\gamma }_{\mathrm{cr}}=16$ and $L/(\,c/{\omega }_{{\rm{p}}})\sim 3360$, we find wc/L ∼ 0.003. Thus, all plasmoids with w/L ≫ 0.003 should be cooled to ${\gamma }_{e}^{{\prime} }\sim 1$.

This analytical estimate is confirmed in Figure 6, where we plot the mean internal energy per particle in plasmoids as a function of the plasmoid transverse size. In the radiative case (blue curve), we observe $\langle {\gamma }_{e}^{{\prime} }\rangle \approx 1.5$ for plasmoids with w/L ≲ 0.003. This explains why the bulk motions of small plasmoids (blue points in Figure 5) stay somewhat below the prediction of Equation (18), which assumed ${\gamma }_{e}^{{\prime} }$ = 1. For plasmoids with w/L ≳ 0.003, $\langle {\gamma }_{e}^{{\prime} }\rangle $ drops toward unity, in agreement with the estimate in Equation (22). We have also verified that wc remains consistent with Equation (22) in a larger simulation (with doubled L), i.e., ${w}_{c}/(\,c/{\omega }_{{\rm{p}}})$ remained unchanged.

Figure 6.

Figure 6. Mean internal energy per particle (in units of ${m}_{{ec}}^{2}$) in plasmoids, as a function of the plasmoid transverse size w, in the radiative (blue) and nonradiative (red) simulations. The filled circles connected by the lines show the median values, whereas the error bars show the 10th and 90th percentiles.

Standard image High-resolution image

In contrast, for nonradiative reconnection we always find $\langle {\gamma }_{e}^{{\prime} }\rangle \gtrsim 2$, with a tendency for hotter particles in larger plasmoids (red curve in Figure 6).

4.3. High-energy Particles

In nonradiative reconnection, a large fraction of the released magnetic energy becomes stored in relativistic, ${\gamma }_{e}^{{\prime} }\gg 1$, quasi-isotropic internal motions of particles in the plasma rest frame. This fact is demonstrated in Figure 7. The figure compares the particle distributions over the bulk motion energy $(\gamma -1){m}_{e}{c}^{2}$ and the actual particle energy $({\gamma }_{e}-1){m}_{e}{c}^{2}$. One can see that the distribution ${dN}/d\mathrm{log}({\gamma }_{e}-1)$ is strongly shifted toward high energies compared with ${dN}/d\mathrm{log}(\gamma -1)$. This reflects the fact that nonradiative reconnection heats the plasma to ${\gamma }_{e}^{{\prime} }\gg 1$, and so particles gain ${\gamma }_{e}\gg \gamma $ (note that ${\gamma }_{e}\sim \gamma {\gamma }_{e}^{{\prime} }$). The hot distribution peaks at ${\gamma }_{e}\simeq 6\sim 0.5\sigma $. The broad bump at ${\gamma }_{e}\sim 0.5\sigma $ is dominated by plasma between plasmoids, in particular by particles in the outflows from X-points (they are energized by the pickup process discussed in detail below). The slope $p\equiv d\mathrm{log}N/d\mathrm{log}{\gamma }_{e}$ on the high-energy side of this bump is steep, making the appearance of a power-law tail with $p\sim 3.5$.11

Figure 7.

Figure 7. Particle energy distribution in the nonradiative simulation (${\gamma }_{\mathrm{cr}}=\infty $), time-averaged in the interval 1.5 ≲ ct/L ≲ 5. The solid red curve shows the particle distribution ${dN}/d\mathrm{log}({\gamma }_{e}-1)$ in the reconnection region (as defined in Section 4.2). The dotted red curve shows the distribution in the entire computational box, including the cold inflow region, where particles move with speed v ≈ ηrecc and form the low-energy peak at ${\gamma }_{e}-1\sim {10}^{-2}$. The dashed orange curve shows the particle distribution in the reconnection region when accounting only for bulk kinetic energy (i.e., we plot ${dN}/d\mathrm{log}(\gamma -1)$ as a function of γ − 1). All three distributions are normalized to the total number of particles in the reconnection region.

Standard image High-resolution image

Strong IC losses change the particle distribution (Figure 8). There are three main changes. First, nonthermal particle acceleration to ${\gamma }_{e}\gg \sigma $ is suppressed (see also Werner et al. 2019). Second, the peak of ${dN}/d\mathrm{log}({\gamma }_{e}-1)$ is shifted from ${\gamma }_{e}\sim 0.5\sigma $ to a much lower, mildly relativistic value ${\gamma }_{e}\sim 1.3$. This change reflects the fact that the plasma loses most of its energy to radiation. Third, the distribution shape at ${\gamma }_{e}$ < 2 becomes controlled by bulk motions, i.e., the particles populating the peak are cold, with ${\gamma }_{e}^{{\prime} }$ ≈ 1 and ${\gamma }_{e}\approx \gamma $.

Figure 8.

Figure 8. Particle distribution in the radiative simulation (${\gamma }_{\mathrm{cr}}$ = 16), time-averaged in the interval 1.5 ≲ ct/L ≲ 5. Similar to Figure 7, we show ${dN}/d\mathrm{log}({\gamma }_{e}-1)$ in the reconnection region (solid blue) and in the entire computational box (dashed blue). The dashed cyan curve shows the particle distribution in the reconnection region when accounting only for bulk kinetic energy (i.e., we plot ${dN}/d\mathrm{log}(\gamma -1)$ as a function of γ − 1). For comparison, the solid black curve shows a Maxwellian distribution with temperature kTb = 100 keV. The dashed green curve is the difference between the solid blue and the dashed cyan curves; it represents the part of the energy distribution that is not accounted for by bulk motions. The two vertical dashed lines approximately show the division into three energy regions where different acceleration processes dominate: bulk acceleration by magnetic tension forces (${\gamma }_{e}$ ≲ 2), particle pickup by outflows from X-points (2 < ${\gamma }_{e}$ < 6), and X-point acceleration (${\gamma }_{e}$ ≳ 6).

Standard image High-resolution image

The bulk motion distribution ${dN}/d\mathrm{log}(\gamma -1)$ steeply declines above γ ∼ 2, and its shape can be approximately described as Maxwellian, with an effective "bulk temperature" ${{kT}}_{{\rm{b}}}\approx 100\,\mathrm{keV}$. At Lorentz factors ${\gamma }_{e}$ ∼ 2 (kinetic energies $\sim {m}_{e}{c}^{2}=511$ keV), the distribution ${dN}/d\mathrm{log}({\gamma }_{e}-1)$ strongly deviates from ${dN}/d\mathrm{log}(\gamma -1)$ and forms an extended tail with a slope $p\sim -4.3$. Thus, the distribution can be rather cleanly separated into two distinct components: the main 100 keV peak shaped by the bulk motions of cold plasma and the high-energy component dominated by individual motions of nonthermal particles. The high-energy component extends from ${\gamma }_{e}$ ∼ 2 to ${\gamma }_{e}$ of a few σ.12

The cooling time of energetic particles is short, and so the observed extension of dN/d log(${\gamma }_{e}$ − 1) from 500 keV to ∼10 MeV must be sustained by a process of nonthermal particle injection in the reconnection layer. Indeed, we observe that nearly impulsive acceleration of particles operates at the X-points and outflows from the X-points, as described below.

The high-energy component receives a significant fraction of the total power dissipated by magnetic reconnection. We can measure this fraction ${f}_{\mathrm{HE}}$ using the fact that the received power is immediately radiated away, with the emission rate proportional to ${\gamma }_{e}^{2}{\beta }_{e}^{2}$. We have computed the ratio of the power emitted by the high-energy particles (${\gamma }_{e}$ > 2) and the net power emitted by all particles in the reconnection region. This ratio is ${f}_{\mathrm{HE}}\approx 35$% for our fiducial model shown in Figure 8 (${\gamma }_{\mathrm{cr}}=16$ and L/( $c/{\omega }_{{\rm{p}}}$) ≃ 3360). We have also explored how ${f}_{\mathrm{HE}}$ depends on the simulation parameters and found that it decreases to ${f}_{\mathrm{HE}}$ ≈ 27% for a larger computational box L/( $c/{\omega }_{{\rm{p}}}$) ≃ 6720. This decrease is moderate, indicating convergence with increasing box size. We obtained similar values for ${f}_{\mathrm{HE}}$ in models with yet stronger IC cooling (lower ${\gamma }_{\mathrm{cr}}$). In particular, in simulations with ${\gamma }_{\mathrm{cr}}$ = 11.3 we find ${f}_{\mathrm{HE}}$ ≈ 26% for L/( $c/{\omega }_{{\rm{p}}}$) ≃ 3360 and ${f}_{\mathrm{HE}}$ ≈ 22% for L/( $c/{\omega }_{{\rm{p}}}$) ≃ 6720. We conclude that ${f}_{\mathrm{HE}}$ ∼ 20% in the strongly radiative regime.

The injection of high-energy particles is important for the radiation spectrum emitted by the reconnection layer. Therefore, we further investigated the mechanism of nearly impulsive acceleration observed in our simulations. We have analyzed the histories of individual particles and identified two mechanisms generating the high-energy component. (1) The distribution at γe ≳ σ = 10 is dominated by particles ejected from the X-points by the accelerating electric field E ∼ 0.1B0. This mechanism was studied in detail in previous works (e.g., Zenitani & Hoshino 2001). (2) At intermediate energies, 2 ≲ ${\gamma }_{e}$ ≲ 6, the distribution is dominated by particles located between the X-points and the neighboring plasmoids. We call these regions "unstructured outflows" from the X-points. We found that the energetic particles appearing in the outflows are injected by a special "pickup" process described below.

When two opposite horizontal magnetic field lines (above and below the reconnection layer) become connected at the X-point, they form a single field line with a cusp. Magnetic tension causes the cusp to snap horizontally toward the neighboring plasmoid. The snapping motion occurs with an increasing Lorentz factor γ, and the moving cusp "picks up" the plasma particles on the field line. This process resembles collision of a moving wall with a static particle, which gives the particle γe ∼ γ2 in the lab frame (${\gamma }_{e}^{{\prime} }\sim \gamma $ in the wall rest frame). The fastest unstructured outflows have $\gamma \approx \sqrt{\sigma }$ and can energize particles up to γmax ∼ σ. In reality, this is a conservative upper limit, since most of the unstructured outflows move slower, especially in the presence of Compton drag. As a result, the kicks of particles by unstructured outflows dominate the observed distribution only up to γe ∼ 6 < σ = 10.

The picked-up particles start to gyrate about ${\boldsymbol{B}}$ with Lorentz factor ${\gamma }_{e}$' ≈ γ in the rest frame of the moving cusp (note that the magnetic field at the cusp point is perpendicular to its velocity, Bx = 0). The gyrating particle possesses a significant oscillating z-momentum in the lab frame, in addition to the x-momentum associated with the cusp motion. This is a distinct feature of particles energized by pickup,13 as well as particles accelerated at X-points. The particle loses its gyration energy on the IC cooling timescale, and its oscillating z-momentum eventually vanishes. Then, the cold particle simply follows the field-line drift along x, with the bulk speed of the unstructured outflow.

The pickup process is quick—it occurs on the gyration timescale $\sim \gamma /({\eta }_{\mathrm{rec}}\sqrt{\sigma }\,{\omega }_{{\rm{p}}})$, comparable to the acceleration timescale at X-points. The sharp rise in energy, with the development of a significant z-momentum, clearly identifies the acceleration episodes when analyzing the histories of particles. For a particle with a Lorentz factor ${\gamma }_{e}$ at a given moment of time, we examine its recent history during one cooling timescale tIC(${\gamma }_{e}$) and measure the evolution of d ${\gamma }_{e}$/dt during the past Δt = tIC. Its median value represents a characteristic gain rate d ${\gamma }_{e}$/dt in the recent past. It is expected to be high for particles with high ${\gamma }_{e}$, as otherwise they would have cooled. In addition, we measure the ratio $| {u}_{e,z}/{u}_{e,x}| $ for the particle during the same time interval Δt = tIC and compute its median value. Particles that have just experienced the pickup or X-point acceleration are expected to have a high $| {u}_{e,z}/{u}_{e,x}| $, as their gyration involves a large z-component. In order to accumulate good statistics, we repeat the measurement of d ${\gamma }_{e}$/dt and $| {u}_{e,z}/{u}_{e,x}| $ for a large number of particles (roughly 10 million), with a temporal cadence of $18\,{\omega }_{{\rm{p}}}^{-1}$.

The results are presented in Figure 9. In the ${\gamma }_{e}-d{\gamma }_{e}/{dt}$ plot we observe a transition at ${\gamma }_{e}\sim 2$, from particles with a low recent acceleration $d{\gamma }_{e}/d({\omega }_{{\rm{p}}}t)\sim {10}^{-3}$ to particles with a strong recent acceleration. A similar transition is seen in the plot of $| {u}_{e,z}/{u}_{e,x}| $. Particles with ${\gamma }_{e}\lt 2$ have a dominant ${u}_{e,x}$, as expected for cold particles moving with the local bulk velocity of the unstructured outflow. The ratio $| {u}_{e,z}/{u}_{e,x}| $ quickly increases at ${\gamma }_{e}\gt 2$ to values comparable to unity. In particular, $| {u}_{e,z}/{u}_{e,x}| \sim 0.5$ for $3\lesssim {\gamma }_{e}\lesssim 7$ and $| {u}_{e,z}/{u}_{e,x}| \sim 1$ for ${\gamma }_{e}\gtrsim 7$. These diagnostics provide further confidence in the separation of the energy distribution in the reconnection region into two distinct components: bulk-motion-dominated (${\gamma }_{e}\lesssim 2$) and high-energy particle injection (${\gamma }_{e}\gtrsim 2$).

Figure 9.

Figure 9. Diagnostics of particle acceleration in the radiative simulation (${\gamma }_{\mathrm{cr}}$ = 16). Panels (a) and (b) show 2D histograms of the tracked particles, time-averaged in the interval 1.5 ≲ ct/L ≲ 5. In panel (a), the vertical axis represents the median acceleration rate e/d(ωpt) during the recent time interval Δt equal to the IC cooling time of the particle tIC(γe). In panel (b), the vertical axis represents the median value of $| {u}_{e,z}/{u}_{e,x}| $ during the recent Δt = tIC(γe). The histogram values are shown by color. They are normalized separately for each bin δ ln(${\gamma }_{e}$ − 1) (the integral over each vertical slice equals unity); the green curve shows the median of the histogram. The vertical dashed line is the boundary (${\gamma }_{e}$ > 2) of one of the two high-energy intervals discussed in the text. (c) Contributions of particles with strong recent accelerations (solid blue curve) or high z-momenta (dashed blue curve) to the particle energy distribution in the reconnection layer (black curve). Strong acceleration/high z-momentum is defined as the region above the horizontal dashed line in panel (a)/(b). One can see that both conditions define particles that dominate the distribution at ${\gamma }_{e}$ > 2.

Standard image High-resolution image

Further insight into the acceleration mechanism at ${\gamma }_{e}\gtrsim 6$ is provided by Figure 10. For a particle with Lorentz factor ${\gamma }_{e}(t)$, we calculated the "parallel" work done during the recent time interval ${\rm{\Delta }}t={t}_{\mathrm{IC}}({\gamma }_{e})$,

Equation (23)

where ${E}_{\parallel }={\boldsymbol{E}}\cdot {\boldsymbol{B}}/B$ and ${v}_{e}^{\parallel }={{\boldsymbol{v}}}_{e}\cdot {\boldsymbol{B}}/B$ are the electric field and the particle velocity components parallel to the local magnetic field ${\boldsymbol{B}}$. The condition of ideal magnetohydrodynamics ${E}_{\parallel }=0$ is maintained almost everywhere in the reconnection layer except the X-points, where a large ${E}_{\parallel }$ develops and nonideal effects become important. Thus, the work ${W}_{\parallel }(t)$ is a proxy for recent particle acceleration at X-points. As one can see from Figure 10, W is large for particles with ${\gamma }_{e}$ > 6 and becomes the main factor responsible for their high energies. We conclude that the high-energy end of the particle spectrum is the result of acceleration by E at X-points. The energy released through the X-point acceleration approximately equals the energy received (and radiated) by the population with ${\gamma }_{e}$ > 6; it accounts for ∼1% of the total energy released by magnetic reconnection.

Figure 10.

Figure 10. Diagnostics of particle acceleration by the nonideal electric field E in the radiative simulation (${\gamma }_{\mathrm{cr}}$ = 16). (a) The vertical axis represents the cumulative work W done on the particle by E during the recent time interval Δt equal to the IC cooling time tIC(γe). The color-coded 2D histogram was constructed and normalized similarly to that in Figure 9, with the green curve showing the median value. For particles with ${\gamma }_{e}$ ≫ 6, W accounts for a large fraction of their energy, approaching the relation ${W}_{\parallel }=({\gamma }_{e}-1){m}_{e}{c}^{2}$ shown by the dotted black line. (b) Contribution of particles with large W (blue curve) to the particle energy distribution in the reconnection layer (black curve). Here large W is defined as the region above the horizontal dashed line in panel (a). One can see that the nonideal acceleration dominates the distribution at ${\gamma }_{e}$ > 6.

Standard image High-resolution image

A more detailed analysis demonstrates two types of X-point acceleration: most particles with 6 ≲ ${\gamma }_{e}$ ≲ 10 are generated in the reconnection plane y ≈ 0, and particles with ${\gamma }_{e}$ > 10 are mainly generated in the secondary (vertical) reconnection layers formed at the interfaces between merging plasmoids. The higher energies achieved by particles in the secondary layers reflect the stronger magnetization: the magnetic field in plasmoids is stronger than in the inflow region, and so a merger-induced reconnection layer effectively has a higher σ than the nominal σ = 10.

In summary, we find that three distinct populations contribute to the spectrum in the reconnection region: (i) at moderate energies, ${\gamma }_{e}$ < 2, the particle distribution is dominated by bulk motions of large plasmoids cooled to nonrelativistic temperatures and pushed by magnetic stresses to mildly relativistic speeds; (ii) at intermediate energies, 2 ≲ ${\gamma }_{e}$ ≲ 6, the distribution is dominated by particles freshly picked up by the unstructured outflows from X-points toward neighboring plasmoids; and (iii) at the highest energies, ${\gamma }_{e}$ ≳ 6, the distribution becomes dominated by particles accelerated by E at X-points, either in the main reconnection layer or in the reconnection layers formed at the interface between merging plasmoids. The contributions of all these acceleration mechanisms are shown in Figure 11. The populations (ii) and (iii) account for ${f}_{\mathrm{HE}}$ ∼ 20% of the total power dissipated by reconnection.

Figure 11.

Figure 11. Energy distribution of particles in the radiative simulation (${\gamma }_{\mathrm{cr}}$ = 16), time-averaged over the interval 1.5 ≲ ct/L ≲ 5. The black solid curve shows the distribution in the reconnection region (same as in Figures 9 and 10), and the dotted black curve shows the distribution in the entire box, including the inflow region. The distribution in the reconnection region is further dissected into four components. This dissection is performed using particle spectra binned in vertical spatial slices of the reconnection region. Each slice has the thickness of 100 grid cells (i.e., 20 $c/{\omega }_{{\rm{p}}}$) along the x-direction, and it extends in the y-direction until the boundaries of the reconnection region (as identified by the mixing criterion described in Section 4.2). We also use the magnetic vector potential to identify the locations of X-points and the contours of plasmoids (see Sironi et al. 2016). If an X-point lies at the boundary between two neighboring plasmoids, it is identified as an X-point formed by plasmoid mergers, and the particle spectrum from that spatial slice contributes to the yellow curve. Otherwise, the X-point resides in the main reconnection layer, and the corresponding spectrum contributes to the green curve. Spatial slices that contain only particles residing inside plasmoids give the red curve. If none of these conditions are met, the slice intersects the unstructured outflow located in between plasmoids (and not containing X-points), which gives the blue curve.

Standard image High-resolution image

4.4. Dependence on the Radiation Field Intensity

In the previous sections we compared the results of two simulations: the fiducial model with strong cooling (${\gamma }_{\mathrm{cr}}$ = 16), and the nonradiative model (${\gamma }_{\mathrm{cr}}=\infty $). We have also run simulations with different values of ${\gamma }_{\mathrm{cr}}$ to study in more detail the effect of increasing radiative losses (decreasing ${\gamma }_{\mathrm{cr}}$ from $\infty $ down to 11.3). The results are presented in Figures 12 and 13.

Figure 12.

Figure 12. Cumulative distribution of plasmoid width w, for different values of ${\gamma }_{\mathrm{cr}}$ as indicated in the figure. The histogram (with Poissonian error bars) is normalized to the overall number of plasmoids Npl. The corresponding differential distribution is given by f(w) = dN(w)/dw. The predictions $N(w)\propto {w}^{-1}$ by Uzdensky et al. (2010) and Loureiro et al. (2012) and $N(w)\propto \mathrm{const}$ by Huang & Bhattacharjee (2012) are plotted as dashed and dotted black lines, respectively.

Standard image High-resolution image

Figure 12 shows the cumulative distribution of plasmoid sizes. It demonstrates that the statistics of plasmoid sizes is not appreciably affected by the level of cooling losses.

Figure 13 shows that the reconnection rate is nearly independent of the degree of IC losses. Furthermore, the plasma compression in the reconnection layer is only moderately increased by cooling: from $\langle {n}^{{\prime} }\rangle /{n}_{0}\sim 4$ in the nonradiative case to $\langle {n}^{{\prime} }\rangle /{n}_{0}\sim 7$ in the case of ${\gamma }_{\mathrm{cr}}$ = 11.3. This is in agreement with the results of Werner et al. (2019), but in contradiction with earlier works (Uzdensky & McKinney 2011; Uzdensky & Spitkovsky 2014; Bégué et al. 2017), which argued that cooling results in strong compression until the plasma pressure in the reconnection layer balances the upstream magnetic pressure. Instead, we find that the dynamic, magnetically dominated plasmoid chain sustains the required pressure with a negligible plasma contribution, close to the force-free regime, and thus avoids strong compressions. Most of the magnetic energy in the layer is contained in the plasmoids, and the resulting mean magnetic energy in the reconnection layer is nearly insensitive to IC losses.

Figure 13.

Figure 13. Dependence of the main properties of the reconnection system on radiative losses. The decreasing parameter ${\gamma }_{\mathrm{cr}}$ from left to right corresponds to increasing the level of radiative losses. Red and blue colors highlight the two simulations discussed in most detail in the text: ${\gamma }_{\mathrm{cr}}=\infty $ (nonradiative) and ${\gamma }_{\mathrm{cr}}$ = 16 (strongly radiative). (a) Reconnection rate, or inflow rate, in units of the speed of light; (b) average comoving particle density in the reconnection region (as defined in Section 4.2), in units of the upstream density n0; (c) comoving magnetic energy fraction ${\epsilon }_{B}^{{\prime} }$ (dashed) and internal energy fraction ${\epsilon }_{\mathrm{int}}^{{\prime} }$ (solid), averaged over the reconnection region (they are both normalized to the upstream rest-mass energy density); (d) rms bulk momentum $\langle {u}^{2}{\rangle }^{1/2}$ in the reconnection region (dashed) and maximum outflow bulk momentum umax (solid), with the horizontal dotted black line corresponding to the Alfvénic limit of $\sqrt{\sigma };$ (e) ratio of the power emitted by the high-energy particles (i.e., the component that cannot be accounted for by bulk motions; in Figure 8, it is shown as the green dashed line) and the net power emitted by all particles in the reconnection region. For the comoving quantities in panels (b) and (c), for each computational cell we transform particle and field quantities to the corresponding fluid rest frame, given the local bulk velocity ${\boldsymbol{\beta }}$. All points are obtained by time averaging over 1.5 ≲ ct/L ≲ 5.

Standard image High-resolution image

By contrast, the internal energy in the reconnection region drops by two orders of magnitude as ${\gamma }_{\mathrm{cr}}$ decreases from infinity to ${\gamma }_{\mathrm{cr}}$ = 11.3. The internal energy fraction ${\epsilon }_{\mathrm{int}}^{{\prime} }$ is defined as the average internal energy density in the reconnection region normalized to the rest-mass energy density of the upstream plasma. The internal energy density was evaluated in each computational cell by boosting the plasma energy density from the lab frame to the plasma rest frame (which has the local velocity ${\boldsymbol{\beta }}$), using the approximation of isotropic plasma stress tensor in the rest frame.14

The bulk motions in the reconnection layer somewhat slow down at low ${\gamma }_{\mathrm{cr}}$, because Compton drag opposes the accelerating magnetic forces. However, this effect is much weaker than the loss of internal energy, in agreement with analytical estimates in Section 4.2. In agreement with expectations, γβ stays below the $\sqrt{\sigma }$ limit, i.e., plasma motions in the layer are moderately sub-Alfvénic, even for the smallest and fastest plasmoids.

Figure 13 also shows the measured ${f}_{\mathrm{HE}}$, the ratio of the power emitted by the high-energy particles (i.e., the component that cannot be accounted for by bulk motions) and the net power emitted by all particles in the reconnection region. It drops from ${f}_{\mathrm{HE}}$ ∼ 90% at ${\gamma }_{\mathrm{cr}}=\infty $ to ${f}_{\mathrm{HE}}$ ≲ 30% at low ${\gamma }_{\mathrm{cr}}$, when strong radiative losses limit the lifetime of accelerated particles.

5. 3D Simulations

Real plasmas are not 2D, i.e., the invariance in the z-direction is not enforced. Then, the kink instability develops in the current sheet (Zenitani & Hoshino 2007; Sironi & Spitkovsky 2014). When the 3D reconnection process is viewed in the x-y cross section, it looks similar to that in 2D simulations—plasmoids are formed by the tearing instability and accelerated by magnetic stresses along the x-axis. The 2D plasmoids are the cross sections of magnetic flux ropes that extend in the z-direction. The kink instability bends the flux ropes, forcing them to entangle and reconnect.

The characteristic coherence length of a dynamic flux rope in the z-direction is comparable to its thickness (the plasmoid size in the x-y plane), and relativistic particles traverse this length on the light-crossing timescale. Thus, unlike the 2D simulations, the particles trapped in a flux rope are not completely buried. They get a chance to become repeatedly energized as they move along the bent flux rope, approach new reconnection sites, and interact with the fast outflows emanating from X-points. The repeated energization events may better counteract radiative losses.

To explore the effects of 3D dynamics on radiative reconnection, we have performed 3D simulations with a setup similar to that of our 2D simulations. Due to computational constraints, we choose a computational box with half-lengths in the z- and x-directions of Lz = Lx = L = 806 $c/{\omega }_{{\rm{p}}}$. We employ periodic boundary conditions in the z-direction. The skin depth is resolved with 2.5 cells, and we initialize the upstream plasma with one particle per computational cell. We keep σ = 10 and Bg/B0 = 0.1, the same as in our fiducial 2D runs. In order to (marginally) preserve the hierarchy of relevant scales and energies discussed in Section 3, we focus on the model with strongest cooling, which has ${\gamma }_{\mathrm{cr}}$ = 11.3 ≳ σ; then, the relation in Equation (11) still marginally holds.

The results are presented in Figures 14 and 15. Similar to our 2D runs (see panel (c) of Figure 1), we find that most of the reconnection region (in particular, the large plasmoids/flux ropes) is efficiently cooled down (Figure 14). The hottest plasma (and hence strongest IC emission) is localized in small regions where individual particles are quickly energized at X-points and in the fast outflows emanating from the X-points. This is clearly observed in the thin parts of the reconnection layer, where the freshly energized plasma has not cooled yet. The same happens at the interfaces between merging flux ropes where secondary current sheets are formed.

Figure 14. 3D structure of the reconnection layer at time t = 4.6(L/c) in our largest 3D run (Lz = L = 806 $c/{\omega }_{{\rm{p}}}$) with very strong IC losses (γcr = 11.3). We show the region $| y| /L\lt 0.45$ where reconnection occurs (the actual extent of the computational box along y grows with time as described in Section 2). We show the local average $\langle {\gamma }_{e}^{2}{\beta }_{e}^{2}\rangle =\langle {\gamma }_{e}^{2}\rangle -1$, obtained for each cell by averaging over the particles in the neighboring 5 × 5 × 5 cells. This quantity is proportional to the local total IC power radiated per particle. The figure shows four isosurfaces of $\langle {\gamma }_{e}^{2}{\beta }_{e}^{2}\rangle $, each colored according to the value of $\langle {\gamma }_{e}^{2}{\beta }_{e}^{2}\rangle $. One can see that large plasmoids are cold (blue and cyan), whereas the thin regions of the reconnection layer, where the freshly energized plasma has not cooled yet, are sites of bright IC emission (green and orange). An animation of this figure is available. The video shows the time evolution t = 0 to t = 5.6 L/c. The real-time duration of the video is 8 s. The animation is also available at http://user.astro.columbia.edu/~lsironi/Site/share/ic_in_reconnection/.

(An animation of this figure is available.)

Video Standard image High-resolution image

In Figure 15 we compare the particle energy distributions found in the 3D and 2D simulations with the same magnetization and cooling rate. One can see that two 3D simulations with different box sizes Lz gave nearly identical results, indicating convergence, and these results are similar to those of the 2D simulation. Particle acceleration is somewhat more efficient in the 3D simulation in the range 3 ≲ γe ≲ 10; however, this effect is modest.15 Regardless of the dimensionality of the computational domain, 2D or 3D, the peak of the energy distribution (γe ≲ 2) is controlled by the bulk motions, i.e., the particles populating the peak are cold.

Figure 15.

Figure 15. Particle distribution ${dN}/d\mathrm{log}({\gamma }_{e}-1)$ in the 3D and 2D simulations with ${\gamma }_{\mathrm{cr}}$ = 11.3, time-averaged in the interval 1.5 ≲ ct/L ≲ 5.8. For each simulation, the distribution was calculated in the reconnection region (solid curves) and in the entire computational box (dotted curves). The dashed curves show the distribution in the reconnection region when accounting only for bulk kinetic energy (i.e., we plot ${dN}/d\mathrm{log}(\gamma -1)$ as a function of γ − 1). The color-coding represents different box half-lengths along the z-direction: Lz = 0 in cyan (i.e., a 2D simulation), Lz = 115 $c/{\omega }_{{\rm{p}}}$ in green, and Lz = 806 $c/{\omega }_{{\rm{p}}}$ in red. The box half-length along the x-direction is L = 806 $c/{\omega }_{{\rm{p}}}$ in all the models. The distributions are normalized to the number of particles in the reconnection region.

Standard image High-resolution image

6. Conclusions

Kinetic plasma simulations provide an efficient tool to study the nonlinear dynamics of magnetic reconnection. This paper focused on the radiative regime, where heated particles radiate their energies on a timescale much shorter than the light-crossing time of the reconnection layer. Our simulation setup was motivated by the hard X-ray activity of accreting black holes, which can be powered by radiative reconnection in magnetically dominated regions of the corona. In this paper, we studied radiative reconnection in a pair-dominated plasma with a high magnetization parameter σ ≫ 1, which is likely found near the black hole spin axis.

For simplicity, we fixed the radiation field during the simulation and assumed that it is composed of low-energy photons that scatter in the Thomson regime (see also Werner et al. 2019). As the main parameter controlling radiative (Compton) cooling we chose ${\gamma }_{\mathrm{cr}}$ (defined in Section 2.1), which is proportional to the cooling timescale and inversely proportional to the radiation density. We have studied reconnection with various levels of radiative losses by performing simulations with different ${\gamma }_{\mathrm{cr}}$, from ${\gamma }_{\mathrm{cr}}=\infty $ (nonradiative) down to ${\gamma }_{\mathrm{cr}}$ = 11.3 (extremely radiative). Running the simulations in 2D allowed us to use large computational boxes, with length 2L up to 13,440 $c/{\omega }_{{\rm{p}}}$. This permits good separation of spatial and energy scales, which is particularly important in the radiative regime. We have also performed a 3D simulation for the case of ${\gamma }_{\mathrm{cr}}$ = 11.3. Its results suggest that our main conclusions from the 2D simulations will hold in 3D models.

Our conclusions are as follows:

  • 1.  
    In radiative magnetic reconnection, plasma particles are energized as they enter the reconnection layer, and then quickly cool and become locked in plasmoids at a low temperature kT ≪ 100 keV. This makes thermal Comptonization unable to generate 100 keV X-rays.
  • 2.  
    The plasma energy in the reconnection layer (and its inverse Compton emission) is dominated by mildly relativistic bulk motions of cold plasmoids, which are pulled against Compton drag by magnetic stresses. This confirms the radiative mechanism described in B17. We find that, for a magnetically dominated pair plasma, the plasmoid bulk motions have an effective temperature of ∼100 keV and their inverse Compton emission mimics thermal Comptonization.
  • 3.  
    Particle acceleration in radiative reconnection is less efficient compared with the nonradiative regime. However, there are accelerated particles in the reconnection layer, which form a high-energy tail of the particle distribution function. We have identified two mechanisms injecting energetic particles: acceleration at X-points and particle pickup by outflows from X-points. Both mechanisms are nearly impulsive (i.e., operate on a timescale much shorter than the cooling timescale). Slower acceleration processes, in particular Fermi reflection and magnetic compression (e.g., Petropoulou & Sironi 2018; Guo et al. 2019), are suppressed in the radiative regime. We found that about 20% of the dissipated reconnection power goes into high-energy particles. Their inverse Compton emission will generate a high-energy tail of the radiation spectrum, which may explain the MeV tail detected in the hard state of Cyg X-1 (McConnell et al. 2002).
  • 4.  
    Radiative losses influence the picture of plasma dynamics in the reconnection layer: radiation efficiently removes thermal plasma pressure, exerts Compton drag on bulk motions, and clears cavities inside moving plasmoids (Figure 1). However, these dynamic effects weakly change the net reconnection rate and magnetic field structure in the reconnection layer. Our results support the conjecture of B17 that the magnetic forces moving plasmoids are approximately the same in the radiative and nonradiative regimes. Furthermore, the statistics of plasmoid sizes are similar (Figure 12). We conclude that the plasmoid chain formation is largely controlled by magnetic stresses, and plasma pressure plays a minor role even in the hot (nonradiative) regime. We also note that similar plasmoid chains were observed in "force-free" simulations that completely neglected plasma pressure and inertia (Parfrey et al. 2015), again illustrating the dominant role of magnetic stresses in shaping the structure of relativistic reconnection.

In summary, our first-principle simulations support radiative magnetic reconnection in the magnetically dominated regime as the mechanism powering the hard state of accreting black holes. Our results confirm the radiative mechanism proposed in B17, give the effective temperature of bulk motions in the reconnection layer kTb ∼ 100 keV, and give the fraction of dissipated reconnection power deposited into high-energy particles ${f}_{\mathrm{HE}}$ ∼ 20% (much lower than in previous nonradiative simulations).

In all our models reconnection occurred in e± plasma with magnetization σ = 10. Recent studies of nonradiative reconnection in electron–proton and electron–positron–proton plasmas have shown that for σ ≳ 1 leptons still pick up a significant fraction of the dissipated magnetic energy (Rowan et al. 2017, 2019; Werner et al. 2018; Petropoulou et al. 2019). Future work will explore how the outcome of radiative reconnection depends on σ and ion abundance.

It may also be useful to investigate how the results change in the presence of a strong guide field, stronger than 0.1B0 assumed in our simulations. We also leave for future work detailed calculations of the X-ray spectrum. Such calculations can be done similarly to the Monte Carlo simulations in B17, but with the particle energy distribution provided by first-principle kinetic plasma simulations.

L.S. acknowledges support from DoE DE-SC0016542, NSF ACI-1657507, NASA ATP NNX17AG21G, and NSF PHY-1903412. A.M.B. is supported by NASA grant NNX17AK37G, Simons Foundation grant No. 446228, and the Humboldt Foundation. The simulations have been performed at Columbia (Habanero and Terremoto) and with NASA-HEC (Pleiades) and NERSC (Cori) resources.

Footnotes

  • For the same reason, we prescribe by hand that the particles belonging to the hot population initialized in the layer are not subject to IC losses.

  • If reconnection were to transfer all of the field energy to the particles, the mean particle energy would be $\sim \sigma {m}_{e}{c}^{2}/2$. So, our definition of ${r}_{0,\mathrm{hot}}$ corresponds, apart from a factor of two, to the Larmor radius of the typical particles heated/accelerated by reconnection.

  • The box size along y increases over time, and at the end it is comparable to or larger than the x extent.

  • In the following, by "dissipated reconnection power" we indicate the energy per unit time transferred by reconnection from the fields to the particles. This is ∼50% of the Poynting flux carried toward the reconnection layer, since half of the incoming energy stays locked in the magnetic field (Sironi et al. 2015).

  • Faster bulk motions in the y-direction are attained for the radiative case ${\gamma }_{\mathrm{cr}}$ = 16 (left panel) than for the uncooled case ${\gamma }_{\mathrm{cr}}=\infty $ (right panel). This is related to the density cavities at the front of cooled plasmoids, where the local magnetization reaches high values (due to the very low density), resulting in fast reconnection outflows in between merging plasmoids.

  • Hereafter primed quantities refer to the plasmoid comoving frame. The transverse size w is invariant under Lorentz boosts along the direction of plasmoid motion.

  • 10 

    This substitution approximately corresponds to measuring the plasmoid speeds at the exits of the computational box (x = ±L). It somewhat overestimates the average ${t}_{\mathrm{age}}^{{\prime} }({w}_{\max })$, as many plasmoids end their lives in mergers, at a smaller ${t}_{\mathrm{age}}^{{\prime} }$. Therefore, Equation (19) somewhat overestimates $u({w}_{\max })$.

  • 11 

    This high-energy slope is significantly steeper than $p\sim 2$ usually reported by PIC simulations of relativistic nonradiative reconnection (e.g., Guo et al. 2014; Sironi & Spitkovsky 2014; Werner et al. 2016). As discussed in Sironi et al. (2016), this difference comes from the choice of boundary conditions at $y=\pm L$, where the plasmoids exit the box. Our simulation employed the outflow boundary conditions, so the growing plasmoids disappear as they exit the box. This leads to a relatively large contribution from plasma between the plasmoids. In contrast, in simulations with periodic boundaries, the plasmoids stay in the box and keep growing. Then, the particle distribution becomes dominated by the largest plasmoids, filled with a broad flat distribution. Sironi et al. (2016) demonstrated that, regardless of the chosen boundary conditions, the particle spectrum inside large plasmoids has the slope $p\sim 2$, and with open boundary conditions the interplasmoid bump around ${\gamma }_{e}\sim 0.5\sigma $ emerges above the plasmoid spectrum.

  • 12 

    The high-energy tail extends up to γe ∼ 30, contrary to the anticipated cutoff at γe ∼ ${\gamma }_{\mathrm{cr}}$ = 16 (see Equation (3)). The cutoff is pushed to higher γe because the reconnection layers between merging plasmoids have a stronger B > B0 and hence a stronger E ≈ 0.1B, capable of pushing particles to higher energies.

  • 13 

    In the simple case with Bg = 0, the electromagnetic fields in the unstructured outflow at y = 0 are dominated by Ez (the reconnection electric field) and By (the reconnected magnetic field), and the drift speed ${\boldsymbol{E}}\times {\boldsymbol{B}}/{B}^{2}$ can reach $| {E}_{z}/{B}_{y}| \sim {v}_{A}/c\approx 1$. The picked-up particle initially follows a cycloid in the x-z plane (starting along z and then getting deflected along x), with a relativistic drift in the x-direction $u\sim \sqrt{\sigma }$ and a comparable 4-velocity of gyration in the drift frame.

  • 14 

    See Sironi et al. (2016) for further details. In a few cases, we have checked that a Lorentz transformation of the full stress energy tensor to the plasma rest frame yields nearly identical results (i.e., that our values are not sensitive to the assumption of isotropic stress tensor).

  • 15 

    At higher energies, γe ≳ 10, the difference between 2D and 3D results may be attributed to the difference in merger dynamics of 2D plasmoids versus 3D flux ropes. However, particles with γe ≳ 10 make a negligible contribution to the overall particle census and the IC power.

Please wait… references are loading.
10.3847/1538-4357/aba622