Black Hole Flares: Ejection of Accreted Magnetic Flux through 3D Plasmoid-mediated Reconnection

Bright flaring from accreting black holes is seen at all wavelengths, but the mechanism powering high-energy flares is still a topic of major debate. Rapid γ-ray flares have been observed from active galactic nuclei, in the form of very-high-energy (>100 GeV) emission (Aharonian et al. 2007; Albert et al. 2007; Aharonian et al. 2009; Aleksić et al. 2014). The variability timescale of the flares can be shorter than the light-crossing time of the event horizon, constraining the emitting region to be of the order of a Schwarzschild radius. Bright TeV flares are also periodically observed from the supermassive black hole M87*, in the center of the Messier 87 galaxy (Aharonian et al. 2006; Acciari et al. 2010; Aliu et al. 2012; Blanch 2021). The flares show a flux rise and decay timescale of 1–3 days, emitting ≳10⁴¹ erg s⁻¹ (Abramowski et al. 2012), which is nonnegligible compared to the total jet power of 10⁴²–10⁴⁴ erg s⁻¹ (e.g., Prieto et al. 2016). High-energy γ-rays originating nearby the horizon can be absorbed by background photons to create electron–positron pairs, preventing their escape. Therefore, it is unclear if there is a mechanism that can produce such flares near the horizon and under which conditions the radiation can freely escape. Furthermore, the black hole in the Galactic Center, Sgr A*, shows intriguing infrared and X-ray flares on similarly short dynamical timescales (Baganoff et al. 2001; Eckart et al. 2004; Neilsen et al. 2015) originating from near the horizon (Gravity Collaboration et al. 2018; Collaboration et al. 2021).

Magnetically arrested disk (MAD; Bisnovatyi-Kogan & Ruzmaikin 1974, 1976; Narayan et al. 2003) accretion is the most plausible scenario for the accretion flow onto active galactic nuclei showing strong jets (see, e.g., Event Horizon Telescope Collaboration et al. 2021 for M87*). Sources fed by stellar winds, like Sgr A*, are also capable of producing MADs (Ressler et al. 2020). General-relativistic magnetohydrodynamics (GRMHD) simulations show that a large amount of poloidal (pointing in the R- and z-directions) magnetic flux (proportional to the square root of the mass accretion rate) is forced into the black hole by the accreting gas, until the flux becomes dynamically important and strong enough to push the accreting gas away (Igumenshchev et al. 2003; Igumenshchev 2008; Tchekhovskoy et al. 2011). The MAD state is accompanied by large-amplitude fluctuations, caused by quasiperiodic accumulation and escape of the magnetic flux bundles in the vicinity of the black hole (Igumenshchev 2008; Tchekhovskoy et al. 2011; Dexter et al. 2020; Porth et al. 2021).

Recently, extreme-resolution two-dimensional (2D) GRMHD simulations showed that escape of magnetic flux bundles from the black hole, resulting in the decay of magnetic flux on the horizon, occurs through plasmoid-mediated reconnection (Ripperda et al. 2020, hereafter RBP20). The magnetic flux decay is accompanied by the ejection of the accretion disk (Proga & Begelman 2003). The ejection results in the formation of a magnetosphere, consisting of an equatorial plasmoid-unstable current sheet of an oppositely directed magnetic field that separates two highly magnetized jet regions. Reconnection in the current sheet releases energy that can power a flare, and the tension of the reconnected flux can push gas away and suppress the mass accretion rate. The jets, which supply matter in the current sheet, are highly magnetized because their large-scale magnetic field serves as a barrier to ions within the accretion disk. Pair discharges can generate ample electron–positron plasma to fill the magnetospheric region (e.g., Crinquand et al. 2020). The collisional mean free path of particles is much larger than the characteristic length scale of the system. As a result, the magnetospheric electron–positron plasma is collisionless and can be accelerated in a reconnecting current sheet into a power-law distribution and subsequently power high-energy flares. In magnetized and collisionless plasma conditions, reconnection occurs in the plasmoid-mediated regime at a universal reconnection rate of v_rec/v_A ∼ 0.1, where v_rec is the inflow velocity into a current sheet, and v_A ∼ c is the Alfvén speed (Guo et al. 2014; Sironi & Spitkovsky 2014; Werner et al. 2015).

In collisional systems, as described by GRMHD, the reconnection rate in the plasmoid-mediated regime at high Lundquist numbers (and at sufficiently high resolution to resolve the spatial scales associated to that Lundquist number) converges to a universal value of v_rec/v_A ∼ 0.01, becoming independent of the resistivity (Bhattacharjee et al. 2009; Uzdensky et al. 2010; Ripperda et al. 2019a; RBP20). ¹² Resolving plasmoid-mediated reconnection, and hence a converged universal reconnection rate, in global black hole simulations requires resolutions higher than ∼2000 cells in the θ-direction to capture thin current sheets liable to the plasmoid instability (RBP20; Bransgrove et al. 2021). The flare timescale is governed by the flux decay, which is directly set by the reconnection rate (Bransgrove et al. 2021); this makes it particularly important to resolve the plasmoid instability in thin current sheets.

Our goal here is to understand if a macroscopic reconnecting current sheet can form and power a flare in 3D GRMHD simulations with a converged universal reconnection rate, v_rec ∼ 0.01c, for the largest current sheets in the system, despite the excitation of nonaxisymmetric effects, like a Rayleigh-Taylor-type instability (RTI) preventing the complete arrest of accretion (Tchekhovskoy et al. 2011; Papadopoulos & Contopoulos 2019). In this letter, we conduct the highest-resolution global 3D GRMHD simulations to date to show that plasmoid-mediated magnetic reconnection in transient, nonaxisymmetric current sheets can power flares from accreting black holes and that the magnetic flux decay on the black hole event horizon is governed by the universal reconnection rate.

Throughout the manuscript, we use geometrized units with gravitational constant, black hole mass, and speed of light G = M = c = 1, such that length scales are normalized to the gravitational radius r_g = GM/c² and times are given in units of r_g c⁻¹. We employ spherical Kerr–Schild coordinates, where r is the radial coordinate, θ and ϕ are the poloidal and toroidal angular coordinates, respectively, and t is the temporal coordinate.

Reconnecting current sheets are plasmoid-unstable for Lundquist numbers (Bhattacharjee et al. 2009; Uzdensky et al. 2010)

$$\begin{eqnarray}&&S={v}_{{\rm{A}}}w/{\eta }_{\mathrm{num}}\geqslant {S}_{\mathrm{crit}}={10}^{4},\end{eqnarray} \tag{ 1 }$$

assuming the Alfvén speed v_A ∼ c, and the length of a current sheet w ∼ r_g. Here, we assume that the numerical resistivity proportional to the cell size is η_num ∝ Δx^p , where p ≈ 2 depending on the details of the second-order accurate algorithm. Thus, the constraint on S (Equation (1)) directly determines the required resolution. In the plasmoid-mediated regime, the reconnection rate converges to the asymptotic v_rec ∼ 0.01c in GRMHD (RBP20; Bransgrove et al. 2021), directly determining the (converged) rate of magnetic flux decay on the horizon. To achieve the resolution required to capture the plasmoid-mediated reconnection and, hence, achieve long-sought convergence in the reconnection rate, we employ our GPU-accelerated GRMHD code H-AMR (Liska et al. 2019). We set the effective numerical resolution to N_r × N_θ × N_ϕ = 5376 × 2304 × 2304 (dubbed "extreme resolution" from here onward) to ensure that we capture thin plasmoid-unstable current sheets (RBP20). To study convergence of the reconnection rate and the rate at which magnetic flux can escape from the black hole, we also conduct three lower resolution runs at N_r × N_θ × N_ϕ = 2240 ×1056 × 1024 ("high resolution"); 580 × 288 × 256 ("standard resolution"); and 288 × 128 × 128 ("low resolution"). The resolution in the r − and θ − dimensions is satisfied throughout the domain. To keep the cell aspects ratio approximately uniform in our spherical grid, we use three internal and four external derefinement levels (Liska et al. 2019) in ϕ to reduce the resolution from the full N_ϕ = 128–2304 at 30° <θ < 150° to N_ϕ = 16–18 within 05–75 of each pole. In all of these runs, we fix the radial domain to [1.2, 2000]r_g , and we use a minimum 10,000r_g c⁻¹ integration time. We use outflow boundary conditions in r, transmissive boundary conditions in θ, and periodic boundary conditions in ϕ, as described in Liska et al. 2018. We initialize our simulation to obtain a prograde MAD around a Kerr black hole with dimensionless spin a = 0.9375, starting from a torus threaded by a single weak poloidal magnetic field loop, defined by the vector potential ${A}_{\phi }\,\propto \max \left[\rho /{\rho }_{\max }{\left(r/{r}_{\mathrm{in}}\right)}^{3}{\sin }^{3}\theta \exp \left(-r/400\right)-0.2,0\right]$, normalized to the gas-to-magnetic-pressure ratio β = 2p/b² = 100. We replenish the gas density ρ in low-density regions to maintain ${\sigma }_{\max }=25$ where the magnetization σ = b²/(4π ρ c²) is defined using the magnetic field strength b in the frame comoving with the fluid, and fluid-frame rest-mass density ρ. We adopt an equation of state for a relativistic ideal gas with an adiabatic index of $\hat{\gamma }=13/9$, in between a fully relativistic gas $\hat{\gamma }=4/3$ and a fully nonrelativistic gas $\hat{\gamma }=5/3$. We employ dimensionless temperature units T = p/ρ with thermal gas pressure p, where T = 1 corresponds to k_B T = m_i c² with ion mass m_i and Boltzmann's constant k_B, such that T > 1 indicates relativistic ion temperatures.

We analyze the flaring mechanism and its properties in the MAD after t ≈ 5000r_g c⁻¹ when the accretion flow has settled into a quasi–steady state of a constant mass accretion rate and magnetic flux on the black hole event horizon (see Figure 8 in Appendix C). The accumulation of magnetic flux on the horizon cannot continue beyond the limit in which the outward magnetic force balances the inward gravitational force. When the magnetic flux reaches this limit in axisymmetry (2D), accretion is halted completely and a low-density magnetosphere with an equatorial current sheet can form transiently (RBP20). In 3D, a large spectrum of RTI modes develops in the turbulent inner edge of the disk, steadily driving accretion. The magnetic flux periodically erupts from the black hole into the disk. These eruptions are made possible by near-event-horizon reconnection, which converts the magnetic energy into the energy of emitting particles and can naturally power a flare. Figures 1 (at ϕ = 0, i.e., the meridional plane) and 2 (at θ = π/2, i.e., the equatorial plane) show the gas temperature T = p/ρ with magnetic field lines plotted as green lines, the gas-to-magnetic-pressure ratio β = 8π p/B², and rest-mass density ρ around the time of one such flare at t ∼ 9500r_g c⁻¹. Namely, we show the quantities in the quiescent period (i.e., a period of quasi-constant magnetic flux at the horizon) before, during, and after the large magnetic flux eruption, respectively, at t = 9122r_g c⁻¹, t = 9422r_g c⁻¹, and t = 9782r_g c⁻¹ (where we zoom out to show large-scale effects). Shortly before and during a flare, accretion only occurs through large-scale (i.e., low azimuthal mode number) spiral RTI modes (see also Takasao et al. 2019 for a very similar scenario explaining protostellar flares) creating a transient, nonaxisymmetric (i.e., over an angle ϕ < 2π), magnetized (i.e., low plasma-β), low-density magnetosphere (top and middle rows in Figures 1 and 2) pushing the accretion disk outward and resulting in a drop in the mass accretion rate. A macroscopic equatorial current sheet forms in the magnetosphere, extending from the horizon to the disk at $x=r\sin \theta \cos \phi \approx -5{r}_{{\rm{g}}}$ at $z=r\cos \theta \approx 0$ shown by the antiparallel magnetic field lines (inset in panel (D); green lines). Reconnection pinches off the horizontal magnetic field in the sheet, transforming it into vertical (z) magnetic field, reminiscent of the 2D results of RBP20. The flux eruption originates from the inner magnetosphere where the highly magnetized plasma in the jet directly feeds the current sheet. The plasma density in the jet is determined by the density floor at ${\sigma }_{\max }=25$ in our simulations, whereas in reality it is much more strongly magnetized ($\sigma \gg {\sigma }_{\max }$) pair plasma. Reconnection occurs locally in X-points where a field line breaks and reconnects to other field line (see insets in Figure 1(D) and 1(E)). In these X-points, reconnection heats the plasma up to $T\sim {\sigma }_{\max }=25$ (left panels) after which it is expelled from the layer at Lorentz factors up to ${\rm{\Gamma }}\propto \sqrt{{\sigma }_{\max }}=5$ (Lyubarsky 2005; see also Appendix B for an exploration of different ${\sigma }_{\max }$ in 2D). The flux is expelled through reconnection into the low-density region in between the large low-mode-number RTI modes accreting spirals. Electrons and positrons accelerated to nonthermal energies through reconnection at the X-points in the macroscopic equatorial current sheet can power high-energy flares that may reach a distant observer during the drop in the mass accretion rate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Small plasmoids are visible close to the horizon, and a larger hot plasmoid is detected at x = −3r_g (middle row in Figure 1) as a result of the merger of smaller escaping plasmoids. The plasmoids that escape the gravitational pull of the black hole interact with the disk and jet sheath resulting in significant heating up to at least z ≳ ±40r_g. The bottom row of Figure 1 shows a large magnetic flux tube at x ≈ 20–30r_g: a low-density region with strong vertical field (low plasma-β) heated to medium temperature T ∼ 0.1 − 1. The flux tube forms as a result of the reconnection that converts the horizontal magnetic field into a vertical field that is ejected from the reconnection layer. Filled with heated plasma, the flux tube can appear as a hot spot. The accumulated vertical magnetic flux in this hot spot can remain coherent for approximately one orbital timescale between 10 and 30r_g (bottom row in Figure 2 between y ≈ − 20r_g and y ≈ 20r_g), while the inner 10r_g is already in the quiescent accretion state at t = 9782r_g c⁻¹. RTIs develop at the boundary of the hot spot, which mixes the hot low-density plasma into the surrounding accreting gas. The hot spots are expected to be filled with positrons and electrons energized by the reconnection, which in this way can end up in the accretion disk. After the flaring episode, magnetic flux builds up on the horizon and the quasi-steady-state accretion cycle develops again. Smaller and less hot current sheets where B^ϕ changes sign also exist in the inner ∼20r_g of the turbulent accretion disk during the quiescent period, indicated by thin high-β layers of antiparallel field lines (top and bottom rows in Figures 1 and 2).

Figure 3(A) visualizes the 3D nature of the hot current sheet by showing the temperature and magnetic field line structure in the inner 10r_g during the flare at t = 9422r_g c⁻¹. The current sheet has a relativistic temperature T > 1, whereas shortly before the flare at t = 9122r_g c⁻¹ (3(B)) there are no structures at T > 1. During the flare, the (green) field lines in the current sheet (i.e., seeded in the T > 1 region in 3(A)) have a clear spiral structure and are separated from the more vertical field lines in the disk (blue). During the quiescence before the flare (Figure 3(B)), no such distinction is visible, and all field lines (green and blue, which are seeded at the same points as in panel 3(A)) are part of the disk. The extreme resolution allows to capture multiple plasmoids identified as 3D helical field line structures in the sheet (Figure 3(C)) during the magnetic flux eruption. We highlight a typical X-point as the manifestation of reconnection, separating an infalling (purple field line) and escaping flux tube (green field line) in the hot current sheet. Similar X-points can be detected in e.g., the inset in Figure 1(D).

Zoom In Zoom Out Reset image size

Figure 4(A)–(D) zooms into the current sheet during large magnetic flux eruptions for the four numerical resolutions employed. The drop in magnetic flux at low and standard resolutions (panels (A), (B)) is not accompanied by a large drop in the mass accretion rate (see panels (E), (G)), due to the large numerical diffusion. The magnetic field diffuses through the thick current sheet and does not reconnect, due to the large numerical resistivity. This results in a too high reconnection rate and a large heated area (see Appendix C and Figure 9 for more properties of the large magnetic flux eruption at low resolution). The current sheet is in these cases not plasmoid-unstable. The high-resolution flux eruption (panel (C)) behaves similarly to the extreme-resolution result (panel (D)) from Figure 1, indicating that the plasmoid instability is resolved on the grid, and that the reconnection rate is converged to a universal value of 0.01c (panel (H)). In Figure 4, we also analyze the magnetic flux ${\dot{\phi }}_{\mathrm{BH}}:= \tfrac{1}{2}{\int }_{0}^{2\pi }{\int }_{0}^{\pi }{| }^{* }{F}^{{rt}}| \sqrt{-g}\,d\theta d\phi$ on the horizon (Figure 4(E)) and the mass accretion rate $\dot{m}:= -{\int }_{0}^{2\pi }{\int }_{0}^{\pi }\rho {u}^{r}\sqrt{-g}d\theta d\phi$ through the inner 5r_g (Figure 4(G)), where g is the metric determinant, u^μ is the fluid four-velocity, ^* F^{μ ν} is the dual of the Faraday tensor, and ρ is the fluid-frame rest-mass density. After ∼5000r_g c⁻¹, the flow sets into a quasi–steady state, which is globally converged for all resolutions. For the extreme-resolution run (magenta line in Figure 4(E)), we observe two major flux decays, which we associate with large flares, at t ≈ 7300r_g c⁻¹ and t ≈ 9300r_g c⁻¹, both lasting for a few ∼100r_g c⁻¹. We also observe a small flux decay at t ≈ 6800r_g c⁻¹, associated with a smaller flare, or "miniflare". For all flares, the magnetic flux on the event horizon decays quasi-exponentially with time with characteristic timescale τ ≈ 500 r_g c⁻¹ (indicated by the black dashed lines), implying that the decay is governed by reconnection at a universal rate of 0.01c, consistent with the decay observed for a split monopole magnetic field on the event horizon (Bransgrove et al. 2021). All three events are accompanied by a large drop in the mass accretion rate (Figure 4(G)) that is related to the ejection of the accretion disk such that the accretion is funneled through a small azimuthal angle ϕ < 2π and nearly halts.

Zoom In Zoom Out Reset image size

For the high-resolution run (red line in Figure 4(E)), similar flare episodes can be observed at t ≈ 7500r_g c⁻¹, t ≈ 8300r_g c⁻¹, and t ≈ 9400r_g c⁻¹, with flux decaying on the same timescale τ ≈ 500 r_g c⁻¹. For lower resolutions (blue and green lines), there is a clearer distinction: large flares show (e.g., at t ≈ 7300r_g c⁻¹ for low resolution, and t ≈ 8300r_g c⁻¹ and 8600r_g c⁻¹ at standard resolution) a faster decay rate τ ≈ 350 r_g c⁻¹, implying a faster reconnection rate > 0.01c. Miniflares (e.g., at t ≈ 9300r_g c⁻¹ for low resolution and t ≈ 7500r_g c⁻¹ for high resolution) instead show a flux decay at a rate of τ ≈ 500 r_g c⁻¹ implying a reconnection rate of ∼0.01c. At low and standard resolutions, these miniflares are typically not accompanied by a clear drop in ${\dot{m}}_{5{r}_{{\rm{g}}}}$ (Figure 4(G)), while large flares are showing a clear drop in ${\dot{m}}_{5{r}_{{\rm{g}}}}$ implying a large (≳5r_g) current sheet. This can be explained by the large numerical diffusion of the thinning current sheet in both the z- and y-directions, resulting in a too broad accretion funnel at low and standard resolution (Figures 4(A), (B)). Miniflares are better captured at lower resolutions than large flares due to the shorter length of the current sheet and the higher effective resolution of the spherical grid at small radii (see Appendix D).

The reconnection rate can be determined directly by selecting a current sheet during a flare episode and measure the inflow speed of the plasma into the reconnection layer. To do so, we first transform the Eulerian velocity and magnetic fields into a locally Minkowski frame (see, e.g., White et al. 2016) to apply standard reconnection analysis in flat spacetime (RBP20). The fields are expressed in minimum variance coordinates (Howes 2016), with B^L projected in the flat frame along the poloidal direction parallel to the current sheet, B^M along the toroidal direction, and B^N perpendicular to the current sheet, to determine the upstream geometry, showing a typical Harris-type sheet structure in Figure 4(F). Both the toroidal and poloidal components switch sign in the sheet, indicating that zero-guide-field reconnection occurs. The total vertical velocity of the flow consists of the inflow of the fluid into the current sheet due to reconnection, v_rec, and the advection of the current sheet with the bulk velocity, v_bulk. In Figure 4(H), we then measure the flow speeds from left and right of the current sheet as v_in,left = (v_bulk + v_rec)/(1 + v_bulk v_rec/c²) and v_in,right = (v_bulk − v_rec)/(1 − v_bulk v_rec/c²) and solve for v_rec, where we account for the relativistic speed of the bulk flow. We determine the profile of the upstream magnetic field projected along the current sheet and find the location where the profile becomes flat (Figure 4(F)). We then select 10 cuts of the current sheet at different radii and consistently find a reconnection rate of ∼0.01c, indicating a Lundquist number of at least $S={v}_{{\rm{A}}}w/{\eta }_{\mathrm{num}}={({v}_{\mathrm{rec}}/c)}^{-2}={10}^{4}$. Reconnection thus occurs in the asymptotic plasmoid-mediated regime where S ≥ S_crit = 10⁴ (Bhattacharjee et al. 2009) for our extreme-resolution run, where the length of the sheet w ≳ r_g, Alfvén speed in the upstream, ${v}_{{\rm{A}}}=\sqrt{{\sigma }_{\mathrm{up}}/({\sigma }_{\mathrm{up}}+1)}c\sim c$, for σ_up = 25, and numerical resistivity η_num. The reconnection rate is consistent with 2D resistive GRMHD simulations of plasmoid-dominated reconnection in MAD flows (RBP20) for Lundquist numbers $S={L}_{\mathrm{sheet}}c/\eta \gtrsim { \mathcal O }({10}^{5})$ with explicit resistivity η = 5 × 10⁻⁵. In Appendix D, we show the same analysis for the lower-resolution simulations, concluding that the extreme- and high-resolution results are in the plasmoid-mediated regime, whereas the standard- and low-resolution runs show reconnection rates a factor 2 to 3 larger than 0.01c and do not display plasmoids. The enhanced reconnection rate due to larger numerical resistivity at lower resolutions manifests itself as an increased flux decay rate and hence directly affects the flare timescale.

In Figure 5, we show the temperature (left column), plasma-β (middle column), and rest-mass density (right column) in both the meridional (top row) and equatorial (bottom row) plane for the miniflare in the extreme-resolution run at t ≈ 6800r_g c⁻¹. In this case, the accretion disk is not ejected far beyond 5r_g, still creating a spiral density gap causing the mass accretion rate to drop significantly. Reconnection occurs in a shorter, ≲5r_g plasmoid-unstable current sheet, very close to the horizon, and this is also the main area that is heated to relativistic temperatures T > 1. These miniflares could potentially result in smaller very-high-energy flares and shorter variability timescales (Abramowski et al. 2012).

Zoom In Zoom Out Reset image size

To probe the radiation emitted by accelerated particles in the reconnection layer, a self-consistent radiative kinetic approach is necessary (Hakobyan et al. 2019; Crinquand et al. 2020, 2021). Motivated by the results in the previous sections, we assume the flaring is associated with the formation of a transient macroscopic current sheet in a magnetospheric region near the event horizon without further relying on the GRMHD results We then use the well-constrained parameters for M87* and Sgr A* to estimate the expected emission properties due to reconnection occurring in the radiative regime.

4.1. M87* Flares Powered by Radiative Reconnection

In our simulations, the current sheet is fed by plasma in the jet at the floor density with a magnetization of ${\sigma }_{\max }=25$. In reality, the reconnection powering the flare close to the event horizon is fed by collisionless pair plasma from the jet with a rate of v_rec/c = 0.1 (Bransgrove et al. 2021) at magnetization ${\sigma }_{\mathrm{up}}={B}_{\mathrm{up}}^{2}/(4\pi {{nm}}_{{\rm{e}}}{c}^{2})=2{U}_{{\rm{B}}}/({{nm}}_{{\rm{e}}}{c}^{2})$, where n is the number density of electrons with mass m_e, B_up is the magnetic field strength upstream from the current sheet, and the magnetic energy density is ${U}_{{\rm{B}}}={B}_{\mathrm{up}}^{2}/8\pi$ ¹³ . The plasma particles are impulsively accelerated by nonideal electric fields at the X-points (Sironi & Spitkovsky 2014). When they encounter plasmoids, they experience strong synchrotron losses. To parameterize the effect of the radiation backreaction, we define the particle Lorentz factor ${\gamma }_{\mathrm{rad}}^{\mathrm{sync}}$ for which the radiation drag force is comparable to the force due to the accelerating electric field E ∼ B_up v_rec/c (Uzdensky et al. 2011):

$$\begin{eqnarray}&&2{\sigma }_{{\rm{T}}}{U}_{{\rm{B}}}{({\gamma }_{\mathrm{rad}}^{\mathrm{sync}})}^{2}={v}_{\mathrm{rec}}{{eB}}_{\mathrm{up}}/c,\end{eqnarray} \tag{ 2 }$$

where ${\sigma }_{{\rm{T}}}=(8\pi /3){r}_{{\rm{e}}}^{2}$ is the Thomson cross section, r_e = e²/(m_e c²) is the classical electron radius, and e is the electron's charge. We then find ${({\gamma }_{\mathrm{rad}}^{\mathrm{sync}})}^{2}=3{v}_{\mathrm{rec}}{B}_{\mathrm{cl}}/(2{{cB}}_{\mathrm{up}})$, where ${B}_{\mathrm{cl}}={m}_{{\rm{e}}}^{2}{c}^{4}/{e}^{3}\simeq 6\times {10}^{15}$ G is the classical magnetic field. The global magnetic field strength at 5r_g is estimated to be 1–30G (Event Horizon Telescope Collaboration et al. 2021), resulting in 5–150 G at the horizon, assuming a 1/r dependence (RBP20). We can compare this to the magnetic field strength in the jet, feeding the current sheet close to the event horizon of M87* by equating the observed limits on the total jet power L_jet ∼ 10⁴²–10⁴⁴ erg s⁻¹ (Prieto et al. 2016) to the Blandford–Znajek jet power ${L}_{\mathrm{BZ}}=\kappa {{\rm{\Omega }}}_{\mathrm{BH}}^{2}{\dot{\phi }}_{\mathrm{BH}}^{2}/(4\pi c)$, where κ ≈ 0.044 for a parabolic field geometry, Ω_BH = ac/2r_H ≃ c/2r_g is the black hole's angular frequency, M ≈ 6 · 10⁹ M_⊙ for M87*, and ${r}_{{\rm{H}}}={r}_{{\rm{g}}}(1+\sqrt{1-{a}^{2}})$ is the horizon radius (Blandford & Znajek 1977; Tchekhovskoy et al.2011), resulting in a range of B_horizon ∼ 20–200 G at the horizon. By normalizing to a fiducial B_up = 100 G in this range, we then obtain

$$\begin{eqnarray}&&{\gamma }_{\mathrm{rad}}^{\mathrm{sync}}\approx 3\cdot {10}^{6}{\left(\displaystyle \frac{{B}_{\mathrm{up}}}{100{\rm{G}}}\right)}^{-1/2}.\end{eqnarray} \tag{ 3 }$$

The magnetization σ_up sets the available magnetic energy per particle and determines the typical particle Lorentz factor, γ ∼ σ_up (which in GRMHD corresponds to the temperature of reconnection-heated fluid, whereas the bulk Lorentz factors of reconnection outflows scale as ${\rm{\Gamma }}\sim \sqrt{{\sigma }_{\mathrm{up}}}$), for the acceleration at X-points, if cooling were negligible (Sironi & Spitkovsky 2014; Guo et al. 2014; Werner et al. 2015). We can rewrite σ_up = ω_B/(2Ω_BH λ), where we plugged in the nominal electron gyrofrequency ω_B = eB_up/(m_e c) and defined the plasma density with respect to the Goldreich–Julian density, n = λ n_GJ = λΩ_BH B_up/(2π ce), where λ is the multiplicity of the pair cascade in the charge-starved gap in the funnel region λ ≲ 10³ (Chen & Yuan 2020; Crinquand et al. 2020) or of collisions of photons from the disk, if that process is more efficient (Moscibrodzka et al. 2011). The typical ratio between the electron gyrofrequency and the angular frequency of M87* is ω_B/Ω_BH ∼ 10¹⁴(M/6 · 10⁹ M_⊙)(B_up/100G), such that σ_up ∼ 10¹⁴(M/6 · 10⁹ M_⊙)(B_up/100G)/2λ. For these parameters, ${\gamma }_{\mathrm{rad}}^{\mathrm{sync}}\ll {\sigma }_{\mathrm{up}}$, so that leptons impulsively accelerated at X-points are quickly cooled in plasmoids (Hakobyan et al. 2019). Thus, the reconnection occurs in the radiative regime (Uzdensky 2011).

To understand the radiative efficiency of reconnection, we determine the magnetic compactness ℓ_B = U_B σ_T w/(m_e c²) (Beloborodov 2017). Using Equation (2) and the ω_B/Ω_BH relation, we can rewrite ${{\ell }}_{{\rm{B}}}={v}_{\mathrm{rec}}w{\omega }_{{\rm{B}}}/({c}^{2}{({\gamma }_{\mathrm{rad}}^{\mathrm{sync}})}^{2})$ and obtain

$$\begin{eqnarray}&&{{\ell }}_{{\rm{B}}}\sim 1\left(\displaystyle \frac{w}{1{r}_{{\rm{g}}}}\right)\left(\displaystyle \frac{M}{6\cdot {10}^{9}{M}_{\odot }}\right){\left(\displaystyle \frac{{B}_{\mathrm{up}}}{100{\rm{G}}}\right)}^{2},\end{eqnarray} \tag{ 4 }$$

so ℓ_B ∼ 1, suggesting potentially efficient pair production but negligible annihilation (Beloborodov 2017). In this regime, the cooling time of accelerated particles, ct_sync/w ∼ 1/(ℓ_B γ), is much shorter than the light-crossing time of the current sheet. Inverse Compton (IC) cooling of accelerated particles on the ∼10⁴¹ erg s⁻¹ low-energy photons with energy density ${U}_{\mathrm{rad}}^{\mathrm{soft}}\sim 0.003$ erg cm⁻³ in the inner 10r_g results in ${\gamma }_{\mathrm{rad}}^{\mathrm{IC}}\sim {\gamma }_{\mathrm{rad}}^{\mathrm{sync}}\sqrt{{U}_{{\rm{B}}}/{U}_{\mathrm{rad}}^{\mathrm{soft}}}\sim {10}^{9}$ (Broderick & Tchekhovskoy 2015; EHT MWL Science Working Group et al. 2021), which is well above ${\gamma }_{\mathrm{rad}}^{\mathrm{sync}}$. The jet's magnetic field reconnects with a rate of 0.1c in the collisionless radiative regime, after which all reconnected power is directly radiated such that the higher energy density of photons produced by accelerated particles, ${U}_{\mathrm{rad}}^{\mathrm{rec}}\sim 0.1{U}_{{\rm{B}}}$ and hence L_rad ∼ 0.1L_jet (Beloborodov 2017; Bransgrove et al. 2021), can lead to very efficient IC cooling. The exact result depends on the spectral shape and reduction by Klein–Nishina effects.

The peak of the synchrotron radiation spectrum is expected to be at the synchrotron burn-off limit ${{ \mathcal E }}_{\mathrm{ph}}\sim {({\gamma }_{\mathrm{rad}}^{\mathrm{sync}})}^{2}{\hslash }{\omega }_{{\rm{B}}}\sim 200\mathrm{MeV}$ (Uzdensky et al. 2011), which is independent of the magnetic field strength. The highest-energy photons will be produced by IC scattering. Conservatively, the characteristic photon energy that can be produced is $\max ({{ \mathcal E }}_{\mathrm{ph}})={m}_{{\rm{e}}}{c}^{2}{\gamma }_{\mathrm{rad}}^{\mathrm{sync}}\sim 0.511\mathrm{MeV}\cdot {\gamma }_{\mathrm{rad}}^{\mathrm{sync}}\sim$ a few TeV. Additionally, particles can be accelerated beyond $\gamma \gt {\gamma }_{\mathrm{rad}}^{\mathrm{sync}}$ because synchrotron cooling is suppressed in X-points (Uzdensky et al. 2011; Cerutti et al. 2014). For photons with energy above the electron rest-mass energy m_e c² = 0.5MeV, e^± pairs are created if there are enough photon–photon collisions with seed photons with low energy ${{ \mathcal E }}_{{\rm{s}}}\sim {({m}_{{\rm{e}}}{c}^{2})}^{2}/{{ \mathcal E }}_{\mathrm{ph}}$. High-energy photons of energy ${{ \mathcal E }}_{\mathrm{ph},\mathrm{TeV}}$ produced in the magnetospheric region around the current sheet will interact most efficiently with seed photons of energy ${{ \mathcal E }}_{{\rm{s}}}\sim (1\mathrm{TeV}/{{ \mathcal E }}_{\mathrm{ph},\mathrm{TeV}})$ eV. Given the uncertainties about the density of a 1 eV photon field near the event horizon during the flaring state, the escape of TeV photons from the region is an open question (Levinson & Rieger 2011; EHT MWL Science Working Group et al. 2021). Conservatively, if ∼1% of the reconnection-dissipated power U_rad ∼ 0.1U_B, L_rad ∼ 0.1L_jet ∼ 10⁴¹–10⁴³, is emitted in very-high-energy photons, a γ-ray flux of 10³⁹–10⁴¹ erg s⁻¹ can be emitted as a flare.

Our extreme-resolution GRMHD simulation shows transient flaring periods where the mass accretion rate (and thus the luminosity of seed photons) drops significantly, by a factor ∼5–10, resulting in large low-density regions, such that opacity constraints for the escape of γ-ray photons from the equatorial current sheet are less strict than during a quiescent state. The decrease in the mass accretion rate and the local soft photon field can also create favorable conditions for the activation of pair discharges on the jet's magnetic field lines and the potential escape of TeV emission from spark gaps, if the opacity becomes prohibitive during the quiescent state (Levinson & Rieger 2011; Crinquand et al. 2020). The flaring state is distinctively different from the quiescent state observed by Event Horizon Telescope Collaboration et al. (2019a), implying that observations during a mass-accretion-rate drop/flare may result in different 230 GHz images (K. Chatterjee et al., 2022, in preparation.). The magnetic flux decay and mass-accretion drop lasts for a period of ∼100r_g c⁻¹ ∼ 1 month for M87*, which is longer than the typical observed ∼1–3 day TeV flux rise and decay timescales (Abramowski et al. 2012). However, in a collisionless plasma, the magnetic flux decay period is typically ∼3–10 times shorter due to the faster reconnection rate of v_rec ≈ 0.1c (Bransgrove et al. 2021) compared to v_rec ≈ 0.01c in GRMHD models, ¹⁴ resulting in a flare timescale of ∼ few days. We find that pair production in the current sheet can efficiently mass load the jet with electrons and positrons with energies γ ∼ 1–1000, which can emit synchrotron photons with energies ranging from the radio to optical wavelengths (see Appendix A).

4.2. Sgr A* Flares Powered by Radiative Reconnection

By conducting extreme-resolution 3D GRMHD simulations, we have shown that, during periods of magnetic flux decay at the horizon, MAD flows form transient and nonaxisymmetric magnetospheres that possess special qualities revealed only at such high resolutions. Namely, these eruptions lead to a substantial, order-of-magnitude drop in the mass accretion rate and the formation of a thin equatorial current sheet that extends from the horizon out to ∼5–10 r_g into the disk and separates the two polar jets. This current sheet is filled with the electron–positron plasma from the jets and reconnects in the plasmoid-mediated regime. The formation of plasmoids is revealed here for the first time in 3D thanks to the unusually high resolutions achieved in this work, N_r × N_θ × N_ϕ = 5376 × 2304 × 2304. Reconnection-heated to relativistic temperatures, the plasma in the current sheet escapes the black hole's gravitational pull through the exhaust of the reconnection layer: this injects magnetic flux tubes filled with the low-density pair plasma into the accretion disk and hot plasma along the jet–disk boundary. This reconnection-heated plasma can produce a multiwavelength flare. Hot flux tubes orbit in the accretion disk and can remain coherent for one to a few orbital periods. The timescales of the flare are directly governed by the reconnection rate in the equatorial current sheet. We have shown that this rate decreases with increasing numerical resolution until the critical resolution, beyond which it reaches the universal converged value that no longer changes when the resolution is increased any further. Importantly, only at such high resolutions, the structure of the current sheet—X-points and plasmoids—are resolved for the first time with our extreme-resolution 3D GRMHD simulations.

The universal reconnection rate directly sets the magnetic flux decay rate at the horizon. Other studies have related flux decay at the horizon with flares (Ball et al. 2018; Dexter et al. 2020; Chashkina et al. 2021; Scepi et al. 2021) or observed orbiting flux tubes in retrograde disks (those rotating in the opposite direction to their black hole; Porth et al. 2021). However, due to limited numerical resolution, they did do not capture plasmoid-mediated reconnection as the power source and did not identify a direct link between the magnetic flux decay at the event horizon and its origin in reconnection in the equatorial magnetospheric current.

We note that the trigger behind such large flux eruption events is still not understood. Large flares occur when the accretion is governed by large, low azimuthal mode-number spiral RTI modes. It is as of yet unclear why the accretion state switches from a large spectrum of RTI modes in quiescence to a single azimuthal spiral RTI mode during the flare.

In reality, the reconnection powering the flare is fed by highly magnetized pair plasma that eventually ends up in the hot flux tube, buoyantly rising in and mixing with the electron-ion plasma that makes up the accretion disk. Additionally, matter originating from the jet that enters into the equatorial current sheet and gets heated by reconnection, can travel along the jet sheath for large distances, during and shortly after a flux eruption. The temperature of this reconnection-heated matter is proportional to the magnetization in the jet, which in GRMHD simulations is set by the density floor. Therefore, the temperature in the parts of the jet sheath that are causally connected to the exhaust of the reconnection layer cannot be used during a flare period to determine its emission properties. The main uncertainty is the electron temperature, which is unknown in GRMHD simulations. Commonly used parametrized relations connecting the temperatures of ions and electrons based on local plasma-β or σ values in the accretion flow (e.g., Moscibrodzka et al. 2016; Davelaar et al. 2019; Event Horizon Telescope Collaboration et al. 2019b; Dexter et al. 2020; Yoon et al. 2020; Chatterjee et al. 2021, and references therein) or two-temperature GRMHD approaches (Ressler et al. 2015; Chael et al. 2019) therefore cannot describe the nonthermal emission from these events, which involves reconnection in high-σ collisionless pair plasma regime, transport and cooling of nonthermal lepton distributions, as well as efficient pair production.

We note that while the reconnection rate in the equatorial current sheet is converged in GRMHD at the extremely high numerical resolutions used in this work, it converges to v_rec/v_A ∼ 0.01, which is 1 order of magnitude lower than the converged value of ∼0.1 in kinetic simulations (Bransgrove et al. 2021). This difference comes from GRMHD simulations being unaware of the collisionless plasma microphysics, which is important at scales where reconnection happens, i.e., electron skin depth. Incorporating nonideal effects beyond scalar resistivity (e.g., Ripperda et al. 2019b) into GRMHD simulations, such as electron inertia and anisotropic electron pressure tensor effects in Ohm's law (Most et al. 2021), holds promise of matching the (collisional) GRMHD and collisionless reconnection rates (Ng et al. 2020). General relativistic (radiative) kinetic simulations (e.g., Parfrey et al. 2019; Crinquand et al. 2020, 2021) are crucial for probing the nonthermal effects and the impact of the higher reconnection rate in collisionless plasma on the flare properties.

In upcoming work, we will investigate the radiative properties of the flares and the consequences for the image variability as observed by the Event Horizon Telescope (K. Chatterjee et al., 2022, in preparation.). During large flux eruptions, the accretion disk is ejected over a large fraction of the azimuthal angle. The very hot current sheet will then emit nonthermal emission powered by reconnection in the inner 10r_g, that may not be in the radio band. This may result in an observable dimming of the radio image, potentially distinguishable by the Event Horizon Telescope concurrent with multiwavelength flare observations for Sgr A*. Flares are less frequent for M87*, and hence observing potential dimming requires much longer Event Horizon Telescope observation windows, or several observations for separated periods.

In this work, we have for the first time reached a numerical Lundquist number above the plasmoid instability threshold for the largest current sheets close to the event horizon in MAD flows with 3D GRMHD simulations. The formation of these macroscopic plasmoid-unstable current sheets is similar in 2D resistive GRMHD simulations with a resolved explicit resistivity (RBP20). We robustly find that reconnection in the largest current sheets in MAD flows can act as the powering mechanism for large, bright, rapid flares originating from near the event horizon. MHD turbulence is known to intermittently form plasmoid-unstable current sheets at smaller scales (Zhdankin et al. 2013, 2017; Dong et al. 2018; Chernoglazov et al. 2021) that are not resolved in our simulations and that may substantially modify the turbulent cascade and dissipation at even higher Lundquist numbers S ≳ 10⁶ (Boldyrev & Loureiro 2017; Comisso et al. 2018), potentially heating the accretion disk. Additionally, the ideal GRMHD approach taken here does not describe the dynamics of nonideal electric fields and resistive dissipation and relies on a numerical resistivity that is only controlled by numerical resolution. To analyze the formation of nonideal electric fields and probe heating through turbulent reconnection in the accretion disk on smaller scales, even higher-resolution resistive GRMHD simulations are required (RBP20; Chernoglazov et al. 2021).

The robust formation of a macroscopic plasmoid-unstable current sheet close to the event horizon that can heat and accelerate plasma and eject flux tubes as low-density hot spots into an orbiting disk in our extreme-resolution GRMHD simulation suggests that flux eruptions powered by magnetic reconnection are a widespread phenomenon that can potentially explain observations of bright, rapid TeV flares from M87* and flaring hot spots from Sgr A*.

We would like to thank Ashley Bransgrove, Alexander Chernoglazov, Luca Comisso, Doosoo Yoon, Hayk Hakobyan, Gabriele Ponti, Benjamin Crinquand, Elias Most, Amir Levinson, and Yuri Levin for useful discussions. B.R. and M.L. contributed equally to this work. This research was enabled by support provided by grant No. NSF PHY-1125915 along with an INCITE program award PHY129, using resources from the Oak Ridge Leadership Computing Facility, Summit, which is a US Department of Energy office of Science User Facility supported under contract DE-AC05-00OR22725, as well as Calcul Quebec (http://www.calculquebec.ca) and Compute Canada (http://www.computecanada.ca). The computational resources and services used in this work were partially provided by facilities supported by the Scientific Computing Core at the Flatiron Institute, a division of the Simons Foundation. This research is part of the Frontera (Stanzione et al. 2020) computing project at the Texas Advanced Computing Center (LRAC-AST20008). Frontera is made possible by National Science Foundation award OAC-1818253. B.R. is supported by a Joint Princeton/Flatiron Postdoctoral Fellowship. M.L. was supported by John Harvard Distinguished Science Fellowship and ITC Fellowship. K.C. is supported by a black hole Initiative Fellowship at Harvard University, which is funded by grants from the Gordon and Betty Moore Foundation, John Templeton Foundation, and the black hole PIRE program (NSF grant OISE-1743747). The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the Moore or Templeton Foundations. G.M. is supported by a Netherlands Research School for Astronomy (NOVA), Virtual Institute of Accretion (VIA) postdoctoral fellowship. A.P. acknowledges support by the National Science Foundation under Grants No. AST-1910248 and PHY-2010145. Research at the Flatiron Institute is supported by the Simons Foundation. K.C. and S.M. are thankful for support by a Dutch Research Council (NWO) VICI award, grant No. 639.043.513. A.T. acknowledges support by Northwestern University and by the National Science Foundation grants AST-1815304, AST-1911080. Z.Y. is supported by a UK Research & Innovation (UKRI) Stephen Hawking Fellowship.

Pair production in the current sheet near M87* can significantly contribute to the mass loading of the jet. The optical depth for photons of energy ${{ \mathcal E }}_{\mathrm{ph}}$ is ${\tau }_{\gamma \gamma }\sim 0.1{\sigma }_{{\rm{T}}}{U}_{\mathrm{rad}}^{{\rm{s}}}w/{{ \mathcal E }}_{{\rm{s}}}$, where ${U}_{\mathrm{rad}}^{{\rm{s}}}$ is the energy density of photons at energy ${{ \mathcal E }}_{{\rm{s}}}\sim {({m}_{{\rm{e}}}{c}^{2})}^{2}/{{ \mathcal E }}_{\mathrm{ph}}$ such that the photon–photon pair production opacity is maximal for ${{ \mathcal E }}_{{\rm{s}}}{{ \mathcal E }}_{\mathrm{ph}}\sim {({m}_{{\rm{e}}}{c}^{2})}^{2}$. Most of the reconnected power, L_rad ∼ 0.1L_jet, is radiated around the burn-off limit, ∼200MeV, and the peak can be quite broad (Hakobyan et al. 2019). Estimating ${U}_{\mathrm{rad}}^{{\rm{s}}}\sim {L}_{\mathrm{rad}}/4\pi {w}^{2}c$, we get ${\tau }_{\gamma \gamma }\sim 0.1{\sigma }_{{\rm{T}}}{L}_{\mathrm{rad}}/(4\pi {cw}{{ \mathcal E }}_{{\rm{s}}})$. For w ∼ r_g and for typical photon energies ${{ \mathcal E }}_{\mathrm{ph}}\sim {{ \mathcal E }}_{{\rm{s}}}\sim \mathrm{MeV}$, we find ${\tau }_{\gamma \gamma }\sim {10}^{-2}{\sigma }_{{\rm{T}}}{L}_{\mathrm{jet}}/(4\pi {{cr}}_{{\rm{g}}}{{ \mathcal E }}_{{\rm{s}}})\sim {10}^{-4}$ for a jet power L_jet ∼ 10⁴³ erg s⁻¹. The rate of pair creation is then ${\dot{N}}_{\mathrm{pair}}\sim {\tau }_{\gamma \gamma }\cdot (0.1{L}_{\mathrm{jet}}/{{ \mathcal E }}_{\mathrm{ph}})\sim {10}^{44}{s}^{-1}$. We can compare this estimate to the Goldreich–Julian number flux ${\dot{N}}_{\mathrm{GJ}}={I}_{\mathrm{GJ}}/e$, where ${I}_{\mathrm{GJ}}\simeq {({{cL}}_{\mathrm{jet}})}^{1/2}$ is the Goldreich–Julian current, such that ${\dot{N}}_{\mathrm{GJ}}\sim {({{cL}}_{\mathrm{jet}})}^{1/2}/e\sim {10}^{36}{s}^{-1}$. The resulting multiplicity $\lambda ={\dot{N}}_{\mathrm{pair}}/{\dot{N}}_{\mathrm{GJ}}\sim {10}^{8}$ indicates that pairs produced in the current sheet with energies γ ∼ 1–1000 can significantly contribute to mass loading of the jet and emit synchrotron photons with energies ∼ℏ ω_B γ² ∼ 10⁴GHz(γ/400)², ranging from the radio to optical wavelengths.

We performed two additional 2D GRMHD simulations to show that reconnection heats the plasma to $T\sim {\sigma }_{\max }$ and ${\rm{\Gamma }}\sim \sqrt{{\sigma }_{\max }}$, for floor magnetizations in the jet of ${\sigma }_{\max }=25,60,100$. Figure 6 shows both the Lorentz factor γ (top row) and temperature T (bottom row) for the three values of ${\sigma }_{\max }$. The plasma in the current sheath is indeed heated to $T\approx {\sigma }_{\max }$ and ${\rm{\Gamma }}\approx \sqrt{{\sigma }_{\max }}$, and the reconnection exhaust deposits hot plasma in the jet sheath up to large distances.

Zoom In Zoom Out Reset image size

We analyze the effect of the resolution on the reconnection rate by showing the magnetic field components in minimum variance coordinates for the high- (left panels), standard- (middle panels), and low-resolution (right panels) runs in Figure 7, similar to Figure 4.

Zoom In Zoom Out Reset image size

One can see that the reconnection layer broadens for decreasing resolutions. This results in an increased reconnection rate, in accordance with an increasing numerical resistivity due to the decreasing resolution. In the high-resolution run, we detect plasmoids, and hence the numerical Lundquist number is still above the threshold S_crit = 10⁴, confirmed by the measured reconnection rate of v_rec ≈ 0.01c. For the standard-resolution run, the reconnection rate increases to v_rec ≈ 0.02c, and for the low-resolution run, v_rec ≈ 0.025c, so that the Lundquist number is of the order of $S={({v}_{\mathrm{rec}}/c)}^{-2}\sim { \mathcal O }({10}^{3})$. No plasmoids are detected in the low- and standard-resolution simulations. Note that the typical resolutions used in Event Horizon Telescope Collaboration et al. (2019b), Dexter et al. (2020), Porth et al. (2021), and Scepi et al. (2021) are below our standard resolution.

The increased reconnection rate enhances the dissipation of magnetic energy and directly affects the flux decay rate that governs the flare timescale. In Figure 8(B), we show this effect for the low- (green line) and standard-resolution (blue line) runs. For large flares, accompanied by a drop in the mass accretion rate (Figure 8(D)), the flux decays at a rate ∝ e^−t/350 (indicated by the dotted black lines, e.g., at t ≈ 11,300r_g c⁻¹, 15,500r_g c⁻¹, and 18,500r_g c⁻¹) instead of the converged ∝ e^−t/500. For miniflares, the current sheet extends less far from the event horizon and is naturally captured by cells of smaller volume due to the spherically logarithmic grid. Additionally, in ideal GRMHD simulations (i.e., relying on numerical resistivity), the reconnection rate is a function of the numerical resolution and hence of the radius due to the spherical grid. Reconnection close to the horizon, e.g., in a miniflare, will therefore occur closer to the asymptotic value of v_rec ∼ 0.01c than at larger radii. The miniflares are not accompanied by a significant drop in the mass accretion rate and show a decay rate e^−t/500 (indicated by black dashed lines) that is similar to the high- and extreme-resolution simulations. Miniflares occur more often at low resolutions than at high resolutions due to the larger numerical viscosity whereas large flares occur less frequently because the funnel region is not cleared out due to diffusion.

Zoom In Zoom Out Reset image size

The left panels show the magnetic flux (Figure 7(A)) and mass accretion rate (Figure 7(C)) for the initial 10, 000r_g c⁻¹ for all resolutions, showing that the simulations are in the quasi–steady state of MAD accretion after ∼5000r_g c⁻¹.

Figure 9 shows a large flare that forms during a magnetic flux decay at an enhanced rate ∝ e^−t/350, and a significant drop in the mass accretion rate, at low resolution. During the large flare, a very broad current sheet forms, indicated by the high plasma-β in panel (B). The magnetic field lines (green lines; left panels) diffuse through the equatorial sheet without reconnecting. Due to the applied floors and the large numerical diffusion at low resolution, the temperature in the jet region is particularly unreliable, showing a wavy pattern that is absent at high- and extreme-resolution. The area that is heated due to reconnection is broader due to numerical diffusion, and due to the large cells that increase in volume with the radius. This results in a large T > 1 area that lies fully in the equatorial plane (D) because of the thickness of the sheet in the z-direction. The heated area does not correspond to a clear floored region, which is visible in the middle and right plots of β and density ρ. After the flare (bottom two rows), the inner 10r_g is in the quiescent state, and an ejected low density flux tube with vertical field orbits at x ≈ −25r_g, y ≈ 0r_g

Zoom In Zoom Out Reset image size

During miniflares, a short and broad current sheet forms close to the horizon, but there is no clear expulsion of the accretion disk. Therefore, it is hard to distinguish the flare state from the quiescent state. The miniflare shows a clear magnetic flux decay at a rate ∝ e^−t/500 in accordance with a reconnection rate ∼0.01c. There is no clear drop in the mass accretion rate.