Radiative Reconnection-powered TeV Flares from the Black Hole Magnetosphere in M87

To quantify the effects of cooling and its dynamical importance with respect to the acceleration by reconnection, it is helpful to define a dimensionless parameter, γがんま_rad, as ${E}_{\mathrm{rec}}| e| \equiv 2{\sigma }_{T}U{\gamma }_{\mathrm{rad}}^{2}$, where the left-hand side is the accelerating force experienced by particles in the reconnection sites, while the right-hand side corresponds to a radiation drag force either due to synchrotron or IC (where we assume the Thomson regime and σしぐま_T is the Thomson cross section) radiation. Here, E_rec is the electric field strength in the so-called x-point regions of the reconnecting current sheet, where the bulk of particle energization and energy dissipation takes place. Its value is controlled by the strength of the reconnecting magnetic field and the rate at which the magnetic energy is dissipated (reconnection rate); it can be written as E_rec ≈ βべーた_rec B. We take βべーた_rec = 0.1, which is characteristic for collisionless relativistic reconnection. The value of U in the definition of γがんま_rad can either be interpreted as the energy density of the magnetic field, U = U_B ≡ B²/8πぱい, or as the energy density of the background (seed) photon field, depending on whether we are considering synchrotron or IC cooling. When the forces are balanced, particles with energies γがんま ≫ γがんま_rad will cool down faster (either via synchrotron or IC) than they are able to accelerate, whereas cooling will be negligible on acceleration timescales for particles with γがんま ≪ γがんま_rad. Thus, comparing γがんま_rad with the magnetization that powers the acceleration, σしぐま, offers an important insight for understanding the cooling efficiency in the context of acceleration by reconnection.

For synchrotron radiation, which for M87 is by far the strongest cooling mechanism, the parameter can be estimated as ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}={(4\pi {\beta }_{\mathrm{rec}}| e| /({\sigma }_{T}B))}^{1/2}\sim 4\cdot {10}^{6}{\left(B/100\,{\rm{G}}\right)}^{-1/2}$. One can already draw two important conclusions from these estimates. First, since $\sigma \gtrsim {\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$, the reconnection proceeds in the radiative regime where most of the dissipated magnetic energy is quickly converted into synchrotron radiation on timescales much shorter than the dynamical time of the system. This suggests that the energy density of the synchrotron photons in the steady state is comparable to the dissipated magnetic energy density, ${U}_{\mathrm{rad}}^{\mathrm{syn}}\approx {\beta }_{\mathrm{rec}}{U}_{B}$. Note that the synchrotron cooling is negligible in x points, where the magnetic field vanishes, meaning that the fast acceleration process is uninhibited by synchrotron radiation; the highest-energy particles will still accelerate to ${E}_{{e}^{\pm }}\approx \sigma {m}_{e}{c}^{2}$. The peak of the synchrotron emission will then be set by the burnoff limit: ${\varepsilon }_{c}\approx {\hslash }{\omega }_{B}{\left({\gamma }_{\mathrm{rad}}^{\mathrm{syn}}\right)}^{2}\sim 20\,\mathrm{MeV}$ (Uzdensky & Spitkovsky 2014). ¹⁰

The optical depth for synchrotron photons produced in reconnection is small τたう_{γがんま γがんま} ≪ 1 for M87 (Ripperda et al. 2022), yielding a very nonlocal pair production in the whole region of size w in the black hole magnetosphere. These pairs, produced by photon–photon collisions, have arbitrary pitch angles to the local magnetic field and thus radiate rapidly, giving rise to a distinct synchrotron spectral component different from the synchrotron continuum from reconnection. As the main fuel for these secondary pairs are the synchrotron photons from reconnection, their characteristic Lorentz factors at the time of birth can be estimated as ${\gamma }_{\sec }\approx (\mathrm{few}\cdot {\varepsilon }_{c})/({m}_{e}{c}^{2})\sim 200$. The synchrotron cooling timescale for the secondary pairs, ${t}_{\sec }$, compared to the light-crossing time, is ${t}_{\sec }c/w\approx {({\gamma }_{\sec }{l}_{B})}^{-1}\sim {10}^{-3}$, where l_B = σしぐま_T U_B w/(m_e c²) ∼ 1 is the magnetic compactness parameter (Beloborodov 2017), which means that a significant fraction of the energy contained in these pairs is radiated away inside the magnetosphere. The synchrotron peak energy associated with these pairs is ${\varepsilon }_{\sec }\approx {\hslash }{\omega }_{B}{\gamma }_{\sec }^{2}\sim 0.01$–0.1 eV, with the energy density of resulting photons being controlled by the optical depth ${U}_{\mathrm{rad}}^{\sec }\approx {U}_{\mathrm{rad}}^{\mathrm{syn}}{\tau }_{\gamma \gamma }\sim 0.006$ erg cm⁻³ and the luminosity being close to ${L}_{\mathrm{rad}}^{\sec }\approx {L}_{\mathrm{rec}}{\tau }_{\gamma \gamma }\sim {10}^{38}$ erg s⁻¹.

Pairs accelerated by reconnection to energies ${E}_{{e}^{\pm }}\gtrsim {10}^{6}$–10⁷ m_e c² can Compton-upscatter lower-energy photons to TeV energies (the upper limit will be determined by the maximum pair energy). There will be several variations of this TeV signal. First of all, the reconnection region is filled with a soft background photon field from the accretion disk with average energies around εいぷしろん_soft ≈ 0.001 eV (300 GHz) and a characteristic energy density ${U}_{\mathrm{rad}}^{\mathrm{soft}}\approx 0.01\,\mathrm{erg}\,{\mathrm{cm}}^{-3}$ (Broderick & Tchekhovskoy 2015; EHT MWL Science Working Group et al. 2021). On the other hand, the region is also filled with synchrotron photons both from the pairs accelerated in reconnection and from the secondary pairs produced in the magnetosphere (as discussed in the previous section). Primary synchrotron photons from reconnection possess a wide range of energies up to a few εいぷしろん_c ≈ 10–100 MeV. The energy density of these photons is comparable to the magnetic field energy density ${U}_{\mathrm{rad}}^{\mathrm{syn}}\approx {\beta }_{\mathrm{rec}}{U}_{B}\sim 40\,\mathrm{erg}\,{\mathrm{cm}}^{-3}$, whereas the energy density of the secondary synchrotron photons is the aforementioned ${U}_{\mathrm{rad}}^{\sec }\sim 0.006$ erg cm⁻³.

A TeV signal can be produced by upscattering either of the photon components populating the reconnection region by the highest-energy pairs with γがんま ≳ 10⁶: the soft background (we will refer to this channel as IC), the primary synchrotron photons (synchrotron self-Compton, SSC), and the secondary synchrotron photons (secondary synchrotron self-Compton, SSC2). In all of these cases, most of the power will be focused in the narrow frequency range corresponding to the highest-energy pairs, i.e., above a few hundred GeV (Blumenthal & Gould 1970).

We will assume that the reconnection produces a power-law distribution of pairs with a cutoff, ${f}_{{e}^{\pm }}(\gamma )\propto {\gamma }^{-p}{e}^{-\gamma /{\gamma }_{c}}$, where the cutoff is γがんま_c ≈ σしぐま, and 1 < p < 1.5 (A. Chernoglazov et al. 2023, in preparation), an assumption that we later verify using kinetic simulations. Further in this section, we will be comparing the Compton luminosities peaking in TeV (IC, SSC, and SSC2) with the synchrotron luminosity, which peaks at few tens of MeV. The power density of this emission is equal to ${P}_{\mathrm{syn}}\approx 2{\sigma }_{T}{{cn}}_{{e}^{\pm }}\langle {\gamma }^{2}{\tilde{B}}_{\perp }^{2}/(8\pi )\rangle$, where the averaging is done over the particle distribution $\langle x\rangle \equiv \int d\gamma {f}_{{e}^{\pm }}(\gamma )x$ and ${\tilde{B}}_{\perp }$ is the effective magnetic field strength perpendicular to the motion of the particle (see Equation (4) for the full definition). As demonstrated in our simulations in Section 3, we can accurately approximate $\langle {\gamma }^{2}{\tilde{B}}_{\perp }^{2}/(8\pi )\rangle \approx {\chi }^{2}{\left({\gamma }_{\mathrm{rad}}^{\mathrm{syn}}\right)}^{2}{U}_{B}$ (as long as the cooling is strong), with U_B being the magnetic energy density in the unreconnected upstream and χかい ≈ 1 being a dimensionless parameter of order unity. Using that, one finds ${P}_{\mathrm{syn}}\,\approx 2{\sigma }_{T}{{cn}}_{{e}^{\pm }}{\chi }^{2}{\left({\gamma }_{\mathrm{rad}}^{\mathrm{syn}}\right)}^{2}{U}_{B}$.

2.3.1. Total Synchrotron Luminosity

In the previous section, we assumed that the radiated synchrotron power during the reconnection is a nonnegligible fraction of the jet power (L_syn ≈ L_rec ≈ f_rec L_jet, where f_rec ∼ 0.1). To obtain this estimate more rigorously, let us multiply the emitted synchrotron power density by the volume to get L_syn = P_syn αあるふぁ πぱい w² δでるた_cs/(2πぱい), where αあるふぁ ∼ 1 is the azimuthal extent of the region undergoing reconnection, and w and δでるた_cs are the radial extent and the thickness of the current sheet. As we demonstrate from the simulations in Section 3, the characteristic temperature (or the mean energy) of particles in the current sheet is dictated by the synchrotron burnoff limit: $3{T}_{{e}^{\pm }}\,=\left(\kappa {\gamma }_{\mathrm{rad}}^{\mathrm{syn}}\right){m}_{e}{c}^{2}$, where κかっぱ is a dimensionless parameter defined as follows: $\langle \gamma \rangle =\kappa {\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$. Given that, we can eliminate the number density, ${n}_{{e}^{\pm }}$, using the pressure balance condition: ${n}_{{e}^{\pm }}\approx {U}_{B}/{T}_{{e}^{\pm }}$. The width of the current sheet in collisionless reconnection is set by the characteristic Larmor radius of particles (e.g., Uzdensky & Spitkovsky 2014): δでるた_cs ≈〈γがんま〉m_e c²/(∣e∣B). The total radiated synchrotron power then reads:

$$\begin{eqnarray}&&{L}_{\mathrm{syn}}\approx \displaystyle \frac{9\alpha }{2\pi }{\beta }_{\mathrm{rec}}\mathop{\underbrace{\displaystyle \frac{\langle {\gamma }^{2}{\tilde{B}}_{\perp }^{2}\rangle }{{\left({\gamma }_{\mathrm{rad}}^{\mathrm{syn}}\right)}^{2}{B}_{0}^{2}}}}\limits_{{\chi }^{2}}{L}_{\mathrm{BZ}}\approx \mathop{\underbrace{1.4\cdot \alpha {\chi }^{2}{\beta }_{\mathrm{rec}}}}\limits_{\sim 0.1}{L}_{\mathrm{BZ}},\end{eqnarray} \tag{ 3 }$$

where ${L}_{\mathrm{BZ}}\equiv {B}_{0}^{2}{r}_{g}^{2}c/24$ is the Blandford & Znajek (1977) luminosity for a maximally spinning black hole, and we further assumed that the magnetic field falls with the distance from the horizon: $B={B}_{0}{(r/{r}_{g})}^{-1}$. Notice that the dependence on w has disappeared, due to the fact that the magnetic field strength decreases linearly with distance and so does the effective value of ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$. Provided that the jet power is close to the BZ value L_jet ≈ L_BZ, we can state that radiatively efficient reconnection has a potential of emitting about 10% of the jet power, i.e., indeed, f_rec ∼ 0.1.

As discussed earlier in this section, the bulk of the radiated power, L_syn, will be carried by the photons with energies close to the burnoff limit, εいぷしろん_c ∼ 20 MeV. The very-high-energy (TeV) signal in our model is produced via Compton upscattering of lower-energy photons by the accelerated pairs and will thus contain only a small fraction of the total synchrotron power. To estimate the power in the TeV component of the observed emission, in Sections 2.3.2, 2.3.3, and 2.3.4 we investigate three different channels, which essentially consist of three populations of photons that can be upscattered to produce the TeV emission.

2.3.2. Upscattering of Soft Photons from the Disk

2.3.3. Upscattering of Synchrotron Photons from Secondary Pairs

2.3.4. Upscattering of Reconnection-produced Synchrotron Photons

Pair production may greatly inhibit the observed luminosity of both the synchrotron MeV and the inverse-Compton TeV signal if the optical depth for these photons is large. As we estimated above, the optical depth for MeV photons interacting with each other is small, τたう_{γがんま γがんま} ∼ 10⁻⁴. Outgoing TeV photons, however, will interact with ∼1 eV photons from the soft background. We may estimate the optical depth to that process similarly: ${\tau }_{\mathrm{TeV}}\sim 0.1{\sigma }_{T}{{wU}}_{\mathrm{rad}}^{\mathrm{eV}}/(1\,\mathrm{eV})$, where ${U}_{\mathrm{rad}}^{\mathrm{eV}}$ is the energy density of background photons with energies around 1 eV. Given the total energy density of the soft radiation, ${U}_{\mathrm{rad}}^{\mathrm{soft}}$, as well as the slope of the distribution beyond εいぷしろん_soft ≈ 10⁻³ eV, $\varepsilon {dn}/d\varepsilon \propto {\varepsilon }^{-1.2}$ (Broderick & Tchekhovskoy 2015), we can take ${U}_{\mathrm{rad}}^{\mathrm{eV}}\approx {U}_{\mathrm{rad}}^{\mathrm{soft}}{\left({\varepsilon }_{\mathrm{soft}}/1\,\mathrm{eV}\right)}^{0.2}$. We then find an upper limit for the optical depth: τたう_TeV ≈ 1, meaning that a large fraction of the TeV signal escapes the inner few r_g .

Our predictions are shown in Figure 2, where we plot the emerging spectra from all the emission components of the reconnecting sheet. In this plot, we fix the main (synchrotron) luminosity to 10⁴³ erg s⁻¹, while other components scale appropriately. The synchrotron component from the secondary pairs (magenta line) is roughly estimated by taking a power-law distribution of pairs with a cutoff around ${\gamma }_{\sec }$ (see Section 2.2). We vary a few parameters of the problem: the effective magnetic field strength, B, in the region where most of the synchrotron emission takes place, ¹³ the power-law index of the accelerated pairs, p, and their cutoff energy, γがんま_c . For comparison, we overplot the spectrum of soft photons from the disk at energies ≲10¹⁶ Hz from Broderick & Tchekhovskoy (2015).

Zoom In Zoom Out Reset image size

The main takeaway of these estimates is that the outgoing TeV signal is comprised of the IC scattering of soft background photons from the disk, synchrotron photons produced in situ by the reconnecting sheet, and from the cooling of secondary pairs outside the reconnecting sheet. Our estimates show that the Compton scattering of soft background photons (IC component) contributes most to the very-high-energy signal, with the total TeV luminosity up to ∼0.01%–0.1% of the jet power, or about 10³⁹–10⁴¹ erg s⁻¹. Synchrotron radiation from secondary pairs appears as a distinct component at energies ∼0.01–0.1 eV, with a luminosity up to 10³⁸–10³⁹ erg s⁻¹.

All of our simulations are initialized with a cold (T_e ≪ m_e c²) uniform background pair plasma of density n₀ and a hot, thin, overdense layer of width Δでるた along the x-coordinate in the middle of the box (at y = 0). The magnetic field is initialized as ${\boldsymbol{B}}={B}_{0}\hat{{\boldsymbol{x}}}\tanh \left(y/{\rm{\Delta }}\right)$ (where B₀ is the magnetic field strength in the far upstream), and E = 0 everywhere. The width, the temperature, and the density of the layer, as well as the in- and out-of-plane velocities of the particles in the layer, are chosen in such a way that the setup is initially in the Harris equilibrium (Harris 1962). Near the x boundaries of the box, we impose absorbing boundary conditions that damp the fields and absorb outgoing particles (Cerutti et al. 2016). Along the x-axis at a fixed distance from the layer, we inject leptons to replenish the upstream plasma. These boundary conditions allow us to run our simulations for many light-crossing times of the box, ensuring the results are insensitive to the initial conditions.

In all of our simulations, the value of the magnetization in the upstream is fixed: $\sigma ={B}_{0}^{2}/\left(4\pi {n}_{0}{m}_{e}{c}^{2}\right)={10}^{3}$. The size of the box is fixed at a value of 160 × 80 in units of σしぐま ρろー₀, where ${\rho }_{0}={m}_{e}{c}^{2}/\left(| e| {B}_{0}\right)$ is the Larmor radius of fiducial particles with a velocity γがんま βべーた = 1. The value σしぐま ρろー₀ can thus be interpreted as the characteristic Larmor radius of particles with γがんま βべーた ∼ σしぐま in the field equal to the background value; the size of the box is chosen to be large enough to contain many gyro-orbits of the most energetic particles. Our simulations have a resolution of d₀ = c/ωおめが_p0 = 5Δでるたx, with d₀ and ωおめが_p0 being the skin depth and the plasma frequency of the background plasma; this means that the size of the box is 24,000 × 13,000 cells. The particle distribution function is sampled with five particles per cell (PPC) in the upstream, with significantly more in the reconnection layer; the results are virtually unchanged for higher PPC. Deposited currents are smoothed with eight passes of a digital filter with a stencil (1/4, 1/2, 1/4) at each time step.

To model synchrotron cooling, our simulations also include the radiative drag force in the Landau & Lifshitz (1975) formulation:

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{f}}}_{\mathrm{rad}} & = & \displaystyle \frac{{\sigma }_{T}}{4\pi }\left({\boldsymbol{\kappa }}-{\gamma }^{2}{\tilde{B}}_{\perp }^{2}{\boldsymbol{\beta }}\right),\,\,\,\mathrm{where}\\ {\boldsymbol{\kappa }} & = & \left({\boldsymbol{E}}+{\boldsymbol{\beta }}\times {\boldsymbol{B}}\right)\times {\boldsymbol{B}}+\left({\boldsymbol{\beta }}\cdot {\boldsymbol{E}}\right){\boldsymbol{E}},\\ {\tilde{B}}_{\perp }^{2} & = & {\left({\boldsymbol{E}}+{\boldsymbol{\beta }}\times {\boldsymbol{B}}\right)}^{2}-{\left({\boldsymbol{\beta }}\cdot {\boldsymbol{E}}\right)}^{2},\end{array}\end{eqnarray} \tag{ 4 }$$

assuming a particle moves with a four-velocity γがんま βべーた in the background electromagnetic fields, E and B . To be able to properly model the effect of the drag force at realistic timescales, we upscale its magnitude introducing the dimensionless parameter, ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$, discussed in Section 2.1. The magnitude of the force is scaled in such a way that ∣ f _rad∣ = 0.1B₀∣e∣ for $\gamma ={\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$, ∣ B ∣ = B₀, ∣ E ∣ = 0, and βべーた ⊥ B (with B₀ being the strength of the upstream magnetic field). The ratio of ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}/\sigma$ determines the regime of the synchrotron cooling, with ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}/\sigma \lesssim 1$ corresponding to the dynamically strong cooling regime.

To start the reconnection process, we remove the pressure balance in the middle of the box, which triggers the formation of two magnetic islands propagating away from the center of the box. After about one light-crossing time, fast plasmoid-mediated magnetic reconnection is triggered across the whole current layer.

In Figure 3 we plot three different quantities for three simulations with varying cooling strength (${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}/\sigma$). Panels (a1)–(a3) show the overall pair-plasma density in units normalized to the upstream density. In these panels, one can clearly see the main features of the reconnecting current sheet: the plasmoids—circular magnetic structures containing the accelerated plasma, and the x points in between them, where the magnetic field vanishes and the main energy dissipation takes place. Plasmoids, which travel along the current sheet, adiabatically grow and merge with each other, are held intact via the balance between the magnetic and plasma pressure. From panels (a1)–(a3), one can clearly see that when the synchrotron cooling "removes" some of that plasma pressure, densities (and sizes) of the plasmoids have to adjust to maintain the proper balance.

Zoom In Zoom Out Reset image size

In panels (b1)–(b3), we plot the mean energy of particles in each simulation cell, normalized to m_e c². From these three panels, it is evident that the dimensionless parameter ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$ controls the relativistic temperature inside the plasmoids. As the cooling strength increases, i.e., the ratio ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}/\sigma$ drops, the temperature inside the plasmoids decreases, while the density increases, maintaining the pressure (the product of the two) roughly constant. We see no evidence of the variation in the plasmoid magnetic field strength for different cooling regimes.

In panels (c1)–(c3) of Figure 3 we show the total synchrotron emissivity (synchrotron power radiated from a unit volume) for the highest-energy photons. This quantity is defined as $\int d\gamma {f}_{e}^{\pm }{\gamma }^{2}{\tilde{B}}_{\perp }^{2}$, where ${f}_{e}^{\pm }$ is the distribution function of high-energy particles (with γがんま > σしぐま/4) and ${\tilde{B}}_{\perp }$ is the effective magnetic field component perpendicular to the direction of the particle velocity, defined in Equation (4). When there is no synchrotron cooling, or when the cooling is dynamically weak (i.e., the regime least applicable to the reconnection layer in the magnetosphere of M87*), plasmoids carry most of the high-energy particles. These particles lose their energies at timescales much longer than the characteristic injection timescale from the current sheet, and the plasmoids get rapidly replenished with fresh accelerated plasma. In this regime, most of the high-energy photons are thus produced inside the largest plasmoids ((c1)–(c2)). In the strong cooling case, on the other hand, the highest-energy particles are rapidly cooled after they leave the acceleration sites (x points). Because of that, in panel (c3) we can only observe the smallest (youngest) plasmoids, as well as the relativistic outflows in the sheet carrying the freshly accelerated pairs, as sources of energetic synchrotron photons.

In the radiatively efficient regime, ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}/\sigma \lesssim 1$, a fraction which we will denote χかい², of the dissipated magnetic energy inside the current layer, p_diss ≈ ∣e∣E_rec c, is radiated via the synchrotron mechanism: ${p}_{\mathrm{syn}}\approx 2{\sigma }_{T}c\langle {\gamma }^{2}{\tilde{B}}_{\perp }^{2}/(8\pi )\rangle ={\chi }^{2}{p}_{\mathrm{diss}}$ (both powers are computed per particle). Here E_rec ≈ βべーた_rec B₀ is the accelerating electric field in the layer. Using the definition of ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$ one can find $\langle {\gamma }^{2}{\tilde{B}}_{\perp }^{2}\rangle \approx {\chi }^{2}{\left({\gamma }_{\mathrm{rad}}^{\mathrm{syn}}\right)}^{2}{B}_{0}^{2}$, a formula that was used in Section 2.3.1 when computing the total radiated synchrotron power. To find the value of χかい , we compare four different simulations with varying cooling strengths: one with weak cooling (${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}/\sigma =2.5$), one with moderate cooling (${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}/\sigma =1$), and the other two with strong cooling (${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}/\sigma =0.2,\,0.5$). The values for the dimensionless parameter χかい measured from simulations are shown in Table 1; the results indicate that for a wide range of values of ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$, our approximate analytic argument of χかい being constant holds rather well (we inspect this more closely in Appendix C). ¹⁴ We also cite the values of $\kappa \equiv \langle \gamma \rangle /{\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$, which were used to estimate the effective temperature of the current sheet. And again, varying the cooling strength within an order of magnitude only marginally affects this parameter, which is close to 0.1–0.3.

Table 1. Characteristic Values of Physical Parameters That Have Observational Significance Measured Directly from Our Simulations with Different Cooling Strengths

${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}/\sigma$	βべーたrec	χかい	κかっぱ	p	γがんまc
∞	0.25	⋯	⋯	1.1	2830
2.5	0.28	0.24	0.10	1.5	1440
1.0	0.28	0.30	0.17	1.5	1100
0.5	0.28	0.33	0.22	1.6, 2.4	1340
0.2	0.28	0.33	0.29	1.6, 2.4	1100

Note. The first column ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}/\sigma$ indicates the cooling strength, with the infinity corresponding to the run without radiation; βべーた_rec is the time-averaged reconnection rate (normalized to c); the ratio $\chi \equiv \langle {\gamma }^{2}{\tilde{B}}_{\perp }^{2}{\rangle }^{1/2}/\left({\gamma }_{\mathrm{rad}}^{\mathrm{syn}}{B}_{0}\right)$ provides an estimate for the total radiated synchrotron power; $\kappa \equiv \langle \gamma \rangle /{\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$ correlates the effective temperature with the cooling burnoff limit; p and γがんま_c are the best-fit power-law index and the cutoff energy for particle distribution functions ${f}_{e}^{\pm }\propto {\gamma }^{-p}{e}^{-\gamma /{\gamma }_{c}}$. In panels where we cite two values of p, we fit a broken power law with an ankle near $\gamma \approx {\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$.

Download table as:  ASCIITypeset image

We also directly measure the reconnection rate by evaluating the component of the E × B /∣ B ∣² velocity upstream, pointing toward the current sheet (in Figure 3 that corresponds to the ±y direction). We find that the strength of the synchrotron cooling is unimportant in determining the reconnection rate, with only marginal variations likely caused by the intermittency. The value of the rate itself is close to βべーた_rec ∼ 0.25–0.28, in agreement with what has been found by Guo et al. (2015), and Sironi et al. (2016) in 2D simulations.

PIC simulations also enable us to directly measure both the energy distribution of pairs and the spectra of both synchrotron and IC photons from first principles (shown in Figure 4).Particle distribution functions for different values of the cooling strength, ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}/\sigma$, are shown in panel (a) with different colors. Deduced power-law indices and cutoff energies for the best-fit power law (or double power law) of the form ${f}_{{e}^{\pm }}\propto {\gamma }^{-p}{e}^{-\gamma /{\gamma }_{c}}$ (or $\propto {\gamma }^{-{p}_{1}}$ for $\gamma \lt {\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$, and $\propto {\gamma }^{-{p}_{2}}{e}^{-\gamma /{\gamma }_{c}}$ for $\gamma \gtrsim {\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$) are shown in Table 1.

Zoom In Zoom Out Reset image size

In the weak-to-marginal cooling regime, when ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}/\sigma \gtrsim 1$, the distribution of pairs is best described by a single power law, with an index varying between 1 and 1.5 (which justifies the values used in our analytical model in the previous section). When the cooling becomes stronger, we see the second power law generated beyond $\gamma \gtrsim {\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$ with a slightly steeper index of ≈2.5. In all of the cases, however, except for the uncooled regime, the energy cutoff is strictly controlled by the value of σしぐま, since the highest-energy pairs are produced in the x points where there is no synchrotron cooling. In the uncooled case, which is unrealistic for our astrophysical scenario, the spectrum extends much further than σしぐま due to various secondary acceleration channels (Petropoulou & Sironi 2018; Hakobyan et al. 2021; Zhang et al. 2021).

In Figure 4(b) we show both the synchrotron spectra computed on the fly, as well as the IC spectra with an assumed soft isotropic background radiation computed in post-processing. The peaks of the synchrotron emission in all of the cases correspond to ${\hslash }{\omega }_{B}{\left({\gamma }_{\mathrm{rad}}^{\mathrm{syn}}\right)}^{2}$, which agrees with the earlier predictions by Uzdensky & Spitkovsky (2014). The peaks for the IC signal are self-consistently dictated by the energies of leptons. There are several important differences in the analytic expectation (shown with a dashed line). First, at lower energies, our simulation underpredicts the synchrotron flux. This is likely caused by a combination of two reasons: (a) limited scale separation leading to a smaller number of low-energy particles in the sheet (we only consider radiation from particles in the current layer), and (b) in 2D, plasmoids, where most of the radiation takes place, are compressed, and magnetic field strengths are stronger than the upstream, which can yield an on average marginally higher radiation frequency. Also notice that at higher frequencies, $\varepsilon \gtrsim {\hslash }{\omega }_{B}{\left({\gamma }_{\mathrm{rad}}^{\mathrm{syn}}\right)}^{2}$, in the synchrotron spectrum, we see slightly elevated luminosity, especially for the strongest cooling cases. This is due to the fact that in the strong cooling regime, particles beyond $\gamma \gtrsim {\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$ cool very rapidly. This makes the spectrum more intermittent, and, when averaged over time, results in more flux at high energies (see Hakobyan et al. 2019). To see this, let us estimate the characteristic maximum photon energy or the cutoff energy. We can use the peak synchrotron energy formula: ${E}_{\mathrm{cut}}\,\,\approx (3/2){\hslash }{\gamma }_{\max }^{2}(| e| {\tilde{B}}_{\perp })/({m}_{e}{c}^{2})$, where ${\gamma }_{\max }\approx \sigma$, and ${\tilde{B}}_{\perp }$ is the perpendicular component of the magnetic field for the highest-energy particles. Since $\sigma \gtrsim {\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$ (strong cooling regime), particles with Lorentz factors close to ${\gamma }_{\max }$ rapidly lose their energy via synchrotron radiation, emitting the perpendicular component of their momenta. Because of this, the effective ${\tilde{B}}_{\perp }$ is smaller than the value of the upstream field and may be approximated as ${\gamma }_{\max }{\tilde{B}}_{\perp }\approx {\gamma }_{\mathrm{rad}}^{\mathrm{syn}}{B}_{0}$ (see also Appendix C). We may then simplify the expression for the spectral cutoff: ${E}_{\mathrm{cut}}\approx (3/2){\hslash }{\omega }_{B}{\left({\gamma }_{\mathrm{rad}}^{\mathrm{syn}}\right)}^{2}\left(\sigma /{\gamma }_{\mathrm{rad}}^{\mathrm{syn}}\right)\approx 16\,\mathrm{MeV}\left(\sigma /{\gamma }_{\mathrm{rad}}^{\mathrm{syn}}\right)$. We can then conclude that the spectral cutoff for the synchrotron emission is at higher energies for the stronger cooling regimes.

In our simulations, we had to downscale some of the physical parameters to be able to tackle the problem computationally. In particular, the expected magnetization parameter of the upstream plasma, σしぐま, during a flare from the magnetosphere of M87* is expected to be close to 10⁷ (see Equation (2)), whereas in our simulations, we employed σしぐま = 10³. At the same time, the synchrotron limit, ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$, was downscaled from about 10⁶ to 200. In light of this, there might be a concern that some of the arguments we have made in this section might not extrapolate when the problem is rescaled to realistic parameters. First of all, it is important to notice that we still retain the strong-cooling hierarchy, ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}\lt \sigma$, which conserves the balance between the competing mechanisms (acceleration and cooling) when the dimensionless values are downscaled. On top of that, employing σしぐま ∼ 10³ allows us to have at least a decade of scale separation between σしぐま, ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$, and ${u}_{{ \mathcal A }}/c\sim \sqrt{\sigma }$ (Alfvénic four-velocity). We also do not expect the overall dynamics of the layer to change drastically when ${u}_{{ \mathcal A }}/c\gg 1$, as for these values, the Alfvénic three-velocity, which dictates the characteristic velocities of the flow in the sheet, approaches the speed of light.

Our simulations are 2D, while the real current sheet in the vicinity of the black hole is inherently 3D. In particular, 3D reconnection has been demonstrated (see Zhang et al. 2021) to enable a fast acceleration channel for particles beyond γがんま ∼ σしぐま, which in 2D operates at a much slower rate (Petropoulou & Sironi 2018; Hakobyan et al. 2021). Ironically, the peculiarity of radiative reconnection helps us in this scenario, since in this regime any secondary acceleration channel outside the x points is essentially forbidden due to the strong cooling that vanishes only in these microscopic regions comparable in size to the plasma skin depth.

In this paper, we have shown that large-scale magnetic reconnection, occurring in the episodically forming current sheets in magnetically arrested accretion flows around black holes, can power bright multiwavelength flares extending to TeV energies. We have shown that the reconnecting sheet can dissipate a significant fraction of the jet power on timescales controlled by the universal reconnection rate: βべーた_rec ≈ 0.1–0.3; see Equation (3). While the efficient dissipation of the Poynting flux during magnetic reconnection is not a novel concept, ¹⁵ plasma and magnetic field parameters that control the dissipation determine the observational footprint of the process in a wide range of wavelengths.

Our findings indicate that synchrotron radiation is the main cooling mechanism for the leptons, through which the dissipated magnetic energy is converted into radiation. The synchrotron photons from the accelerated plasma, which peak at around a few tens of MeV, have finite optical depth and are thus susceptible to abundant pair production. This process in turn loads the current sheet with electron–positron pairs; the expected resulting pair multiplicity with respect to the Goldreich–Julian density is given by Equation (1), and its estimated value is rather large: ∼10⁷.

We have shown that the expected magnetization parameter upstream of the reconnecting current sheet, which determines the available magnetic energy per lepton, is large: σしぐま ≈ 10⁷; see Equation (2). Despite strong synchrotron cooling, reconnection can accelerate pairs to energies close to a few to tens of TeV, ¹⁶ or ${E}_{{e}^{\pm }}\approx \sigma {m}_{e}{c}^{2}$. The highest-energy pairs can then Compton-scatter both the soft photons from the disk and the synchrotron photons from the sheet, producing the very-high-energy TeV component of the emission. Our estimations suggest that the upscattering of the soft thermal photons produced in the bulk of the accretion flow is the most efficient way of producing a TeV signal, with a predicted peak luminosity of 0.01%–0.1% of the jet power, or about 10³⁹–10⁴¹ erg s⁻¹.

We have also verified some of our assumptions using first-principles 2D PIC simulations of reconnecting current sheets with synchrotron cooling included self-consistently. In particular, our simulations support the claim that the high-energy cutoff of the distribution of particles is not sensitive to the synchrotron cooling strength (see Table 1), since the particle acceleration sites (x points) have zero magnetic field strength and are thus not susceptible to synchrotron cooling. On top of that, both the power-law index (which we found to be close to p ≈ 1.5) and the reconnection rate (βべーた_rec ≈ 0.28) are also insensitive to the synchrotron cooling, as long as the cooling is strong.

Using particle distributions from our simulations, we have generated both synchrotron and IC spectra (shown in Figure 4(b)), with particle Lorentz factors rescaled to realistic values. For the weak cooling regime, our generated spectra overall match analytic predictions; both the spectral peak and the cutoff of the synchrotron component are determined by the value of ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$ and are thus close to ${E}_{\mathrm{peak}}\approx {\hslash }{\omega }_{B}{\left({\gamma }_{\mathrm{rad}}^{\mathrm{syn}}\right)}^{2}\approx 20\,\mathrm{MeV}$. However, in the strong cooling regime, while the peak energy is still set by the value of ${\gamma }_{\mathrm{rad}}^{\mathrm{syn}}$, the spectral cutoff is higher ${E}_{\mathrm{cut}}\approx {E}_{\mathrm{peak}}(\sigma /{\gamma }_{\mathrm{rad}}^{\mathrm{syn}})$. This intermittent synchrotron excess at higher energies leads to a much less eminent gap between the synchrotron component (between a few to ten MeV and GeV) and the IC component (above 100 GeV), with the overall spectrum being more continuous.

GRMHD simulations indicate that during the formation of the current sheet, the accretion flow in the magnetospheric region where reconnection takes place is repelled by the magnetic stresses beyond roughly w ∼ 10r_g . As a consequence, the radio to near-IR flux from the disk is expected to dim significantly during the flare (H. Jia & et al. 2023, in preparation). In Section 2.2 we discussed the possible emergence of a millimeter-to-infrared synchrotron signal from the secondary pairs produced in the reconnection upstream (shown in pink in Figure 2). The luminosity of this signal is unlikely to exceed 0.1%–0.01% of the primary synchrotron counterpart (or about 10³⁸–10³⁹ erg s⁻¹), as the pair-production optical depth, τたう_{γがんま γがんま} ∼ 10⁻⁴, is rather small. Distinguishing it from the dimmed background emission of the accretion flow will thus be challenging; however, a possible feature one should be looking for is enhanced values of the polarization degree in near-IR during the flare. The synchrotron radiation from secondary pairs is produced in the reconnection upstream (i.e., at the jet base), where the magnetic field is coherent, as the turbulent accretion flow would have been repelled. ¹⁷ To summarize, we expect that during the TeV flare, the emission of the disk below eV should dim, as the disk is repelled and the reconnection region is evacuated, while the polarization degree in the same frequency range can increase, owing to the contribution from the synchrotron emission of secondary pairs.

Let us also briefly discuss the temporal characteristics of the flare at different wavelengths. The dimensionless reconnection rate, βべーた_rec, controls both the timescale of the eruption during which the magnetic energy dissipation takes place, as well as its peak synchrotron energy; εいぷしろん_c ≈ (3βべーた_rec/αあるふぁ_F )m_e c²; see Equation (3). In the relativistic collisionless regime, studied using PIC simulations for almost a decade, the value of the reconnection rate varies between 0.1 and 0.3 depending on the magnetization parameter, the dimensionality of the simulation, and the boundary conditions (our 2D simulations find values close to 0.25–0.3, as shown in Table 1; 3D simulations tend to have smaller values; Guo et al. 2015). This prediction is an order of magnitude faster than what has been found by relativistic MHD simulations (Ripperda et al. 2019). For the largest current sheets, 3D GRMHD simulations show a flux decay period, i.e., a reconnection event that can power a flare, of ∼100r_g/c (Ripperda et al. 2022), which is approximately 30 days for M87*. Due to an order of magnitude faster reconnection rate, we expect the flux decay period, and hence the characteristic flare timescale, in collisionless plasma to be shorter (Bransgrove et al. 2021), which aligns well with observational flare durations of the order of a few days for M87*. ¹⁸

The magnetic compactness of the inner region near the horizon is significant: l_B = σしぐま_T U_B w/(m_e c²) ≈ few (Beloborodov 2017), with w ≈ 10r_g being the length of the current sheet. This corresponds to the cooling timescale of t_cool ≈ (w/c)/(γがんま l_B ) for particles with a Lorentz factor γがんま. Since the cooling time of even the marginally relativistic particles is short, the variability on timescales comparable to ∼r_g /c can only be expected from global effects, such as the beaming due to the orientation of the current sheet. At lower energies, the flare can also have shorter intermittency caused by the dynamics of relativistic outflows, the growth and collisions of plasmoids, and particle acceleration and cooling. For example, in the radio range, which is expected to be produced by the particles with Lorentz factors of the order of γがんま ∼ 100, the variability timescale can be comparable to 0.1r_g /c (a few hours for M87*). At higher energies, the same estimate yields ${t}_{\mathrm{int}}\approx {t}_{\mathrm{cool}}{({\hslash }{\omega }_{B}{\gamma }^{2}={\varepsilon }_{s})\approx 10\,\min \cdot ({\varepsilon }_{s}/\mathrm{eV})}^{-1/2}$. While this short-timescale intermittency will be virtually undetectable at the peak energies of MeV (around a second), its detection is possible in the optical to X-ray band, where we expect the characteristic cooling timescale to be around a few tens of seconds to minutes.